1 Introduction↩︎

Quantum simulation is a promising approach to continuum quantum dynamics, where a useful simulation must represent not only a Hilbert space but also the Hamiltonian that defines the physical problem [1]–[6]. In this setting, relativistic quantum simulation is a natural target because the kinetic energy is no longer simply quadratic in momentum. Even for a single positive-energy particle, the relativistic kinetic operator is a nonlinear function of the momentum, and therefore a lattice simulator must specify how the corresponding momentum moments are represented on a finite domain [7]–[11].

The relativistic kinetic term considered here is the positive-energy single-particle operator \[\begin{align} \hat{T}_{\rm rel} &=& mc^2 \left( \sqrt{\hat{\mathbb{1}} + \frac{\hat{p}^{2}}{m^2 c^2}} - \hat{\mathbb{1}} \right) = - m c^2 \sum_{k=1}^{\infty} \chi_k \, \hat{p}^{2k}, \label{eq:intro95positive95energy95operator} \end{align}\tag{1}\] for \(\chi_k = \frac{(-1)^{k} (2k-3)!!}{2^k k! \, \left( m\, c \right)^{2k}}\). In the weakly relativistic regime, has the expansion \[\begin{align} \hat{T}_{\rm rel} = \frac{\hat{p}^{2}}{2m} - \frac{\hat{p}^{4}}{8 m^3 c^2} + O\left( \frac{\hat{p}^{6}}{m^5 c^4} \right), \label{eq:intro95relativistic95expansion} \end{align}\tag{2}\] for continuum momentum operator \(\hat{p}\). Thus, a leading-order relativistic quantum simulation must estimate not only the non-relativistic momentum moment \(\langle\hat{p}^{2}\rangle\), but also \(\langle\hat{p}^{4}\rangle\) [12]. Therefore, the main question addressed in this work is how to define and estimate those moments consistently when the finite-domain relativistic Hamiltonian is subject to either periodic boundary conditions (PBC) or Dirichlet boundary conditions (DBC) [9], [13]–[15].

The boundary conditions are part of the definition of the finite-domain Hamiltonian, rather than a detail to be added after discretization. Under PBC, the spatial domain is closed into a circle, so the finite-difference link connecting the last grid point back to the first is physical. Under DBC, the wavefunction vanishes at the endpoints of an interval, and the physical finite-difference operator is an open-chain operator with no wrap-around link. A relativistic quantum simulator for PBC and DBC therefore corresponds to different lattice momentum operators and, consequently, to different moments. Thus, here we treat PBC and DBC as two boundary-condition-resolved relativistic quantum-simulation problems.

First-quantized grid encodings provide the common qubit-register language. For an \(L\)-qubit position register, a 1D wavefunction is represented as \(\left|\psi\right> = \sum_{j=0}^{2^L - 1} c_j\left|j\right>\), where the computational basis state \(\left|j\right>\) labels a grid point. In this representation, diagonal potentials are evaluated from the position-basis probabilities, while the derivative operators are expressed through finite-difference shift operators on the register [13]–[25]. The qubit-encoding issue is that the natural shift operation on an \(L\)-qubit register is cyclic. This is exactly the right unitary operation for PBC, but for DBC the same cyclic shift contains an unphysical wrap-around coupling. In the methodology developed below, this implementation issue is handled after the target PBC or DBC relativistic Hamiltonian has been fixed: the cyclic translation moments are used to evaluate the bulk finite-difference contribution, while DBC is enforced by boundary-local correction terms which can remove the wrap-around link.

With the above-described perspectives, we develop the first-quantized relativistic quantum-simulation framework under PBC and DBC. We first define the weakly relativistic lattice Hamiltonian for each boundary condition, with the fourth-order kinetic moment understood as the square of the same boundary-consistent momentum-squared operator. We then derive measurement estimators for the energy terms. The PBC case requires translation moments up to second order, and the DBC one uses the same translation moments together with a finite set of endpoint probabilities and endpoint or near-endpoint coherences. The diagonal potentials are estimated by position sampling and classical post-processing. In this way, the qubit-register treatment of boundary conditions supports our relativistic quantum simulation under PBC and DBC.

Our framework for relativistic quantum simulation is validated analytically and numerically. A periodic free particle tests the PBC lattice dispersion and its convergence to the continuum relativistic energy. A Dirichlet infinite square well tests the open-chain DBC spectrum and the boundary-local reconstruction of the kinetic moments that enter the non-relativistic and leading relativistic kinetic-energy terms. The smooth-potential ground states verify the full energy-estimation workflow with diagonal potentials, and finite-shot sampling tests confirm the expected inverse-square-root statistical scaling with the number of measurement shots. These benchmarks keep distinct the finite-grid discretization error, the weak-relativistic truncation error, and the statistical error of the estimator.

The contributions of this work are threefold. First, we formulate first-quantized weakly relativistic lattice Hamiltonians for finite-domain quantum simulation under PBC and DBC. Second, we derive observable estimators that reconstruct the second-order and fourth-order lattice momentum moments required by the weak-relativistic expansion from translation measurements and, for DBC, boundary-local overlap measurements. Third, we validate the construction with analytically solvable PBC and DBC systems, smooth-potential benchmarks, and finite-shot simulations. Together, these results establish a compact methodology for relativistic quantum simulation with PBC and DBC.

2 Relativistic Lattice Hamiltonians under Periodic and Dirichlet Boundary Conditions↩︎

We now define the lattice Hamiltonians that serve as the targets for relativistic simulations with PBC and DBC. The same first-quantized position register is used in both cases, but the finite-difference momentum operator is chosen to match the physical boundary condition. We denote the continuum momentum operator by \(\hat{p}\) and the corresponding lattice momentum operator by \(\hat{P}_{\tau}\) with \(\tau \in \{{\rm P}, {\rm D}\}\) specifying PBC or DBC.

2.1 First-quantized grids and weakly relativistic kinetic terms↩︎

We consider a single particle in one spatial dimension on a finite interval of physical length \(R\). An \(L\)-qubit position register represents \(N=2^L\) grid amplitudes, with a general lattice state written as \[\begin{align} \left|\psi\right> = \sum_{j=0}^{N-1} c_j \left|j\right>, \label{eq:lattice95state} \end{align}\tag{3}\] where \(\sum_{j=0}^{N-1} \left|c_j\right|^2 = 1\). Here, the computational basis state \(\left|j\right>\) labels a spatial grid point. The map from \(j\) to the physical coordinate is part of the finite-domain simulation specification and differs for PBC and DBC.

For PBC, the grid points are taken as, for \(j=0,\ldots,N-1\), \[\begin{align} x_j^{({\rm P})} = j\Delta_{\rm P}, \quad \Delta_{\rm P} = \frac{R}{N}. \label{eq:pbc95grid} \end{align}\tag{4}\] The last point is connected back to the first point, so that the finite difference is naturally represented by a cyclic translation on the register.

For DBC, we instead use \(N\) interior grid points, for \(j=0,\ldots,N-1\), \[\begin{align} x_j^{({\rm D})} = (j+1)\Delta_{\rm D}, \quad \Delta_{\rm D} = \frac{R}{N+1}, \label{eq:dbc95grid} \end{align}\tag{5}\] together with the virtual boundary values, \[\begin{align} \psi(0)=0, \quad \psi(R)=0 . \label{eq:dbc95virtual95boundaries} \end{align}\tag{6}\] This convention implements DBC as an open-chain finite-difference operator. The endpoints of the computational register, \(\left|0\right>\) and \(\left|N-1\right>\), therefore represent the first and last interior grid points, not the physical boundary points themselves. The corresponding physical-grid and dimensional conventions are summarized in 6.

In the weakly relativistic regime, where the relevant states have momentum support small compared with \(mc\), we use the expansion in Eq. (2 ). On the lattice, the corresponding boundary-dependent kinetic estimator through the leading relativistic correction is given by \[\begin{align} \hat{T}_{\tau} = \frac{\hat{P}^{2}_{\tau}}{2m} - \frac{\hat{P}^{4}_{\tau}}{8 m^3 c^2}, \quad \tau \in \{{\rm P}, {\rm D}\}, \label{eq:T95rel95lattice95tau} \end{align}\tag{7}\] where \(\tau={\rm P}\) and \(\tau={\rm D}\) denote PBC and DBC, respectively. Here and below, \[\begin{align} \hat{P}^{4}_{\tau} \equiv \left(\hat{P}^{2}_{\tau}\right)^2. \label{eq:P495definition} \end{align}\tag{8}\] Thus, the relativistic correction is built from the same boundary-consistent lattice momentum-squared operator used for the non-relativistic kinetic term.

For a scalar potential \(V(x)\), the lattice potential is diagonal in the position basis, \[\begin{align} \hat{V}_{\tau} = \sum_{j=0}^{N-1} V\bigl(x_j^{(\tau)}\bigr) \left|j\right>\!\!\left<j\right|. \label{eq:lattice95potential95sec2} \end{align}\tag{9}\] The corresponding total lattice Hamiltonian used for energy estimation is therefore \[\begin{align} \hat{H}^{tot}_{\tau} &=& \hat{T}_{\tau} + \hat{V}_{\tau}. \label{eq:lattice95hamiltonian95tau} \end{align}\tag{10}\]

2.2 Relativistic kinetic lattice Hamiltonian in PBC↩︎

For PBC, we introduce the unitary cyclic translation operator \(\hat{A}\) known as a quantum adder/subtractor [26], \[\begin{align} \hat{A}\left|j\right> &=& \left|j+1 \;{\rm mod} \;N\right>, \nonumber \\ \hat{A}^{\dagger}\left|j\right> &=& \left|j-1 \;{\rm mod} \;N\right>. \label{eq:cyclic95translation} \end{align}\tag{11}\] It satisfies \[\begin{align} \hat{A}^{\dagger}\hat{A} = \hat{A}\hat{A}^{\dagger} = \hat{\mathbb{1}}, \quad \hat{A}^{N} = \hat{\mathbb{1}}. \label{eq:cyclic95translation95unitarity} \end{align}\tag{12}\] The standard second-order finite-difference representation of the momentum-squared operator is given by a lattice Laplacian operator then \[\begin{align} \hat{P}^{2}_{\rm P} = - \frac{\hbar^2}{\Delta_{\rm P}^{2}} \left( \hat{A} + \hat{A}^{\dagger} - 2\hat{\mathbb{1}} \right). \label{eq:P295pbc95operator} \end{align}\tag{13}\] For a lattice state \(\left|\psi\right>\), this gives \[\begin{align} \langle \hat{P}^{2}_{\rm P} \rangle = \frac{2 \hbar^2}{\Delta_{\rm P}^{2}} \left( 1- \mathrm{Re} \langle\hat{A}\rangle \right), \label{eq:P295pbc95expectation} \end{align}\tag{14}\] where \(\mathrm{Re} \langle\hat{A}\rangle = \left( \langle \hat{A} \rangle + \langle \hat{A}^{\dagger} \rangle\right)/2\). The fourth-order momentum moment is obtained by squaring Eq. (13 ). Since \(\hat{A}\) is unitary, \[\begin{align} \hat{P}^{4}_{\rm P} = \frac{\hbar^4}{\Delta_{\rm P}^{4}} \left( \hat{A}^{2} + \hat{A}^{\dagger 2} - 4\hat{A} - 4\hat{A}^{\dagger} + 6\hat{\mathbb{1}} \right), \label{eq:P495pbc95operator} \end{align}\tag{15}\] and, therefore, \[\begin{align} \langle \hat{P}^{4}_{\rm P} \rangle = \frac{2 \hbar^4}{\Delta_{\rm P}^{4}} \left( \mathrm{Re} \langle\hat{A}^2\rangle - 4 \mathrm{Re} \langle\hat{A}\rangle + 3 \right). \label{eq:P495pbc95expectation} \end{align}\tag{16}\] Eqs (14 ) and (16 ) show that the leading relativistic correction under PBC requires only translation moments up to second order. In particular, \(\langle\hat{P}^{2}_{\rm P}\rangle\) depends on \(\langle\hat{A}\rangle\), while \(\langle\hat{P}^{4}_{\rm P}\rangle\) depends on both \(\langle\hat{A}\rangle\) and \(\langle\hat{A}^{2}\rangle\).

The PBC kinetic energy through order \(\hat{P}^{4}\) is thus \[\begin{align}\langle \hat{T}_{\rm P} \rangle = \frac{1}{2m} \langle \hat{P}^{2}_{\rm P} \rangle - \frac{1}{8 m^3 c^2} \langle \hat{P}^{4}_{\rm P} \rangle. \label{eq:T95pbc95expectation} \end{align}\tag{17}\] This expression is the PBC kinetic component of the relativistic simulation framework and will be used in the next Sec. 3 to define the translation-moment measurement protocol.

2.3 Relativistic kinetic lattice Hamiltonian in DBC↩︎

For DBC, the target finite-difference Hamiltonian is the open-chain Hamiltonian associated with the interior grid in Eq. (5 ). In a qubit-register implementation, it is convenient to express this open-chain operator using the same unitary cyclic translation \(\hat{A}\), supplemented by boundary-local terms that remove the wrap-around link. On the \(N\)-point interior grid, we define \[\begin{align} \hat{E}_0 = \left|N-1\right>\!\!\left<0\right| + \left|0\right>\!\!\left<N-1\right|. \label{eq:E095definition} \end{align}\tag{18}\] This operator \(\hat{E}_0\) identifies the coupling between the first and last interior grid points that would be present in a cyclic register translation but is absent from the DBC Hamiltonian.

The DBC momentum-squared operator is therefore \[\begin{align} \hat{P}^{2}_{\rm D} = -\frac{\hbar^2}{\Delta_{\rm D}^{2}} \left( \hat{A} + \hat{A}^{\dagger} - 2\hat{\mathbb{1}} - \hat{E}_0 \right). \label{eq:P295dbc95operator} \end{align}\tag{19}\] Equivalently, if \[\begin{align} \hat{P}^{2}_{\rm cyc}(\Delta_{\rm D}) = -\frac{\hbar^2}{\Delta_{\rm D}^{2}} \left( \hat{A} + \hat{A}^{\dagger} - 2\hat{\mathbb{1}} \right) \label{eq:P295cyc95dgrid} \end{align}\tag{20}\] denotes the cyclic finite-difference operator evaluated with the DBC grid spacing, then \[\begin{align} \hat{P}^{2}_{\rm D} = \hat{P}^{2}_{\rm cyc}(\Delta_{\rm D}) + \frac{\hbar^2}{\Delta_{\rm D}^{2}} \hat{E}_0. \label{eq:P295dbc95correction95operator} \end{align}\tag{21}\] Hence, \[\begin{align} \langle \hat{P}^{2}_{\rm D} \rangle = \langle \hat{P}^{2}_{\rm cyc}(\Delta_{\rm D}) \rangle + \frac{\hbar^2}{\Delta_{\rm D}^{2}} \langle \hat{E}_0 \rangle, \label{eq:P295dbc95expectation} \end{align}\tag{22}\] where \[\begin{align} \langle \hat{E}_0 \rangle = c_{N-1}^{\ast} c_0 + c_0^{\ast} c_{N-1}. \label{eq:E095expectation} \end{align}\tag{23}\]

The fourth-order DBC moment is defined as \[\begin{align} \hat{P}^{4}_{\rm D} = \left( \hat{P}^{2}_{\rm D} \right)^2. \label{eq:P495dbc95definition} \end{align}\tag{24}\] Using Eq. (19 ), one obtains \[\begin{align} \langle \hat{P}^{4}_{\rm D} \rangle &=& \langle \hat{P}^{4}_{\rm cyc}(\Delta_{\rm D}) \rangle + \frac{\hbar^4}{\Delta_{\rm D}^{4}} \left( 4\langle \hat{E}_0\rangle - \langle \hat{E}_1\rangle - \langle \hat{E}_2\rangle + \langle \hat{E}_0^2\rangle \right), \label{eq:P495dbc95expectation} \end{align}\tag{25}\] where \[\begin{align} \hat{P}^{4}_{\rm cyc}(\Delta_{\rm D}) = \frac{\hbar^4}{\Delta_{\rm D}^{4}} \left( \hat{A}^{2} + \hat{A}^{\dagger 2} - 4\hat{A} - 4\hat{A}^{\dagger} + 6\hat{\mathbb{1}} \right), \label{eq:P495cyc95dgrid} \end{align}\tag{26}\] and the additional boundary operators are \[\begin{align} \hat{E}_1 = \hat{A}\hat{E}_0 + \hat{E}_0\hat{A}^{\dagger}, \quad \hat{E}_2 = \hat{E}_0\hat{A} + \hat{A}^{\dagger}\hat{E}_0. \label{eq:E195E295definitions} \end{align}\tag{27}\] The derivation of Eq. (25 ) follows by expanding \((\hat{A}+\hat{A}^{\dagger}-2\hat{\mathbb{1}}-\hat{E}_0)^2\) and collecting the cyclic part and the boundary-local part. The full algebra is given in 7.

For later use, we record the explicit endpoint form of the boundary operators. For \(N>2\), \[\begin{align} \hat{E}_0^2 = \left|0\right>\!\!\left<0\right| + \left|N-1\right>\!\!\left<N-1\right|, \label{eq:E095squared95explicit} \end{align}\tag{28}\] and the compact definitions in Eq. (27 ) give \[\begin{align} \hat{E}_1 &=& 2\left|0\right>\!\!\left<0\right| + \left|1\right>\!\!\left<N-1\right| + \left|N-1\right>\!\!\left<1\right|, \nonumber \\ \hat{E}_2 &=& 2\left|N-1\right>\!\!\left<N-1\right| + \left|0\right>\!\!\left<N-2\right| + \left|N-2\right>\!\!\left<0\right|.~~~ \label{eq:E1295explicit} \end{align}\tag{29}\] Accordingly, \[\begin{align} \langle \hat{E}_0^2\rangle &=& \left|c_0\right|^2 + \left|c_{N-1}\right|^2, \nonumber \\ \langle \hat{E}_1\rangle &=& 2\left|c_0\right|^2 + c_1^{\ast} c_{N-1} + c_{N-1}^{\ast} c_1, \nonumber \\ \langle \hat{E}_2\rangle &=& 2\left|c_{N-1}\right|^2 + c_0^{\ast} c_{N-2} + c_{N-2}^{\ast} c_0. \label{eq:E01295expectation} \end{align}\tag{30}\] Thus, the DBC part of the implementation is boundary-local: it involves only endpoint probabilities and a small number of endpoint or near-endpoint coherences.

The DBC kinetic energy through order \(\hat{P}^{4}\) is \[\begin{align} \langle \hat{T}_{\rm D} \rangle = \frac{1}{2m} \langle \hat{P}^{2}_{\rm D} \rangle - \frac{1}{8m^3c^2} \langle \hat{P}^{4}_{\rm D} \rangle. \label{eq:T95dbc95expectation} \end{align}\tag{31}\] Combining Eqs. (22 ) and (25 ), this can be written as a cyclic translation-moment contribution plus a boundary correction: \[\begin{align}\langle \hat{T}_{\rm D} \rangle &=& \langle \hat{T}_{\rm cyc}(\Delta_{\rm D}) \rangle + \frac{\hbar^2}{2m\Delta_{\rm D}^{2}} \langle \hat{E}_0\rangle - \frac{\hbar^4}{8m^3c^2\Delta_{\rm D}^{4}} \left( 4\langle \hat{E}_0\rangle - \langle \hat{E}_1\rangle - \langle \hat{E}_2\rangle + \langle \hat{E}_0^2\rangle \right).~~~~~~~ \label{eq:T95dbc95correction95form} \end{align}\tag{32}\] This form completes the DBC Hamiltonian construction. The bulk contribution is evaluated by the same translation moments used for the PBC simulator, while the difference between DBC and a cyclic register translation is captured by a finite set of boundary-local expectation values.

3 Quantum Estimators for PBC and DBC Relativistic Simulation↩︎

We now convert the PBC and DBC relativistic lattice Hamiltonians of Sec. 2 into observable estimators. The aim is to specify the energy-evaluation layer of the quantum simulation, independent of any particular state-preparation ansatz. Once a first-quantized lattice state \(\left|\psi\right>\) as in Eq. (3 ) has been prepared, the weakly relativistic energy through \(\hat{P}^{4}\) is reconstructed from a small set of expectation values. PBC uses translation moments of the cyclic shift operator. DBC uses the same translation moments and adds the boundary-local overlap probabilities required by the open-chain Hamiltonian.

Throughout this section, an estimator obtained from finitely many shots will be denoted by a wide-tilde, for example \(\widetilde{m}_l\). The statistical analysis of the finite-shot fluctuations is given in 8.

3.1 PBC translation-moment estimator↩︎

For PBC, the kinetic moments in Eqs. (14 ) and (16 ) depend on the real parts of \(\langle\hat{A}\rangle\) and \(\langle\hat{A}^2\rangle\). Here, we define \[\begin{align} m_l = \mathrm{Re}\left<\psi\right|\hat{A}^{l}\left|\psi\right>, \quad l=1,2. \label{eq:ml95definition} \end{align}\tag{33}\] Since \(\langle \hat{A}^{l} + \hat{A}^{\dagger l} \rangle = 2\mathrm{Re} \langle\hat{A}^{l}\rangle = 2 m_l\), the PBC momentum moments can be written as \[\begin{align} \langle \hat{P}_{\rm P}^{2}\rangle &=& \frac{2\hbar^2}{\Delta_{\rm P}^{2}} \left( 1 - m_1 \right), \tag{34} \\ \langle \hat{P}_{\rm P}^{4}\rangle &=& \frac{\hbar^4}{\Delta_{\rm P}^{4}} \left( 2m_2 - 8m_1 + 6 \right). \tag{35} \end{align}\] Thus, the leading relativistic kinetic energy under PBC in Eq. (17 ) is reconstructed as \[\begin{align}\langle \hat{T}_{\rm P} \rangle &=& \alpha \left( 1 - m_1 \right) - \frac{\beta}{4} \left( m_2 - 4m_1 + 3 \right) = \left( \alpha - \frac{3\beta}{4} \right) + \left( \beta - \alpha \right) m_1 - \frac{ \beta}{4} m_2 . \label{eq:T95pbc95moments} \end{align}\tag{36}\] where \(\alpha = \hbar^2/(m \Delta^2_{\rm P})\) and \(\beta = \hbar^4/(m^3 c^2 \Delta^4_{\rm P})\).

The moments \(m_l\) are obtained from the standard controlled-unitary interferometric primitive [27], [28]. We state the identity explicitly because it is the basic observable estimator used throughout this work.

The quantum circuit of the translation-moment estimator is drawn in Fig. 1. Here, if the imaginary part of a translation moment is required, one may equivalently measure the ancilla in the \(Y\) basis, or insert the appropriate phase rotation before the final Hadamard measurement. The present kinetic estimators require only the real parts because \(\hat{P}^{2}\) and \(\hat{P}^{4}\) are Hermitian and depend on \(\hat{A}^{l}+\hat{A}^{\dagger l}\).

With \(M_l\) shots allocated to the \(l\)th translation moment, the empirical estimator is \[\begin{align} \widetilde{m}_l = \frac{1}{M_l} \sum_{r=1}^{M_l} z_r^{(l)}, \quad z_r^{(l)} \in \{-1,+1\}, \label{eq:empirical95translation95moment} \end{align}\tag{39}\] and Eqs. (34 )–(36 ) are evaluated by replacing \(m_l\) with \(\widetilde{m}_l\).

3.2 DBC boundary-overlap estimator↩︎

The DBC correction terms in Eqs. (22 ) and (25 ) are boundary-local. They contain endpoint probabilities and coherences between a small number of endpoint or near-endpoint basis states. These coherences have the generic form \(\left\langle \bigl( \left|f\right>\!\!\left<g\right| + \left|g\right>\!\!\left<f\right| \bigr) \right\rangle, \label{eq:generic95boundary95coherence}\) with computational basis states \(\left|f\right>\) and \(\left|g\right>\). The following identity gives a direct probability estimator for such terms.

The factor of two in Eq. (?? ) is essential: it follows from the normalisation of \(|s_{fg}^{+}\rangle=(\left|f\right>+\left|g\right>)/\sqrt{2}\). The identity holds for general complex amplitudes and returns the real coherence \(2\mathrm{Re}(c_f^{\ast}\, c_g)\), which is precisely the quantity needed for the Hermitian boundary operators. Schematic of the boundary-overlap estimator is given in Fig. 2. The corresponding phase-dependent identities are collected in 9.

For compactness, we define the pair-coherence estimator \[\begin{align} B_{f, g} \equiv 2P_{fg}^{+} - P_f - P_g. \label{eq:Bfg95definition} \end{align}\tag{42}\] Then, using the explicit DBC boundary operators in Eqs. (28 )–(29 ), their expectation values are \[\begin{align} b_0 &\equiv& \langle\hat{E}_0\rangle = B_{0,N-1}, \nonumber \\ b_{00} &\equiv& \langle\hat{E}_0^2\rangle = P_0 + P_{N-1}, \nonumber \\ b_1 &\equiv& \langle\hat{E}_1\rangle = 2P_0 + B_{1,N-1}, \nonumber \\ b_2 &\equiv& \langle \hat{E}_2\rangle = 2P_{N-1} + B_{0,N-2}. \label{eq:b95estimators} \end{align}\tag{43}\] Here, \(P_j=\left|\left<{j}|{\psi}\right>\right|^2\). We assume \(N>2\) so that the endpoint and near-endpoint labels appearing in Eq. (43 ) are distinct. The small grids \(N=1,2\) are degenerate as finite-difference representations of a Dirichlet interval and are not used in the numerical benchmarks.

The DBC momentum moments are therefore reconstructed from the same translation moments \(m_1, m_2\) and the four boundary quantities \(b_0, b_1, b_2, b_{00}\): \[\begin{align} \langle \hat{P}_{\rm D}^{2}\rangle &=& \frac{\hbar^2}{\Delta_{\rm D}^{2}} \left( 2 - 2 m_1 + b_0 \right), \tag{44} \\ \langle \hat{P}_{\rm D}^{4}\rangle &=& \frac{\hbar^4}{\Delta_{\rm D}^{4}} \left( 2m_2 - 8m_1 + 6 + 4b_0 - b_1 - b_2 + b_{00} \right). \tag{45} \end{align}\] The corresponding relativistic kinetic-energy estimator is \[\begin{align} \langle\hat{T}_{\rm D} \rangle = \frac{1}{2m} \langle\hat{P}_{\rm D}^{2}\rangle - \frac{1}{8 m^3 c^2} \langle\hat{P}_{\rm D}^{4}\rangle . \label{eq:T95dbc95estimator95moments} \end{align}\tag{46}\]

For reference, Table 1 lists the boundary terms needed through the leading relativistic correction. The table is written in terms of measurable probabilities rather than amplitudes, and therefore directly specifies the estimator inputs.

3.3 Potential-energy estimator and total-energy workflow↩︎

The potential term is simpler than the kinetic term because the lattice potential is diagonal in the computational basis. To compute the potential energy \(\langle \hat{V}_{\tau} \rangle\), we may use a quantum circuit introduced in [23], [24]. This scheme requires the intial preparation of a control qubit, wavefunction qubits and a diagonal density matrix \(\hat{\rho}_V= \sum_{j = 0}^{2^L-1} {\cal V}_j \left|j\right>\left<j\right|\), because \(\hat{V} = {\cal S}\, \hat{\rho}_V = {\cal S}\,\sum_{j=0}^{2^L-1} {\cal V}_j \left|j\right>\left<j\right|\) for the scale factor of the potential \({\cal S}\). Then, one performs a controlled block-SWAP gate between the wavefunction qubits and the density matrix.

However, we here examine a direct and practical method of position-sampling estimation for the potential energy. For the boundary conditions \(\tau \in \{{\rm P},{\rm D}\}\), Eq. (9 ) gives \[\begin{align} \langle\hat{V}_{\tau}\rangle = \sum_{j=0}^{N-1} V\bigl( x_j^{(\tau)} \bigr) \left|c_j\right|^2 . \label{eq:potential95position95average} \end{align}\tag{47}\] Thus, \(\langle\hat{V}_{\tau}\rangle\) is obtained by measuring the position register in the computational basis and applying a classical post-processing function \(j \mapsto V\bigl(x_j^{(\tau)}\bigr)\).

The total energy through the leading relativistic correction is then \[\begin{align} \langle\hat{H}_{\tau}^{tot}\rangle = \langle \hat{T}_{\tau} \rangle + \langle\hat{V}_{\tau}\rangle, \quad \tau \in \{{\rm P}, {\rm D}\}. \label{eq:total95energy95estimator95sec3} \end{align}\tag{49}\] The workflow is summarised in Algorithm 3. The same procedure can be used either as a post-processing estimator for a fixed state or as the energy-evaluation subroutine inside a variational optimisation loop [29]–[32]. In the latter case, the variational parameters determine the prepared state \(\left|\psi(\boldsymbol{\theta})\right>\), but the estimator identities themselves are unchanged.

The important structural point is that PBC and DBC fit into the same relativistic energy-estimation workflow. The translation moments encode the cyclic finite-difference contribution to the kinetic energy, while DBC adds only the boundary-local information needed to realize the open-chain Hamiltonian. Consequently, the leading relativistic correction does not require a qualitatively new measurement primitive beyond those already needed for first-quantized kinetic-energy estimation; it extends the measured observable set from \(\langle\hat{A}\rangle\) to \(\langle\hat{A}\rangle, \langle\hat{A}^{2}\rangle\) and, for DBC, a finite set of boundary-overlap terms.

4 Validation of the PBC/DBC Relativistic Simulation Framework↩︎

We now validate our PBC and DBC relativistic simulation framework developed in Secs. 2 and 3. The validation has two roles. First, analytically solvable systems test whether the PBC and DBC lattice Hamiltonians reproduce the expected spectra and continuum trends. Second, the estimator reconstruction is compared with direct matrix evaluation on the same lattice to verify the measurement formulas. This separation is important: the finite-difference lattice introduces a discretization error, the weakly relativistic expansion introduces a truncation error, and finite measurements introduce a statistical error. The matrix construction, estimator reconstruction from state vectors, benchmark parameters, and error metrics used for these validations are collected in 10.

Throughout the numerical examples, we use natural units, i.e., \(\hbar=m=1\), unless stated otherwise. The grid spacing is chosen according to the boundary condition, namely \(\Delta_{\rm P}=R/N\) for PBC and \(\Delta_{\rm D}=R/(N+1)\) for DBC. The relativistic energy reported below is the perturbative estimator \[\begin{align} \langle\hat{T}_{\tau}\rangle = \frac{1}{2m}\langle\hat{P}^{2}_{\tau}\rangle - \frac{1}{8 m^3 c^2}\langle\hat{P}^{4}_{\tau}\rangle, \quad \tau \in \{{\rm P}, {\rm D}\}. \label{eq:validation95T4} \end{align}\tag{50}\]

4.1 PBC benchmark: free-particle dispersion↩︎

For PBC, the natural analytic benchmark is the free particle. Define the lattice Fourier state: for \(n=0,1,\ldots,N-1\), \[\begin{align} \left|q_n\right> = \frac{1}{\sqrt{N}} \sum_{j=0}^{N-1} \exp\left(\frac{2\pi i n j}{N}\right) \left|j\right>. \label{eq:pbc95fourier95mode} \end{align}\tag{51}\] This state diagonalises the cyclic translation operator and therefore also the PBC finite-difference momentum operator.

For this benchmark, we compare three kinetic energies, \[\begin{align} T_{\rm cont} &=& mc^2\left[ \sqrt{1+\frac{p^{2}_{\rm cont}}{m^2c^2}} - 1 \right], \nonumber \\ T_{\rm lat} &=& mc^2\left[ \sqrt{1+\frac{p^{2}_{\rm lat}}{m^2c^2}} - 1 \right], \nonumber \\ T^{(4)}_{\rm pert} &=& \frac{p^{2}_{\rm lat}}{2m} - \frac{p^{4}_{\rm lat}}{8m^3c^2}. \label{eq:T95validation} \end{align}\tag{54}\] Fig. 4 shows the result for \(R=10\), \(c=2\), and \(n=1\), for which \(p_{\rm cont}^2/(m^2c^2) \simeq 9.87 \times 10^{-2}\). The square-root lattice energy \(T_{\rm lat}\) converges to the continuum value as \(L\) increases. The perturbative value \(T_{\rm pert}^{(4)}\) follows the same discretization trend but remains separated from \(T_{\rm cont}\) by the controlled truncation error of the weakly relativistic expansion.

4.2 DBC benchmark: infinite square well↩︎

For DBC, the corresponding analytic benchmark is the infinite square well on the interval \([0,R]\). With the interior-grid convention of Eq. (5 ), the normalised discrete sine modes are, for \(s=1,\ldots,N\), \[\begin{align} \psi_s(j) = \sqrt{\frac{2}{N+1}} \sin\left[ \frac{\pi s(j+1)}{N+1} \right]. \label{eq:dbc95sine95modes} \end{align}\tag{55}\] They diagonalise the open-chain finite-difference Laplacian.

The estimator-level validation is stronger than the spectral check alone. For each sine mode, we compute \(\langle\hat{P}^{2}_{\rm D}\rangle\) and \(\langle\hat{P}^{4}_{\rm D}\rangle\) in three independent ways: direct open-chain matrix evaluation, the boundary-corrected cyclic estimator in Eqs. (44 ) and (45 ), and the analytic eigenvalues above.

Fig. 5 shows the agreement for \(N=64\). The maximum relative discrepancy between the boundary-corrected estimator and the open-chain matrix evaluation is below \(3.1 \times 10^{-14}\) for \(\langle\hat{P}^{2}_{\rm D}\rangle\) and below \(7.4 \times 10^{-11}\) for \(\langle\hat{P}^{4}_{\rm D}\rangle\) over the displayed modes. This confirms that the DBC correction is not an approximation to the open-chain operator; it is an estimator decomposition of the same lattice operator.

4.3 Smooth-potential benchmark↩︎

The previous two tests validate the kinetic operators in exactly solvable settings. We next include smooth potentials to test the full energy workflow on nontrivial lattice ground states. For PBC, we use the periodic potential \[\begin{align} V_{\rm P}(x) = V_0\left[ 1 - \cos\left(\frac{2\pi x}{R}\right) \right], \label{eq:pbc95smooth95potential95validation} \end{align}\tag{57}\] while for DBC, we use a smooth confining potential in the box, \[\begin{align} V_{\rm D}(x) = \frac{1}{2}m\omega^2\left(x-\frac{R}{2}\right)^2. \label{eq:dbc95smooth95potential95validation} \end{align}\tag{58}\] For each boundary condition and each grid size, we diagonalise the non-relativistic lattice Hamiltonian \[\begin{align} \hat{H}_{\tau,{\rm nr}} = \frac{\hat{P}^{2}_{\tau}}{2m} + \hat{V}_{\tau} \label{eq:Hnr95validation} \end{align}\tag{59}\] and denote its ground state by \(\left|\psi_{0,\tau}\right>\). On this fixed state, we then evaluate \[\begin{align} E_{{\rm nr},\tau} &=& \left<\psi_{0,\tau}\right| \hat{H}_{\tau,{\rm nr}} \left|\psi_{0,\tau}\right>, \nonumber \\ \Delta E_{{\rm rel},\tau} &=& - \frac{1}{8m^3c^2} \left<\psi_{0,\tau}\right| \hat{P}^{4}_{\tau} \left|\psi_{0,\tau}\right>, \nonumber \\ E_{{\rm rel},\tau}^{(4)} &=& E_{{\rm nr},\tau} + \Delta E_{{\rm rel},\tau}. \label{eq:E95validations} \end{align}\tag{60}\] Each quantity is evaluated both by direct matrix multiplication and by the estimator reconstruction of Sec. 3.

Fig. 6 shows the results for \(R=10\), \(c=2\), \(V_0=0.5\), and \(\omega=0.4\). The grid sizes range from \(L=4\) to \(L=10\). The lines show direct matrix evaluation and the markers show the estimator reconstruction. The two are visually indistinguishable on the scale of the plot; the largest absolute difference in \(E_{\rm rel}^{(4)}\) over all displayed PBC and DBC data points is \(3.9 \times 10^{-7}\). The maximum value of the diagnostic ratio \(\langle\hat{P}^{2}_{\tau}\rangle/(m^2 c^2)\) is approximately \(5.1 \times 10^{-2}\), so the benchmark remains in the weakly relativistic regime for which Eq. (50 ) is intended.

4.4 Estimator-level finite-shot validation↩︎

Finally, we test the statistical scaling of the measurement primitives. The translation-moment measurement produces a binary outcome \(z^{(l)} \in \{-1,+1\}\) with mean \(m_l={\rm Re}\langle\hat{A}^{l}\rangle\), whereas each boundary-overlap component is estimated from Bernoulli projection measurements. Therefore, for a fixed state, the measurement contribution to the energy error is expected to scale as \[\begin{align} \epsilon_{\rm meas}(M) = O\bigl(M^{-1/2}\bigr), \label{eq:shot95noise95scaling95expected} \end{align}\tag{61}\] where \(M\) is the number of shots allocated to each measurement primitive.

We quantify this scaling using the root-mean-square error \[\begin{align} {\rm RMSE}(M) = \sqrt{ \mathbb{E}\left[ \left( \widetilde{T}_{\tau}(M) - T_{\tau} \right)^2 \right]} \label{eq:rmse95definition95validation} \end{align}\tag{62}\] for the kinetic estimator. The expectation is estimated by Monte Carlo sampling of the corresponding binary or Bernoulli measurement outcomes. The states are the smooth-potential ground states at \(L=5\), using the same physical parameters as in Sec. 4.3. For each value of \(M\), we use \(700\) independent Monte Carlo repetitions.

Fig. 7 shows the resulting RMSE curves. The fitted slopes are \(-0.51\) for the PBC kinetic estimator and \(-0.50\) for the DBC kinetic estimator, in agreement with the expected \(M^{-1/2}\) behaviour. The DBC curve has a slightly larger prefactor because it includes additional boundary-overlap probabilities, but it follows the same asymptotic statistical scaling as the translation-moment estimator.

Table 2 summarises the observable set used in the validation. The table emphasises that the leading relativistic correction enlarges the measured translation moments from \(l=1\) to \(l=1,2\) while DBC adds only boundary-local probability estimators.

5 Discussion↩︎

We have presented a first-quantized methodology for relativistic quantum simulation under periodic and Dirichlet boundary conditions, denoted here by PBC and DBC. The main object of the framework was the finite-domain relativistic Hamiltonian. Starting from the positive-energy kinetic operator, the weakly relativistic simulator required the boundary-consistent moments \(\langle\hat{P}^{2}\rangle\) and \(\langle\hat{P}^{4}\rangle\). PBC and DBC were handled within the same position-register encoding, but they corresponded to different lattice momentum operators and therefore to different estimator formulas.

For PBC, the finite-domain structure was naturally represented by the unitary cyclic translation on the qubit register. The non-relativistic kinetic term was reconstructed from \(\mathrm{Re}\langle\hat{A}\rangle\), and the leading relativistic correction additionally required \(\mathrm{Re}\langle\hat{A}^{2}\rangle\). For DBC, the target Hamiltonian was the open-chain finite-difference Hamiltonian on the interior grid. We implemented this target by decomposing it into a cyclic translation contribution plus boundary-local terms that cancel the unphysical wrap-around coupling. The additional DBC information was therefore limited to endpoint probabilities and endpoint or near-endpoint coherences associated with the expectation values of \(\hat{E}_0\), \(\hat{E}_1\), \(\hat{E}_2\), and \(\hat{E}_0^2\), rather than an extensive set of measurements over the whole grid.

The validation results support the interpretation. The PBC free-particle benchmark reproduced the lattice dispersion relation and its convergence toward the continuum relativistic energy. The DBC infinite-square-well benchmark showed that the boundary-local reconstruction gives the same \(\hat{P}_{\rm D}^{2}\) and \(\hat{P}_{\rm D}^{4}\) moments as direct open-chain matrix evaluation and as the analytic sine spectrum. The smooth-potential benchmarks confirmed that the full energy workflow, including diagonal potential sampling, agrees with direct matrix evaluation for nontrivial ground states under both boundary conditions. Finally, the finite-shot tests showed the expected \(M^{-1/2}\) scaling for the kinetic-energy estimator, with the DBC case carrying a larger but controlled prefactor because of the additional boundary-overlap measurements.

Our formulation is intended for weakly relativistic regimes in which the positive-energy square-root kinetic operator is accurately represented by the low-momentum expansion through \(\hat{P}^{4}\). In practical applications, the prepared lattice state should therefore have momentum support small compared with \(mc\), and cutoff-scale grid modes should not dominate the energy estimate. This is a physical limitation of the perturbative relativistic approximation, not of the estimator identities themselves. The identities remain exact for the lattice operators defined in Sec. 2; the expansion determines when the reconstructed quantity is a quantitatively accurate relativistic energy. Higher-order relativistic corrections and higher-order finite-difference stencils are summarized in 11.

We should also highlight the importance of treating the boundary conditions as an essential part of quantum simulation. This perspective is particularly relevant for finite-domain systems, where the imposed boundary can directly shape the physical observables being simulated. Our framework showed that such the boundary awareness can be incorporated within a first-quantized quantum-simulation setting in a systematic and practical way. Beyond the specific examples considered here, the same viewpoint may be useful for more complex potentials, higher-dimensional geometries, and refined relativistic approximations. We therefore expect this boundary-aware formulations to play an important role in developing reliable relativistic quantum simulations of continuum systems.

Acknowledgments↩︎

This work was supported by the Ministry of Science, ICT and Future Planning (MSIP) by the National Research Foundation of Korea (RS-2024-00432214, RS-2025-03532992, RS-2023-00281456, and RS-2023-NR119931) and the Institute of Information and Communications Technology Planning and Evaluation grant funded by the Korean government (RS-2019-II190003, “Research and Development of Core Technologies for Programming, Running, Implementing and Validating of Fault-Tolerant Quantum Computing System”). This work is also supported by the Grant No. K25L5M2C2 at the Korea Institute of Science and Technology Information (KISTI).

Data availability statement↩︎

The data that support the findings of this study are openly available at the following URL/DOI:

6 Grid conventions and dimensional analysis↩︎

This appendix fixes the relation between the physical coordinate, the dimensionless coordinate, and the lattice parameters used in the main text. We include these details because the finite-difference momentum operators contain explicit powers of the grid spacing, and a dimensionally consistent convention is essential when relativistic parameters are introduced.

Let \(x_{\rm phys} \in [0,R]\) denote the physical coordinate and let \(y = {x_{\rm phys}}/{R}\) be the corresponding dimensionless coordinate. Then, we have \[\begin{align} && \frac{d}{d x_{\rm phys}} = \frac{1}{R}\frac{d}{dy}, \nonumber \\ && \hat{p} = -i\hbar\frac{d}{d x_{\rm phys}} = -i\frac{\hbar}{R}\frac{d}{dy} . \label{eq:app95derivative95rescaling} \end{align}\tag{63}\] Thus, a finite-difference expression written in the dimensionless variable must carry the physical factor \(R^{-1}\) in the momentum and \(R^{-2}\) in the momentum squared.

For PBC, we use \(N=2^L\) grid points on the circle, \[\begin{align} y_j^{({\rm P})} &=& \frac{j}{N}, \nonumber \\ x_j^{({\rm P})} &=& R y_j^{({\rm P})} = \frac{j R}{N}, \nonumber \\ \Delta_{\rm P} &=& \frac{R}{N}. \label{eq:app95pbc95grid95physical} \end{align}\tag{64}\] For DBC, we use \(N\) interior grid points in the interval, \[\begin{align} y_j^{({\rm D})} &=& \frac{j+1}{N+1}, \nonumber \\ x_j^{({\rm D})} &=& R y_j^{({\rm D})} = \frac{(j+1)R}{N+1}, \nonumber \\ \Delta_{\rm D} &=& \frac{R}{N+1}. \label{eq:app95dbc95grid95physical} \end{align}\tag{65}\] The virtual boundary values are located at \(y=0\) and \(y=1\), or equivalently, at \(x_{\rm phys}=0\) and \(x_{\rm phys}=R\), and satisfy \[\begin{align} \psi(0)=0, \quad \psi(R)=0. \label{eq:app95dbc95boundary95values} \end{align}\tag{66}\] The computational states \(\left|0\right>\) and \(\left|N-1\right>\) in the DBC register therefore correspond to the first and last interior points, not to the physical boundary points.

The physical second-derivative finite difference is \[\begin{align} \hbar^2\frac{d^2}{d x_{\rm phys}^2} ~\longrightarrow~ \frac{\hbar^2}{\Delta_\tau^2}, \quad \tau \in \{{\rm P},{\rm D}\}. \label{eq:app95physical95finite95difference} \end{align}\tag{67}\] Equivalently, if a dimensionless grid spacing \(\delta_\tau = \frac{\Delta_\tau}{R}\) is used, then \[\begin{align} \frac{\hbar^2}{\Delta_\tau^2} = \frac{\hbar^2}{R^2\delta_\tau^2}. \label{eq:app95delta95conversion} \end{align}\tag{68}\] This is the point at which dimensional consistency is most easily lost: a grid spacing in the dimensionless coordinate must not be used as if it were a physical length.

The dimensionless parameter controlling the size of the lattice momentum scale relative to \(mc\) is \[\begin{align} \mu_\tau = \frac{\lambda_m}{\Delta_\tau} = \frac{\hbar}{mc\Delta_\tau}, \quad \tau \in \{{\rm P},{\rm D}\}, \label{eq:app95mu95tau95definition} \end{align}\tag{69}\] where \(\lambda_m\) is the reduced Compton wavelength: \(\lambda_m = \hbar /(mc)\).

Then, we have \[\begin{align} \frac{\hbar^2}{\Delta_\tau^2} &=& m^2 c^2 \mu_\tau^2. \label{eq:app95mu95tau95powers} \end{align}\tag{70}\] If one instead works with a dimensionless coordinate \(y\) and defines a notation, such as, \(L_m=\lambda_m/\delta_\tau\), then \(\lambda_m\) in that expression must itself be dimensionless, namely, \[\begin{align} \lambda_m^{({\rm dimless})} = \frac{\lambda_m}{R}. \label{eq:app95dimensionless95lambda} \end{align}\tag{71}\] With this convention, \[\begin{align} L_{m,\tau} = \frac{\lambda_m^{({\rm dimless})}}{\delta_\tau} = \frac{\lambda_m/R}{\Delta_\tau/R} = \frac{\lambda_m}{\Delta_\tau} = \mu_\tau. \label{eq:app95Lm95conversion} \end{align}\tag{72}\] Thus, the dimensionless lattice parameter that appears in the momentum operator is the ratio of the physical Compton wavelength to the physical grid spacing.

The positive-energy relativistic kinetic operator can be written in terms of \[\begin{align} \hat{Z}_\tau = \frac{\hat{P}_\tau^2}{m^2 c^2} \label{eq:app95Z95tau95definition} \end{align}\tag{73}\] as \[\begin{align} \hat{T}_{\rm rel, \tau} = mc^2\left[ (\hat{\mathbb{1}} + \hat{Z}_\tau)^{1/2} - \hat{\mathbb{1}} \right]. \label{eq:app95exact95lattice95relativistic95operator} \end{align}\tag{74}\] The expansion used in the main text is \[\begin{align} \sqrt{1+z} - 1 = \frac{z}{2} - \frac{z^2}{8} + \frac{z^3}{16} + O(z^4), \label{eq:app95scalar95rel95expansion} \end{align}\tag{75}\] and hence, \[\begin{align} \hat{T}_\tau = \frac{\hat{P}_\tau^2}{2m} - \frac{\hat{P}_\tau^4}{8 m^3 c^2}. \label{eq:app95T495tau95again} \end{align}\tag{76}\] For \(z \geq 0\), Taylor’s theorem gives \[\begin{align} 0 \leq (1+z)^{1/2} - 1 - \frac{z}{2} + \frac{z^2}{8} \leq \frac{z^3}{16}. \label{eq:app95remainder95scalar95bound} \end{align}\tag{77}\] Therefore, by functional calculus for the positive semidefinite operator \(\hat{Z}_\tau\), \[\begin{align} 0 \leq \hat{T}_{\rm rel,\tau}-\hat{T}_\tau \leq \frac{1}{16 m^5 c^4}\hat{P}_\tau^6, \label{eq:app95operator95remainder95bound} \end{align}\tag{78}\] where \(\hat{P}_\tau^6 \equiv (\hat{P}_\tau^2)^3\). Consequently, for any normalised lattice state \(\left|\psi\right>\), \[\begin{align} 0 \leq \left<\psi\right|\hat{T}_{\rm rel,\tau}\left|\psi\right> - \left<\psi\right|\hat{T}_\tau\left|\psi\right> \leq \frac{\left<\psi\right|\hat{P}_\tau^6\left|\psi\right>}{16 m^5 c^4}. \label{eq:app95state95remainder95bound} \end{align}\tag{79}\] This bound is not used as a numerical error estimate in the main text, but it makes explicit the low-momentum character of the perturbative approximation.

7 Derivation of the PBC and DBC momentum moments↩︎

Here we derive the momentum-moment identities used in Secs. 2 and 3. We keep the derivation algebraic, because the same manipulations extend directly to higher-order moments.

7.1 PBC translation moments↩︎

Let \(\hat{K} = \hat{A}+\hat{A}^{\dagger}-2\hat{\mathbb{1}}\), where \(\hat{A}\) is the cyclic translation satisfying \(\hat{A}^{\dagger}=\hat{A}^{-1}\) and \(\hat{A}^{N}=\hat{\mathbb{1}}\). The PBC momentum-squared operator is \[\begin{align} \hat{P}_{\rm P}^{2} = -\frac{\hbar^2}{\Delta_{\rm P}^{2}}\hat{K}. \label{eq:app95P2P95K} \end{align}\tag{80}\] Thus, for any positive integer \(r\), \[\begin{align} \hat{P}_{\rm P}^{2r} \equiv (\hat{P}_{\rm P}^{2})^r = \left( -\frac{\hbar^2}{\Delta_{\rm P}^{2}} \right)^r \hat{K}^r. \label{eq:app95general95PBC95moment95K} \end{align}\tag{81}\] Since \(\hat{A}\) and \(\hat{A}^{\dagger}\) commute, the binomial expansion gives \[\begin{align} \hat{K}^r &=& \sum_{q=0}^{r} \binom{r}{q} (-2)^{r-q} (\hat{A}+\hat{A}^{\dagger})^q \nonumber \\ &=& \sum_{q=0}^{r} \binom{r}{q} (-2)^{r-q} \sum_{s=0}^{q} \binom{q}{s} \hat{A}^{q-2s}. \label{eq:app95general95K95expansion} \end{align}\tag{82}\] Here, the negative powers are interpreted as powers of \(\hat{A}^{\dagger}\), i.e., \(\hat{A}^{-\ell}=\hat{A}^{\dagger\ell}\). Taking the expectation in a state \(\left|\psi\right>\) therefore reduces every PBC moment to translation moments \(\langle \hat{A}^{\ell}\rangle\).

For \(r=1\), \[\begin{align} \hat{P}_{\rm P}^{2} = -\frac{\hbar^2}{\Delta_{\rm P}^{2}} (\hat{A}+\hat{A}^{\dagger}-2\hat{\mathbb{1}}), \label{eq:app95PBC95P295explicit} \end{align}\tag{83}\] which is Eq. (13 ). For \(r=2\), \[\begin{align} \hat{K}^2 &=& (\hat{A}+\hat{A}^{\dagger}-2\hat{\mathbb{1}})^2 = \hat{A}^2 + \hat{A}^{\dagger 2} - 4\hat{A} - 4\hat{A}^{\dagger} + 6\hat{\mathbb{1}}, \label{eq:app95K95squared} \end{align}\tag{84}\] and therefore \[\begin{align} \hat{P}_{\rm P}^{4} = \frac{\hbar^4}{\Delta_{\rm P}^{4}} \left( \hat{A}^2 + \hat{A}^{\dagger 2} - 4\hat{A} - 4\hat{A}^{\dagger} + 6\hat{\mathbb{1}} \right), \label{eq:app95PBC95P495explicit} \end{align}\tag{85}\] which is Eq. (15 ).

It is useful to introduce \[\begin{align} m_\ell = {\rm Re}\left<\psi\right|\hat{A}^{\ell}\left|\psi\right>. \label{eq:app95m95ell95definition} \end{align}\tag{86}\] Then, \[\begin{align} \langle\hat{P}_{\rm P}^{2}\rangle &=& \frac{2\hbar^2}{\Delta_{\rm P}^{2}}(1 - m_1), \nonumber \\ \langle\hat{P}_{\rm P}^{4}\rangle &=& \frac{\hbar^4}{\Delta_{\rm P}^{4}}(2 m_2 - 8m_1 + 6), \label{eq:app95PBC95m195m2} \end{align}\tag{87}\] as used in Sec. 3.

7.2 DBC as an open-chain correction to the cyclic shift↩︎

Define the non-unitary open-chain forward shift \[\begin{align} \hat{B} = \sum_{j=0}^{N-2}\left|j+1\right>\!\!\left<j\right|. \label{eq:app95B95open95shift} \end{align}\tag{88}\] It shifts all interior grid points forward except that it does not connect \(\left|N-1\right>\) back to \(\left|0\right>\). The cyclic shift can be decomposed as \[\begin{align} \hat{A} &=& \hat{B}+\left|0\right>\!\!\left<N-1\right|, \nonumber \\ \hat{A}^{\dagger} &=& \hat{B}^{\dagger}+\left|N-1\right>\!\!\left<0\right|. \label{eq:app95A95B95decomposition} \end{align}\tag{89}\] With \(\hat{E}_0 = \left|N-1\right>\!\!\left<0\right|+\left|0\right>\!\!\left<N-1\right|\), we have \[\begin{align} \hat{A}+\hat{A}^{\dagger} - \hat{E}_0 = \hat{B} + \hat{B}^{\dagger}. \label{eq:app95A95minus95E095open95shift} \end{align}\tag{90}\] Therefore, we verify \[\begin{align} \hat{A} + \hat{A}^{\dagger} - 2\hat{\mathbb{1}} - \hat{E}_0 = \hat{B} + \hat{B}^{\dagger} - 2\hat{\mathbb{1}}, \label{eq:app95open95laplacian95identity} \end{align}\tag{91}\] which is exactly the open-chain second-difference operator. The DBC momentum operator is thus \[\begin{align} \hat{P}_{\rm D}^{2} &=& -\frac{\hbar^2}{\Delta_{\rm D}^{2}} (\hat{B}+\hat{B}^{\dagger}-2\hat{\mathbb{1}}) = -\frac{\hbar^2}{\Delta_{\rm D}^{2}}(\hat{K}-\hat{E}_0), \label{eq:app95P2D95open95chain95identity} \end{align}\tag{92}\] where \(\hat{K}=\hat{A}+\hat{A}^{\dagger}-2\hat{\mathbb{1}}\). Equivalently, \[\begin{align} \hat{P}_{\rm D}^{2} = \hat{P}_{\rm cyc}^{2}(\Delta_{\rm D}) + \frac{\hbar^2}{\Delta_{\rm D}^{2}}\hat{E}_0. \label{eq:app95P2D95cyclic95plus95E0} \end{align}\tag{93}\] By taking the expectation values, we have Eq. (22 ).

7.3 Full derivation of the DBC fourth moment↩︎

The fourth DBC moment is the square of \(\hat{P}_{\rm D}^{2}\): \[\begin{align} \hat{P}_{\rm D}^{4} = \left( \hat{P}_{\rm D}^{2}\right)^2 = \left(\frac{\hbar^2}{\Delta_{\rm D}^{2}}\right)^2(\hat{K} - \hat{E}_0)^2. \label{eq:app95P4D95start} \end{align}\tag{94}\] Expanding the non-commuting product gives \[\begin{align} (\hat{K}-\hat{E}_0)^2 = \hat{K}^2 - \hat{K}\hat{E}_0 - \hat{E}_0\hat{K} + \hat{E}_0^2. \label{eq:app95Q95squared95expansion} \end{align}\tag{95}\] The first term is the cyclic fourth-moment operator evaluated with the DBC grid spacing. The mixed terms are \[\begin{align} \hat{K}\hat{E}_0 + \hat{E}_0\hat{K} &=& (\hat{A}+\hat{A}^{\dagger}-2\hat{\mathbb{1}})\hat{E}_0 + \hat{E}_0(\hat{A}+\hat{A}^{\dagger}-2\hat{\mathbb{1}}) \nonumber\\ &=& \hat{A}\hat{E}_0 + \hat{A}^{\dagger}\hat{E}_0 + \hat{E}_0\hat{A} + \hat{E}_0\hat{A}^{\dagger} - 4\hat{E}_0. \label{eq:app95mixed95terms95expanded} \end{align}\tag{96}\] Using \(\hat{E}_1 = \hat{A}\hat{E}_0 + \hat{E}_0\hat{A}^{\dagger}\) and \(\hat{E}_2 = \hat{E}_0\hat{A} + \hat{A}^{\dagger}\hat{E}_0\), we find the following: \[\begin{align} \hat{K}\hat{E}_0 + \hat{E}_0\hat{K} = \hat{E}_1 + \hat{E}_2 - 4\hat{E}_0. \label{eq:app95mixed95terms95E195E2} \end{align}\tag{97}\] Therefore, we have \[\begin{align} \hat{P}_{\rm D}^{4} = \hat{P}_{\rm cyc}^{4}(\Delta_{\rm D}) + \frac{\hbar^4}{\Delta_{\rm D}^{4}} \left( 4\hat{E}_0 - \hat{E}_1 - \hat{E}_2 + \hat{E}_0^2 \right), \label{eq:app95P4D95full95operator} \end{align}\tag{98}\] and taking expectation values gives Eq. (25 ).

For \(N>2\), the explicit endpoint forms follow directly. First, we have \[\begin{align} \hat{E}_0^2 &=& \bigl( \left|N-1\right>\!\!\left<0\right| + \left|0\right>\!\!\left<N-1\right| \bigr)^2 = \left|N-1\right>\!\!\left<N-1\right| + \left|0\right>\!\!\left<0\right|. \label{eq:app95E095square95derivation} \end{align}\tag{99}\] Next, \[\begin{align} \hat{A}\hat{E}_0 &=& \left|0\right>\!\!\left<0\right| + \left|1\right>\!\!\left<N-1\right|, \nonumber \\ \hat{E}_0\hat{A}^{\dagger} &=& \left|N-1\right>\!\!\left<1\right| + \left|0\right>\!\!\left<0\right|, \label{eq:app95AE095E0Adag95derivation} \end{align}\tag{100}\] which gives \[\begin{align} \hat{E}_1 = 2\left|0\right>\!\!\left<0\right| + \left|1\right>\!\!\left<N-1\right| + \left|N-1\right>\!\!\left<1\right|. \label{eq:app95E195explicit95derivation} \end{align}\tag{101}\] Similarly, \[\begin{align} \hat{E}_0\hat{A} &=& \left|N-1\right>\!\!\left<N-1\right| + \left|0\right>\!\!\left<N-2\right|, \nonumber \\ \hat{A}^{\dagger}\hat{E}_0 &=& \left|N-2\right>\!\!\left<0\right| + \left|N-1\right>\!\!\left<N-1\right|, \label{eq:app95E0A95AdagE095derivation} \end{align}\tag{102}\] and hence, we have \[\begin{align} \hat{E}_2 &=& 2\left|N-1\right>\!\!\left<N-1\right| + \left|0\right>\!\!\left<N-2\right| + \left|N-2\right>\!\!\left<0\right|. \label{eq:app95E295explicit95derivation} \end{align}\tag{103}\] These expressions show explicitly why the DBC correction through \(\hat{P}^{4}\) is boundary-local.

8 Shot-noise and measurement-cost analysis↩︎

This appendix derives the leading finite-shot scaling of the kinetic-energy estimators. All estimates below assume the independent shot sets for different measurement primitives. The correlated or shared-shot estimators can reduce the constant prefactors but do not change the \(M^{-1/2}\) scaling of the root-mean-square error.

8.1 Translation-moment variance↩︎

The Hadamard estimator for \(m_l = {\rm Re}\langle \hat{A}^l \rangle\) produces a binary random variable \(Z_l \in \{-1,+1\}\) with \[\begin{align} \mathbb{E}[Z_l] = m_l, \quad {\rm Var}(Z_l) = 1 - m_l^2. \label{eq:app95binary95translation95variance} \end{align}\tag{104}\] With \(M_l\) shots, \[\begin{align} \widetilde{m}_l = \frac{1}{M_l}\sum_{r=1}^{M_l}Z_{l,r}, \nonumber \\ {\rm Var}(\widetilde{m}_l) = \frac{1 - m_l^2}{M_l} \leq \frac{1}{M_l}. \label{eq:app95m95l95variance} \end{align}\tag{105}\]

For PBC, define \[\begin{align} \alpha_{\rm P} = \frac{\hbar^2}{m\Delta_{\rm P}^{2}}, \quad \beta_{\rm P} = \frac{\hbar^4}{m^3 c^2\Delta_{\rm P}^{4}}. \label{eq:app95aP95bP} \end{align}\tag{106}\] Then, Eq. (36 ) can be written as \[\begin{align} T_{\rm P} = \alpha_{\rm P}(1-m_1) - \frac{ \beta_{\rm P}}{4} (m_2 - 4m_1 + 3). \label{eq:app95TP95linear95moments} \end{align}\tag{107}\] The linear error propagation formula gives \[\begin{align} {\rm Var}\bigl(\widetilde{T}_{\rm P} \bigr) &=& \left(-\alpha_{\rm P}+\beta_{\rm P}\right)^2{\rm Var}(\widetilde{m}_1) + \frac{1}{16} \beta_{\rm P}^2{\rm Var}(\widetilde{m}_2), \label{eq:app95var95TP95exact95linear} \end{align}\tag{108}\] assuming independent estimates of \(m_1\) and \(m_2\). Hence, \[\begin{align} {\rm Var}\bigl(\widetilde{T}_{\rm P} \bigr) \leq \frac{\left(-\alpha_{\rm P}+\beta_{\rm P}\right)^2}{M_1} + \frac{\beta_{\rm P}^2}{16 M_2}. \label{eq:app95var95TP95bound} \end{align}\tag{109}\] If \(M_1\) and \(M_2\) are both proportional to a common shot budget \(M\), then \({\rm RMSE}(\widetilde{T}_{\rm P})=O(M^{-1/2})\).

8.2 Boundary-overlap variance↩︎

A probability \(P\) estimated from \(M\) independent projection measurements has Bernoulli variance \[\begin{align} {\rm Var}(\widetilde{P}) = \frac{P(1-P)}{M} \leq \frac{1}{4M}. \label{eq:app95bernoulli95probability95variance} \end{align}\tag{110}\] For \(B_{fg} = 2P_{fg}^{+}-P_f-P_g\), using the independent estimates of \(P_{fg}^{+}\), \(P_f\), and \(P_g\) with \(M\) shots each gives \[\begin{align} {\rm Var}(\widetilde{B}_{fg}) &=& 4{\rm Var}(\widetilde{P}_{fg}^{+}) + {\rm Var}(\widetilde{P}_{f}) + {\rm Var}(\widetilde{P}_{g}) \leq \frac{3}{2M}. \label{eq:app95Bfg95variance95bound} \end{align}\tag{111}\] The bound is conservative because the endpoint probabilities can often be obtained from the same computational-basis samples used for the potential estimator.

For DBC, define \[\begin{align} a_{\rm D} = \frac{\hbar^2}{m\Delta_{\rm D}^{2}}, \quad d_{\rm D} = \frac{\hbar^2}{2m\Delta_{\rm D}^{2}}, \quad e_{\rm D} = \frac{\hbar^4}{8 m^3 c^2\Delta_{\rm D}^{4}}. \label{eq:app95aD95dD95bD} \end{align}\tag{112}\] Combining Eqs. (44 ) and (45 ), the DBC kinetic estimator can be written as \[\begin{align}T_{\rm D} = a_{\rm D}(1-m_1) + d_{\rm D}b_0 - e_{\rm D}(2m_2 - 8m_1 + 6 + 4b_0 - b_1 - b_2 + b_{00}) \nonumber\\= a_{\rm D}-6 e_{\rm D} + (-a_{\rm D} + 8 e_{\rm D})m_1 - 2 e_{\rm D}m_2 + (d_{\rm D} - 4 e_{\rm D})b_0 + e_{\rm D}b_1 + e_{\rm D} b_2 - e_{\rm D}b_{00}. \label{eq:app95TD95linear95form} \end{align}\tag{113}\] The first term in the second line is a constant shift with respect to the sampled quantities, and therefore it does not contribute to the variance. Assuming independent estimates of the non-constant terms, \[\begin{align}{\rm Var}\left(\widetilde{T}_{\rm D}^{(4)}\right) &=& (-a_{\rm D} + 8e_{\rm D})^2{\rm Var}(\widetilde{m}_1) + (2e_{\rm D})^2{\rm Var}(\widetilde{m}_2) \nonumber\\&&+ (d_{\rm D}-4e_{\rm D})^2{\rm Var}(\widetilde{b}_0) + e_{\rm D}^2{\rm Var}(\widetilde{b}_1) + e_{\rm D}^2{\rm Var}(\widetilde{b}_2) + e_{\rm D}^2{\rm Var}(\widetilde{b}_{00}). \label{eq:app95var95TD95linear} \end{align}\tag{114}\] Each variance on the right-hand side is \(O(M^{-1})\) when \(M\) shots are allocated to each primitive, so \[\begin{align} {\rm RMSE}(\widetilde{T}_{\rm D}^{(4)}) = O(M^{-1/2}). \label{eq:app95TD95RMSE95scaling} \end{align}\tag{115}\] The finite-shot simulations in Sec. 4.4 confirm this scaling.

8.3 Potential-energy sampling variance↩︎

To keep the notation consistent with Proposition 3, we reserve \(\hat{V}_{\tau}\) for the potential operator and use \(\widetilde{V}_{\tau}\) for the Monte Carlo estimator. Let \(J_1,\ldots,J_M\) be independent position-basis samples with \({\rm Pr}(J_r=j)=\left|c_j\right|^2\). Define the sample-level random variables and their average by \[\begin{align} X_{r,\tau} := V\bigl(x_{J_r}^{(\tau)}\bigr), \quad \widetilde{V}_{\tau} := \frac{1}{M}\sum_{r=1}^{M}X_{r,\tau}, \label{eq:app95potential95random95variable} \end{align}\tag{116}\] where \(r=1,\ldots,M\). Then, \[\begin{align} \mathbb{E}\bigl[\widetilde{V}_{\tau}\bigr] &=& \frac{1}{M}\sum_{r=1}^{M}\mathbb{E}[X_{r,\tau}] = \sum_{j=0}^{N-1}V\bigl(x_j^{(\tau)}\bigr)\left|c_j\right|^2 = \langle\hat{V}_{\tau}\rangle, \nonumber \\ {\rm Var}(\widetilde{V}_\tau) &=& \frac{1}{M^2}\sum_{r=1}^{M}{\rm Var}(X_{r,\tau}) = \frac{{\rm Var}(X_{\tau})}{M}. \label{eq:app95potential95variance95exact} \end{align}\tag{117}\] Here, the last equality uses the iid assumption, and \(X_{\tau}\) denotes a generic copy of the single-shot variable \(X_{r,\tau}\). If \(V(x_j^{(\tau)}) \in [V_{\min}, V_{\max}]\), then \[\begin{align} {\rm Var}(\widetilde{V}_\tau) \leq \frac{(V_{\max} - V_{\min})^2}{4M}. \label{eq:app95potential95variance95bound95again} \end{align}\tag{118}\] The total-energy variance is obtained by adding the kinetic and potential variance contributions when independent shot sets are used.

8.4 Measurement settings and circuit scaling↩︎

Through the leading relativistic correction, the PBC kinetic estimator requires two translation-moment settings: \(l=1\) and \(l=2\). The DBC kinetic estimator uses the same two translation settings and adds the boundary-overlap settings for the three pairs \[\begin{align} (0,N-1), \quad (1,N-1), \quad (0,N-2), \label{eq:app95boundary95pairs} \end{align}\tag{119}\] plus the endpoint probabilities, which can be obtained from computational-basis samples. Thus, up to \(\hat{P}^{4}\), the number of boundary-specific DBC settings is independent of \(N\).

The controlled implementation of \(\hat{A}^l\) is a controlled modular increment by \(l\) on an \(L\)-qubit register [26]. For fixed \(l\), a ripple-carry construction has gate count scaling polynomially and typically linearly in \(L\) up to the architecture-dependent constants and ancilla choices. The QFT-based adders give an alternative implementation with a different constant and connectivity profile. The estimator framework does not depend on a particular adder design. It only requires that controlled cyclic translations can be implemented or otherwise estimated as unitary register operations.

9 Boundary-overlap measurement identities↩︎

This appendix gives the probability identities underlying the boundary-overlap estimators. Let \[\begin{align} \left|\psi\right> = \sum_{j=0}^{N-1}c_j\left|j\right>, \quad c_j=\left<{j}|{\psi}\right>. \label{eq:app95overlap95state} \end{align}\tag{120}\] For two distinct computational basis states \(\left|f\right>\) and \(\left|g\right>\), define \[\begin{align} P_f &=& \left|c_f\right|^2, \nonumber \\ P_g &=& \left|c_g\right|^2, \nonumber \\ P_{fg}^{+} &=& \left|\langle{s_{fg}^{+}}|{\psi}\rangle\right|^2, \label{eq:app95Pf95Pg95Pplus} \end{align}\tag{121}\] where \(|s_{fg}^{+}\rangle = ({\left|f\right>+\left|g\right>})/{\sqrt{2}}\). Then, we have \[\begin{align} P_{fg}^{+} &=& \frac{1}{2}\left|c_f + c_g\right|^2 = \frac{1}{2} \left( \left|c_f\right|^2 + \left|c_g\right|^2 + c_f^{\ast} c_g + c_g^{\ast} c_f \right). \label{eq:app95Pplus95expanded} \end{align}\tag{122}\] Consequently, \[\begin{align} 2P_{fg}^{+}-P_f-P_g &=& c_f^{\ast} c_g + c_g^{\ast} c_f = \left<\psi\right| \bigl( \left|f\right>\!\!\left<g\right| + \left|g\right>\!\!\left<f\right| \bigr) \left|\psi\right>. \label{eq:app95real95overlap95identity} \end{align}\tag{123}\] The factor of two multiplying \(P_{fg}^{+}\) is a direct consequence of the normalisation of \(|{s_{fg}^{+}}\rangle\).

The same construction can access imaginary coherences when needed. Define the phase-dependent superposition \[\begin{align} P_{fg}^{\phi} = \left|\langle{s_{fg}^{\phi}}|{\psi}\rangle\right|^2, \label{eq:app95phase95superposition} \end{align}\tag{124}\] where \(|{s_{fg}^{\phi}}\rangle = {(\left|f\right>+e^{i\phi}\left|g\right>)}/{\sqrt{2}}\). Then, we have \[\begin{align} && 2P_{fg}^{\phi} - P_f - P_g = e^{-i\phi}c_f^{\ast}c_g + e^{i\phi} c_g^{\ast} c_f = 2{\rm Re}\bigl(e^{-i\phi} c_f^{\ast} c_g\bigr). \label{eq:app95phase95overlap95identity} \end{align}\tag{125}\] The choice \(\phi=0\) gives the real coherence in Eq. (123 ), while \(\phi=\pi/2\) gives \[\begin{align} && 2P_{fg}^{\pi/2} - P_f - P_g = 2{\rm Im}(c_f^{\ast}c_g) = \left<\psi\right|\bigl( -i\left|f\right>\!\!\left<g\right| + i\left|g\right>\!\!\left<f\right| \bigr)\left|\psi\right>. \label{eq:app95imag95overlap95identity} \end{align}\tag{126}\] The DBC kinetic estimators in the main text require only the real coherence because the boundary operators \(\hat{E}_0\), \(\hat{E}_1\), and \(\hat{E}_2\) are Hermitian combinations of off-diagonal projectors.

Using the shorthand \[\begin{align} B_{f, g} = 2P_{fg}^{+} - P_f - P_g, \label{eq:app95Bfg95again} \end{align}\tag{127}\] the DBC boundary terms through \(\hat{P}^{4}\) are \[\begin{align} \langle\hat{E}_0\rangle &=& B_{0,N-1}, \nonumber \\ \langle\hat{E}_0^2\rangle &=& P_0 + P_{N-1}, \nonumber \\ \langle\hat{E}_1\rangle &=& 2P_0 + B_{1,N-1}, \nonumber \\ \langle\hat{E}_2\rangle &=& 2P_{N-1} + B_{0,N-2}. \label{eq:app95Es95Bfg} \end{align}\tag{128}\] Thus, the boundary correction can be obtained from endpoint probabilities and three real two-state coherences.

10 Numerical methods↩︎

This appendix describes the numerical procedures used for the validation figures in Sec. 4. All numerical simulations use dense linear algebra because the purpose is operator validation rather than large-scale classical performance.

10.1 Matrix construction↩︎

For a given \(L\), \(N=2^L\). The cyclic translation matrix is \[\begin{align} (\hat{A})_{j,k} = \delta_{j,k+1 \;{\rm mod} \;N}, \label{eq:app95A95matrix95elements} \end{align}\tag{129}\] where row and column indices run from \(0\) to \(N-1\). The PBC momentum matrices are \[\begin{align} \hat{P}_{\rm P}^{2} = -\frac{\hbar^2}{\Delta_{\rm P}^{2}} (\hat{A}+\hat{A}^{\dagger}-2\hat{\mathbb{1}}), \quad \hat{P}_{\rm P}^{4} = (\hat{P}_{\rm P}^{2})^2. \label{eq:app95PBC95matrix95construction} \end{align}\tag{130}\] For DBC, the open-chain shift is \[\begin{align} (\hat{B})_{j,k} = \delta_{j,k+1} \quad (0 \leq k \leq N-2), \label{eq:app95B95matrix95elements} \end{align}\tag{131}\] with no wrap-around element. The DBC matrices are \[\begin{align} \hat{P}_{\rm D}^{2} = -\frac{\hbar^2}{\Delta_{\rm D}^{2}}(\hat{B}+\hat{B}^{\dagger}-2\hat{\mathbb{1}}), \quad \hat{P}_{\rm D}^{4} = (\hat{P}_{\rm D}^{2})^2. \label{eq:app95DBC95matrix95construction} \end{align}\tag{132}\] Equivalently, the DBC matrices can be constructed from the cyclic shift using \(\hat{E}_0\) as in Eq. (93 ). Both constructions are numerically identical up to floating-point precision.

The lattice potential is the diagonal matrix \[\begin{align} \hat{V}_{\tau} = \sum_{j=0}^{N-1}V\bigl(x_j^{(\tau)}\bigr)\left|j\right>\!\!\left<j\right|, \quad \tau\in\{{\rm P},{\rm D}\}. \label{eq:app95potential95matrix95construction} \end{align}\tag{133}\] The non-relativistic Hamiltonian used to generate the smooth-potential ground states is \[\begin{align} \hat{H}^{tot}_{\tau,{\rm nr}} = \frac{\hat{P}_{\tau}^{2}}{2m}+\hat{V}_{\tau}. \label{eq:app95Hnr95matrix} \end{align}\tag{134}\] The ground state is obtained by exact diagonalisation of the Hermitian matrix \(\hat{H}_{\tau,{\rm nr}}\).

10.2 Estimator reconstruction from a state vector↩︎

Given a normalised state vector \(\vec{c}=(c_0,\ldots,c_{N-1})^{T}\), the exact translation moments are computed as \(m_l={\rm Re}(\vec{c}^{\,T}\hat{A}^l \vec{c}\,)\). To avoid ambiguity in implementation, this can be written component-wise as \[\begin{align} m_l = {\rm Re}\left( \sum_{j=0}^{N-1} c_j^{\ast}c_{j-l \;{\rm mod} \;N} \right) \quad (l=1,2), \label{eq:app95numerical95ml} \end{align}\tag{135}\] for the convention \(\hat{A}\left|j\right>=\left|j+1 \;{\rm mod} \;N\right>\). The PBC estimator reconstruction is then \[\begin{align} \langle\hat{P}_{\rm P}^{2}\rangle_{\rm est} &=& \frac{2\hbar^2}{\Delta_{\rm P}^{2}}(1-m_1), \nonumber \\ \langle\hat{P}_{\rm P}^{4}\rangle_{\rm est} &=& \frac{2\hbar^4}{\Delta_{\rm P}^{4}}(m_2-4m_1+3). \label{eq:app95PBC95estimator95reconstruction} \end{align}\tag{136}\] For DBC, the boundary quantities are reconstructed as \[\begin{align} b_0 &=& 2{\rm Re}(c_0^{\ast}c_{N-1}), \nonumber \\ b_{00} &=& \left|c_0\right|^2 + \left|c_{N-1}\right|^2, \nonumber \\ b_1 &=& 2\left|c_0\right|^2 + 2{\rm Re}(c_1^{\ast}c_{N-1}), \nonumber \\ b_2 &=& 2\left|c_{N-1}\right|^2 + 2{\rm Re}(c_0^{\ast} c_{N-2}). \label{eq:app95bs95numerical} \end{align}\tag{137}\] The DBC estimator reconstruction is \[\begin{align} \langle\hat{P}_{\rm D}^{2}\rangle_{\rm est} &=& \frac{2\hbar^2}{\Delta_{\rm D}^{2}}(1-m_1) + \frac{\hbar^2}{\Delta_{\rm D}^{2}}b_0, \nonumber \\ \langle\hat{P}_{\rm D}^{4}\rangle_{\rm est} &=& \frac{\hbar^4}{\Delta_{\rm D}^{4}} \left(2m_2-8m_1+6+4b_0-b_1-b_2+b_{00}\right). \label{eq:app95P24D95estimator95reconstruction} \end{align}\tag{138}\] These expressions are compared with direct matrix expectations \[\begin{align} \langle\hat{P}_{\tau}^{2}\rangle_{\rm mat} = \vec{c}^{\,T} \hat{P}_{\tau}^{2} \vec{c}, \quad \langle\hat{P}_{\tau}^{4}\rangle_{\rm mat} = \vec{c}^{\,T} \hat{P}_{\tau}^{4} \vec{c}. \label{eq:app95matrix95expectations} \end{align}\tag{139}\]

10.3 Benchmark parameters↩︎

The PBC free-particle benchmark uses the Fourier mode \(n=1\) with \(R=10\), \(\hbar=m=1\), \(c=2\), and \(L=3,4,\ldots,10\). For each \(L\), the lattice momentum eigenvalue is \[\begin{align} p_{\rm lat}^{2}(n) = \frac{4\hbar^2}{\Delta_{\rm P}^{2}} \sin^2\left(\frac{\pi n}{N}\right), \label{eq:app95pbc95free95plat} \end{align}\tag{140}\] and the continuum value is \[\begin{align} p_{\rm cont}^{2}(n) = \left(\frac{2\pi n\hbar}{R}\right)^2. \label{eq:app95pbc95free95pcont} \end{align}\tag{141}\] The energies plotted in Fig. 4 are \(T_{\rm cont}\), \(T_{\rm lat}\), and \(T_{\rm pert}\) as defined in Eq. (54 ).

The DBC square-well benchmark uses \(N=64\), \(R=10\), \(\hbar=m=1\), and compares modes \(s=1,\ldots,20\). The analytic sine eigenvectors are \[\begin{align} \psi_s(j) = \sqrt{\frac{2}{N+1}} \sin\left(\frac{\pi s(j+1)}{N+1}\right), \label{eq:app95dbc95sine95numerical} \end{align}\tag{142}\] with eigenvalues \[\begin{align} p_{{\rm D},{\rm lat}}^{2}(s) = \frac{4\hbar^2}{\Delta_{\rm D}^{2}} \sin^2\left(\frac{\pi s}{2(N+1)}\right). \label{eq:app95dbc95square95eigenvalue} \end{align}\tag{143}\]

The smooth-potential benchmark uses \(R=10\), \(\hbar=m=1\), \(c=2\), and \(L=4,5,\ldots,10\). For PBC, the potential is \[\begin{align} V_{\rm P}(x) = V_0\left(1-\cos\left(\frac{2\pi x}{R}\right)\right), \label{eq:app95pbc95smooth95potential} \end{align}\tag{144}\] for \(V_0=0.5\). For DBC, the potential is \[\begin{align} V_{\rm D}(x) = \frac{1}{2} m \omega^2\left(x-\frac{R}{2}\right)^2, \label{eq:app95dbc95smooth95potential} \end{align}\tag{145}\] with \(\omega=0.4\). For each grid, we diagonalise \(\hat{H}^{tot}_{\tau,{\rm nr}}\) and evaluate \(E_{{\rm nr},\tau}\), \(\Delta E_{{\rm rel},\tau}\), and \(E_{{\rm rel},\tau}^{(4)}\) on the non-relativistic ground state, as in Eq. (60 ).

The finite-shot benchmark uses the smooth-potential ground states at \(L=5\) with the same physical parameters, allocates \(M\) shots to each measurement primitive, and estimates the RMSE over independent Monte Carlo repetitions. The discretized ground states \(\left|\psi_{0,\tau}\right>\) are represented by trigonometric and Gaussian wavefunctions for PBC and DBC, respectively.

10.4 Error metrics↩︎

For a momentum moment, we use the relative discrepancy \[\begin{align} \epsilon_{P^{2r}} = \frac{ \left|\langle\hat{P}_{\tau}^{2r}\rangle_{\rm est} - \langle\hat{P}_{\tau}^{2r}\rangle_{\rm mat}\right|}{\max\left\{ \left|\langle\hat{P}_{\tau}^{2r}\rangle_{\rm mat}\right|, \; \epsilon_{\rm floor}\right\}}, \quad r=1,2, \label{eq:app95relative95moment95error} \end{align}\tag{146}\] where \(\epsilon_{\rm floor}\) is a small numerical floor used only to avoid a zero denominator. For total energies, we use the absolute reconstruction error \[\begin{align} \epsilon_E = \left|E_{\rm est} - E_{\rm mat}\right|. \label{eq:app95absolute95energy95error} \end{align}\tag{147}\] For the weakly relativistic diagnostic, we monitor \[\begin{align} \eta_{\tau} = \frac{\langle\hat{P}_{\tau}^{2}\rangle}{m^2c^2}. \label{eq:app95eta95tau95diagnostic} \end{align}\tag{148}\] This state-dependent quantity does not bound all higher moments, but it is a useful indicator that the low-momentum expansion is being applied in the intended regime. A more conservative diagnostic is the spectral-support condition \[\begin{align} \eta_{\tau, {\rm max}} = \max_{\lambda \in {\rm supp}(\psi)}\frac{\lambda}{m^2 c^2}, \label{eq:app95spectral95eta} \end{align}\tag{149}\] where \(\lambda\) runs over eigenvalues of \(\hat{P}_{\tau}^{2}\) that have non-negligible weight in the state \(\left|\psi\right>\).

11 Higher-order extensions and additional benchmarks↩︎

This appendix records two natural extensions of the framework: higher-order relativistic moments and higher-order finite-difference stencils. These extensions are not required for the main validation, but they clarify how the present construction generalises.

11.1 Including the \(\hat{P}^{6}\) relativistic correction↩︎

The next term in the positive-energy expansion is \[\begin{align} \hat{T}^{(6)}_\tau = \frac{\hat{P}_\tau^2}{2m} - \frac{\hat{P}_\tau^4}{8 m^3 c^2} + \frac{\hat{P}_\tau^6}{16 m^5 c^4}, \label{eq:app95T695definition} \end{align}\tag{150}\] where \(\hat{P}_\tau^6=(\hat{P}_\tau^2)^3\). For PBC, using \(\hat{K}=\hat{A}+\hat{A}^{\dagger}-2\hat{\mathbb{1}}\) and \(\hat{P}_{\rm P}^{2}=-(\hbar^2/\Delta_{\rm P}^{2})\hat{K}\), one finds \[\begin{align} \hat{P}_{\rm P}^{6} = \frac{\hbar^6}{\Delta_{\rm P}^{6}} \left[ 20\hat{\mathbb{1}} - 15\left(\hat{A}+\hat{A}^{\dagger}\right) + 6\left(\hat{A}^{2}+\hat{A}^{\dagger 2}\right) - \left( \hat{A}^{3}+\hat{A}^{\dagger 3} \right) \right]. \label{eq:app95PBC95P695operator} \end{align}\tag{151}\] Hence, \[\begin{align} \langle\hat{P}_{\rm P}^{6}\rangle = \frac{\hbar^6}{\Delta_{\rm P}^{6}} \left( 20 - 30m_1 + 12m_2 - 2m_3 \right), \label{eq:app95PBC95P695moments} \end{align}\tag{152}\] where \(m_3={\rm Re}\langle\hat{A}^{3}\rangle\). Thus, the \(\hat{P}^{6}\) term requires one additional translation moment under PBC.

For DBC, \[\begin{align} \hat{P}_{\rm D}^{6} = -\frac{\hbar^6}{\Delta_{\rm D}^{6}}\left(\hat{K}-\hat{E}_0\right)^3. \label{eq:app95DBC95P695formal} \end{align}\tag{153}\] Because \(\hat{K}\) and \(\hat{E}_0\) do not commute, the expansion contains all ordered products of \(\hat{K}\) and \(\hat{E}_0\). The resulting correction terms remain local near the endpoints, but they involve a wider boundary stencil than the \(\hat{P}^{4}\) correction. In particular, terms can connect endpoint basis states to sites up to three lattice steps away. This is the natural higher-order analogue of the \(\hat{E}_0\), \(\hat{E}_1\), and \(\hat{E}_2\) terms appearing at fourth order.

11.2 Higher-order finite-difference stencils↩︎

The main text uses the standard second-order finite-difference approximation to \(-\partial_x^2\). A fourth-order accurate PBC stencil for the momentum squared is \[\begin{align} \hat{P}_{\rm P,4th}^{2} = \frac{\hbar^2}{\Delta_{\rm P}^{2}} \left[ \frac{5}{2}\hat{\mathbb{1}} - \frac{4}{3}\left(\hat{A}+\hat{A}^{\dagger}\right) + \frac{1}{12}\left(\hat{A}^{2}+\hat{A}^{\dagger 2}\right) \right]. \label{eq:app95fourth95order95stencil} \end{align}\tag{154}\] This reduces the finite-grid discretization error for smooth states but requires additional translation moments even for the non-relativistic kinetic term. Squaring this operator to obtain \(\hat{P}^{4}\) introduces moments up to \(\hat{A}^{4}\). Under DBC, the same stencil requires the boundary corrections that remove wrap-around couplings associated with both \(\hat{A}\) and \(\hat{A}^{2}\). Thus, higher-order stencils trade discretization accuracy for a larger but still structured observable set.

References↩︎

References↩︎