If we consider only the diagonal terms, where $i = j$ in Eq. (<ref>), we need to perform a linear combination of $O(2d)$ terms, for the each pair of particles in the summand of Eq. (<ref>). Therefore, the total gate cost for a diagonal convective operator is $O(d^2p^2 N\log(N)^2)$ Toffoli or simpler gates. The circuit cost to implement the diagonal terms of the convective operator for a scalar pair-wise potential is an additional $O(\eta^2)$ controlled applications of the block encoding of the convective operator, which results in the final gate cost that is

\begin{equation}
O(\eta^2 d^2p^2 N\log(N)^2)
\end{equation}

Toffoli or simpler gates to block encode the diagonal terms of the convective operator. Now, if we express $N$ in terms of $\epsilon$ as $N \sim \epsilon^{-1/p}$, we obtain

\begin{equation}
O(\eta^2 d^2 \epsilon^{-1/p}\log(1/\epsilon)^2)
\end{equation}

Toffoli or simpler gates needed to implement the desired block encoding.

§  DISCUSSION AND OUTLOOK

The simulation of physical systems is

one of the most promising use cases of future quantum computers. In this work, we have addressed the block encoding of a class of elliptic operators that arise in such simulations, when discretized on a Cartesian grid of $N$ grid points per dimension in $d$ spatial dimensions and satisfying. We provide explicit block encodings for the standard Laplacian operator in $d$ dimensions with Dirichlet, Neumann, and Robin boundary conditions that require only $O(d\log(N))$ Toffoli gates.
However, the error convergence is restricted to low order, i.e., $\epsilon \sim N^{-2}$(Dirichlet/Periodic) or $N^{-1}$(Neumann/Robin), using standard approaches. We provide these block encodings as linear combinations of unitaries, where the unitaries are simply the modular inc(dec)rement operators and reflections and their controlled versions, and can be implemented with $O(\log(N))$ Toffolis. For the Dirichlet and Neumann boundary value problems, we show how to obtain high-order convergence rates through the use of the periodic extension, so that for any $p>3$ we can bound the global truncation error by $\epsilon \sim O(N^{-p+1})$ for these boundary value problems. We also address the block encoding of operators with boundary conditions on an irregular geometry (e.g., not given as the simple Cartesian product of closed intervals), assuming that we can efficiently compute a membership oracle marking to which part of the computational domain a specific grid point belongs providing the same overall complexity for generalized rectangular domains in $\mathbb{R}^d$.

We then show how to block encode a many-body differential operator by first block encoding a multi-particle potential $V$ given as the sum of pair-wise interactions. We focus on the evaluation of the Lennard-Jones potential, which involves calculating terms of the form $V(r) = \frac{1}{r^{12}} - \frac{1}{r^6}$, where $r$ is the interparticle distance. We provide efficient quantum circuits to block encode a radially symmetric “discrete position operator" in <ref>, which computes the squared Euclidean distance for every tuple of grid points in $d$-dimensions. We address the singular behavior of the potential near $r=0$ through the use of a cutoff parameter $\delta$ and polynomial approximation over the whole domain. We show an alternative technique using a piecewise polynomial approximation in conjunction with an inequality test, we can find faster polynomial approximations to the function than with uniform polynomial approximations. This technique improves the dependence on $\delta$ in the polynomial approximation and thereby reduces the total number of queries to the block encoding. Finally, we show how to block encode the convective operator given as $\nabla V \cdot \nabla$, which describes particles interacting with a force field $\nabla V$ where $\nabla$ is the $\eta d$-dimensional gradient operator.

Although the block encodings we construct in this work are provided as linear combinations of unitaries (LCUs), additional work is still needed to convert these LCUs into 1- and 2-qubit gates. In particular, we do not specify the construction of circuits that prepare the LCU coefficients, nor the construction of the controlled applications of the unitaries. More work is needed to address the constant factors in the complexity resulting from the controlled applications of the 1 sparse unitaries we use in the LCU.
Preparation of the LCU coefficients are not a bottleneck for our block encoding, and concrete circuit implementations will be addressed in future work.
We note that generation of software to automatically compile and optimize the quantum circuit implementations of these block encodings would

be a useful development.
More generally, it may also be worthwhile to explore Cartesian embedding methods <cit.> in connection with the methods presented in Sec. <ref>, to construct efficient circuits for representing differential operators on curved domains. However, we point out that if one wishes to obtain generic high-order accuracy on irregular geometries, the $hp$-FEM approach is a much more robust framework than embedding methods.
Block encoding the FE matrix is much more involved, and is the subject of future work.

This work is one step toward the goal of developing efficient quantum subroutines for realistic simulation tasks involving key differential operators on future quantum computers. The resulting circuits are conceptually simple and rely on the use of well-known primitives such as shift and reflection gates to block encode the Laplacian matrix, as well as on more modern methods based on inequality testing and quantum signal processing to efficiently block encode a more general elliptic operator such as $\nabla V(\mathbf{x})\cdot \nabla$, with $V$ a scalar multi-particle potential. This work thus provides a blueprint for the block encoding of many-body differential equations.

As such, we expect these constructions will be useful for analyzing the quantum resources to simulate high-dimensional many-body systems of interest in the physical and applied sciences.

T.D.K., A.M.A., K.M., and K.B.W. were supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research under Award Number DE-SC0023273.
T.D.K was also supported by a Siemens FutureMakers Fellowship.
J.-P.L. was supported by the National Science Foundation (PHY-1818914, CCF-1729369), the NSF Quantum Leap Challenge Institute (QLCI) program (OMA-2016245, OMA-2120757), and a Simons Foundation award (No. 825053). T.D.K would also like to acknowledge Ian Convy, Di Fang, Cory Hargus, and Torin Stetina for many insightful discussions. We also thank Philipp Schleich and Torin Stetina for helpful comments on an early version of the manuscript.

§  CIRCUIT IMPLEMENTATIONS

§.§ Modular increment and decrement

The modular left and right shift operators implement the following transformation on $N = 2^n$ states

\begin{equation}
\begin{split}
S^{-1}\ket{l} &= \ket{(l-1)\mod N}\\
S^1\ket{l} &= \ket{(l+1)\mod N}.
\end{split}
\end{equation}

In this section, we will construct the quantum circuit for the incrementer operation to provide some understanding of the underlying circuit structure.

Let's first consider the simplest example, which is just a single qubit. In this case the shift-by-one operator is just the bitflip operator, or the Pauli-$X$ gate

\begin{equation}
\ket{b_0} \mapsto \ket{b_0\oplus1}.
\end{equation}

which is realized as the circuit

$\ket{b_0}$     X

The two qubit case $\ket{b_1,b_0}\mapsto \ket{b_1\oplus b_0, b_0\oplus 1}$, which flips the state of qubit $b_1$ if the state of $b_0$ is 1 and flips $b_0$. Therefore $\ket{00} \mapsto \ket{01} \mapsto \ket{10} \mapsto \ket{11} \mapsto \ket{00}$. This circuit is realized via a controlled-X gate and an X gate:

$\ket{b_1}$     01

$\ket{b_0}$     0     X

The three qubit case is similar. This is given by the mapping

\begin{equation*}
\ket{b_2,b_1,b_0}\mapsto \ket{b_2\oplus b_1 b_0 , b_1\oplus b_0, b_0\oplus 1}
\end{equation*}

, which flips the 2nd qubit if the $b_0 = b_1 = 1$, flips the first qubit if $b_0 = 1$, and flips $b_0$. We can see that
$\ket{000}\mapsto \ket{001}\mapsto \ket{010} \mapsto \ket{011}\mapsto \ket{100} \mapsto \ket{101} \mapsto \ket{110} \mapsto\ket{111}$
which is realized via the circuit:

$\ket{b_2}$     1 1

$\ket{b_1}$     1             01

$\ket{b_0}$     0             0           X

The four qubit case is similar, the algebraic representation is

\begin{equation*}
\ket{b_3,b_2,b_1,b_0}\mapsto\ket{b_3\oplus b_2b_1b_0, b_2\oplus b_1b_0, b_1\oplus b_0, b_0\oplus 1}
\end{equation*}

and the circuit is

$\ket{b_3}$    11

$\ket{b_2}$     2    1 1

$\ket{b_1}$     1    1             01

$\ket{b_0}$     0    0             0           X

Carrying on in a similar fashion, for larger numbers of qubits will prepare the $R$ operator for $n-$qubits. The implementation for the left-shift-by-one (i.e. decrement) operator is similar, and can be obtained by replacing the control on 1 with a control on 0. The total circuit depth is $\mathcal{O}(n)$. However, the multi-control Toffoli gates used in this scheme need to be compiled down to Toffoli gates, which can then be further compiled to one and two qubit gates. In <cit.>, it is shown that the direct compilation of this circuit will require $O(n^2)$ Toffoli or simpler gates. However, we can use the scheme provided in <cit.> to implement this quantum circuit using $O(n)$ Toffoli gates, at the cost of an ancilla qubit.

§.§ Reflection Operators

Let $A \subset [N]$ mark states in the standard computational basis. Let $R_A$ be reflection about the states in $A$, that is

\begin{equation}
R_A = -2\sum_{j\in A}\ket{j}\bra{j} + I
\end{equation}

which is just

\begin{equation}
R_A\ket{i} = \begin{cases}
1 & i \notin A\\
-1 & i \in A
\end{cases}
\end{equation}

Then the number of elementary gates to implement $R_A$ is $O(|A| \mathcal{C}(C^n(Z)))$ where $\mathcal{C}(C^n(Z))$ is the cost to implement the $n$-qubit controlled $z$-gate.

Let $a \in A$ and $a\in\{0,1\}^n$ is an $n-$bit string so that we may write $a = \sum_{j=0}^{n-1} a_j 2^j \simeq \{a_0, a_1, \ldots, a_{n-1}\} \equiv b$. Let $S[a] = \{j; a_j = 0 \forall a_j \in a\}$ be the set of all unset bits of $a$. Then we let $X_{S[a]} = \bigotimes_{i\in[n]}X^{\neg a_i}$ be the operator that performs an $X-$gate on every qubit that corresponds to an unset bit in the bitstring $a$. Additionally, let $C^{n}(Z)$ be the n-fold controlled $Z$ operation,

\begin{equation}
C^{n}(Z) = \sum_{i \in [N-1]}\ket{i}\bra{i} - \ket{N-1}\bra{N-1}.
\end{equation}

Then $\bra{d_0}\bra{d_1}\cdots\bra{d_{n-1}}X_{S[a]}C^{n}(Z)X_{S[a]}\ket{d_0}\ket{d_1}\cdots\ket{d_{n-1}}$ satisfies

\begin{align*}
&=\bra{d_0\oplus \neg a_0}\bra{d_1\oplus \neg a_1}\cdots\bra{d_{n-1}\oplus \neg a_{n-1}}C^{n}(Z)\ket{d_0\oplus \neg a_0} \ket{d_1\oplus \neg a_1}\cdots \ket{d_{n-1}\oplus \neg a_{n-1}}\\
\end{align*}

But then, since $\bra{i}C^n(z)\ket{i} = -1 \iff i = N-1$, implies that in order for the input state $\ket{d} = \ket{d_0}\cdots\ket{d_{n-1}}$ to be in the marked subspace, it must satisfy:

\begin{equation}
\begin{aligned}
d_0\oplus \neg a_0 &= 1\\
d_1 \oplus \neg a_1 &= 1\\
d_{n-1}\oplus \neg a_{n-1} &= 1
\end{aligned}
\end{equation}

which implies

\begin{equation}
\begin{aligned}
d_0 &= a_0\\
d_1 &= a_1\\
d_{n-1} &= a_{n-1}
\end{aligned}
\end{equation}

therefore, the operations above prepared the reflection about the marked state $a$.

We can form a product of these reflections to get a reflection over all the elements in a subspace marked by elements of a set. Since the number of elements in the set is $|A|$, we expect this process needs to be performed for $O(|A|)$ times. Resultantly, the cost to implement this operation is $O(|A| \mathcal{C}(C^n(z)))$ where $\mathcal{C}(C^n(z)) = O(n)$ is the cost to implement the controlled-Z on n-qubits in terms of Toffoli and T gates.

§  BLOCK ENCODING COSTS

§.§ Verifying LCUs for differential operators

Let $L$ be the $d-$dimensional finite difference Laplacian with a $p = 2a+1$ point scheme on $N$ grid points. And let $A = A_0 \cup A_1 \cdots \cup \cdots A_{d-1}$ be the set of all boundary points over the $d$-dimensional domain with each $A_i \subset [N]$. We additionally assume that each $A_i$ is a connected set in the sense that $\forall x,y \in A_i$, $\exists i \leq |A_i|$ such that $x\pm i \equiv y \mod N $ and we also assume that each $A_i$ can be specified without dependence on the coordinates in $A_j$ for every $j\neq i$.

We define the shifted set

\begin{equation}
A_{i+c} = \{j + c \mod N: j\in A_i\}.
\end{equation}

Then the $1d$ Laplacian $L_d$ with boundary nodes $A_d$ can be specified via the LCU

\begin{equation}
\frac{1}{2}\sum_{j=-a}^a r_j\left(S^j + S^jR_{A_d-j}\right) + \frac{1}{2}(I + R_{A_d}).
\label{eq:Laplacian LCU appdx}
\end{equation}

And the $d$-dimensional Laplacian can be encoded using the direct sum

\begin{equation}
L = \bigoplus_{i=0}^{d-1} L_i,
\end{equation}

where $L_i$ is the Laplacian formed by separable connected Dirichlet boundary conditions in each direction $i$.

Our goal is to show that Eq. (<ref>) results in a matrix acting as the identity on the boundary points and the Laplacian on the interior points. Then, we enforce the assumed condition that the boundary nodes can be specified in each dimensions independently to combine these one-dimensional operators into the appropriate $d-$dimensional operator.

The finite-difference, periodic $L$ operator on $N$ nodes with is expressible as the sum

\begin{equation}
L = \sum_{i=-a}^{a} r_i S^i
\end{equation}

where the coefficients $r_i$ are given by Eq. (<ref>) and $S^i$ are the shift operators. Now, consider $A_d$ a set of boundary points on the $d$th spatial coordinate. Then,

\begin{align*}
&\sum_{j=-a}^a(r_jS^j + r_jS^jR_{A_d-j})\\
&=\sum_{j=-a}^a(r_jS^j + r_jS^j(-2\sum_{k\in A_d} \ket{k-j}\bra{k-j} -I))\\
&=\sum_{j=-a}^ar_j \sum_{i\in[N]}(\ket{i+j}\bra{i} + \ket{i+j}\bra{i}(-2\sum_{k\in A_d} \ket{k-j}\bra{k-j} + I))\\
\end{align*}

\begin{align*}
&\sum_{i\in[N]}(-2\sum_{k\in A_d} \ket{i+j}\bra{i}\ket{k-j}\bra{k-j} + \ket{i+j}\bra{i}) \\
&=-2\sum_{i\in[N]}\sum_{k\in A_d} \ket{i+j}\delta_{i,k-j}\bra{k-j} + \sum_{i\in[N]}\ket{i+j}\bra{i} \\
&=-2\sum_{k\in A_d} \ket{k}\bra{k-j} + \sum_{i\in[N]}\ket{i+j}\bra{i}\\
\end{align*}

\begin{align*}
&\sum_{i\in[N]}\ket{i+j}\bra{i} +  \ket{i+j}\bra{i}(-2\sum_{k\in A_d} \ket{k}\bra{k-j} + I))\\
&= \sum_{i\in[N]}\ket{i+j}\bra{i} - 2\sum_{k\in A_d} \ket{k}\bra{k-j} + \sum_{i\in[N]}\ket{i+j}\bra{i}\\
&= 2\sum_{i\in[N]}\ket{i+j}\bra{i} - 2\sum_{k\in A_d} \ket{k}\bra{k-j}\\
&\propto \sum_{i\in[N]}\ket{i+j}\bra{i} - \sum_{k\in A_d} \ket{k}\bra{k-j}\\
\begin{cases}
0 & \forall \text{ } k \in A_d :\hspace{.2cm} \ket{i+j}\bra{i} = \ket{k}\bra{k-j}\\
1 & \text{else}
\end{cases}
\end{align*}

\begin{align*}
&\bra{m}\sum_{i\in[N]}\ket{i+j}\bra{i} - \sum_{k\in A_d} \ket{k}\bra{k-j}\ket{n}\\
&=\sum_{i\in[N]}\delta_{i+j,m}\delta_{i,n} - \sum_{k\in A_d} \delta_{m,k}\delta_{n,k-j}\\
&=\delta_{n+j,m} - \delta_{m,A_d} \delta_{n,m-j}\\
&=\delta_{n+j,m} - \delta_{m,A_d} \delta_{n+j,m}\\
0 & m \in A_d \text{ and } n+j = m\\
1 & m \notin A_d \text{ and } n+j=m
\end{cases}
\end{align*}

which is the desired behavior.

So, we may now write:

\begin{align*}
&\sum_{j=-a}^ar_j\left(\sum_{n,m\in[N]}\delta_{n+j,m} - \sum_{n\in[N]}\sum_{m\in A_d}\ \delta_{n+j,m}\right)\\
&=\sum_{j=-a}^ar_j\left(S^j - \sum_{n\in[N]}\sum_{m\in A_d}\ \delta_{n+j,m}\right)\\
&=\sum_{j=-a}^ar_j\left(S^j - \sum_{n\in[N]}\sum_{m\in A_d} \delta_{n+j,m}\right)
\end{align*}

which more plainly now, we can see that this serves to zero out the entries on the rows supported in $A_d$. Now we turn to the particular case of the diagonal:

\begin{align*}
r_0[I - \sum_{m\in A_d}\delta_{m,m}]_{jj} = \begin{cases}
r_0 & m \notin A_d\\
0 & m \in A_d.
\end{cases}
\end{align*}

Now, if we add another projector onto $A_d$ we find,

\begin{equation}
[r_0(I - \sum_{m\in A_d}\delta_{m,m}) + \Pi_{A_d}]_{ii} = \begin{cases}
r_0 & i \notin A_d\\
1 & i \in A_d
\end{cases}.
\end{equation}

Then, we can decompose the projector $\Pi_{A_d} = I - \Pi_{A_d^c}$ as

\begin{equation}
\Pi_{A_d} = \frac{1}{2}\left(I - R_{A_d^c}\right) = \frac{1}{2}\left(I + R_{A_d}\right)
\end{equation}

to obtain the final LCU

\begin{equation}
L_d =\frac{1}{2}\sum_{j=-a}^a r_j\left(S^j + S^jR_{A_d-j}\right) + \frac{1}{2}(I + R_{A_d}).
\label{eq:LCU for BCs}
\end{equation}

Resultantly, the final matrix has the form

\begin{equation}
L_d =
\begin{pmatrix}
{L}^{(l)}&  & \\
& I_{|A_d|} & \\
& & {L}^{(r)}
\end{pmatrix}
\end{equation}

as desired.

The number of elementary gates needed to implement the LCU given by Eq. (<ref>) is $O(p|A|\mathcal{C}(C^n(Z)))$ Toffoli. Furthermore, the above implements an $(O(1), \log(2p),0)$ block encoding of the desired matrix.

The proof of this theorem can be found in appendix <ref> <ref>

The linear combination provided via Eq. (<ref>) requires $2p$ shift matrices and $p$ reflection matrices. Using the construction from <cit.> gates, the shift operators can be encoded using $O(n)$ Toffoli-or-smaller gates with an additional ancilla. We assume these to be the bottleneck for implementing these operations, so we will count the complexity in terms of the total number of Toffoli gates needed. Therefore, the cost to implement these shift gates is $O(np)$ Toffoli or simpler gates.

The reflection matrices are diagonal matrices with $\pm 1$ on the diagonal, with $-1$'s on rows who's indices correspond to states of the subspace being reflected over,

\begin{equation}
R_{A}\ket{j} =\begin{cases}
\ket{j} & j \notin A\\
-\ket{j} & j \in A
\end{cases}
\end{equation}

as given in Lemma <ref>, the cost to implement the quantum circuit for the reflection is $O(|A| \mathcal{C}(C^n(Z))p)$. So that the total number of Toffoli operations needed is $O(p|A|\mathcal{C}(C^n(Z))) = O(np)$.

The subnormalization factor is the $1-$norm of the coefficients in the LCU.

\begin{equation}
\lambda = \sum_{j=-a}^a |r_j| + 1 \leq \frac{2\pi^2}{3} + 1 = O(1)
\end{equation}

where the first inequality can be found in Lemma 6 of Ref <cit.>. Then, in $d$-dimensions, using $O(dp + \log(d))$ ancilla qubits to form the linear combination

\begin{equation}
\frac{1}{d}\sum_{i\in[d]}U_{L_j}
\end{equation}

where $U_{L_j}$ is the block encoding of the direction $j$ Laplacian.

So in total, this requires $O(dp+\log(d))$ ancilla, $\alpha = O(\frac{1}{d})$, a number of Toffoli or smaller gates that scales as $O(dnp\mathcal{C}(C^n(Z)) |A|)$ and in time $O(np\mathcal{C}(C^n(Z)))$. And similarly for $\eta$ particles another multiplicative factor of $\eta$; $\alpha = O(\frac{1}{\eta d}) < \frac{\pi^2}{\eta d}$, $T = O(\eta dnp\mathcal{C}(C^n(Z)) |A|) $, but add no cost to overall simulation time as they can be performed in parallel.

§  SUCCESS PROBABILITY FOR DIRICHLET BLOCK ENCODING

In constructing the system of equations resulting from discretization of the affiliated differential operators, one must encode in the right hand side $\mathbf{b}$ the initial data from the boundary conditions, and additionally place a $1$ on the diagonal of the matrix corresponding to the node where the boundary data is known.

On the quantum computer, the workflow begins by preparing a quantum state $\ket{b}$ on the system register that corresponds to the right hand side of the differential equaiton. For the Dirichlet problem, one might have

\begin{equation}
\begin{aligned}
L[u] &= 0\\
u(0) &= 0\\
u(1) &= 1,
\end{aligned}
\end{equation}

where $\ket{b}$ is simply the basis vector $\ket{N-1}$ after appropriate normalization. Then, the next step is to prepare the block-encoding and apply it to $\ket{b}$. Now we seek to obtain a bound on the success probability in terms of $N$. The first step is the application of $A$, the matrix representation of the problem, onto $\ket{b}$. The success probability of this first step is crucial to the efficiency of the algorithm.

We seek to obtain a lower bound to

\begin{equation}
||\mathcal{L}\ket{b}||^2 = |\bra{b}\mathcal{L}^\dagger \mathcal{L}\ket{b}|,
\label{eq: success prob eqn}
\end{equation}

which we then divide by $d$ times the one dimensional subnormalization factor

\begin{equation}
\alpha_d = d\sum_{j=-a}^a|r_j|
\end{equation}

We will have $d^2$ terms in (<ref>). $d$ of these terms are given by

\begin{equation}
L_0^\dagger L_0 \otimes I^{d-1}, I\otimes L_1^\dagger L_1 \otimes I^{d-2}, \ldots, I^{\otimes d-1}\otimes L_d^\dagger L_d.
\end{equation}

Since each of these operators are $1$-local it suffices to consider only the the action of $L_j^\dagger L_j$ on the subspace corresponding to the $d$th dimension, as the other terms in the tensor product contribute a multiplicative factor of $1$.

We have for the $j$th term of this form

\begin{align*}
\frac{1}{N_j}\sum_{k,l \in B_j}\bra{k}L_j^\dagger L_j\ket{l}\ket{l}\\
\end{align*}

but since each $k \in A_j$ is by definition either on the boundary or interior of the boundary region, $L$ is diagonal in that region so that $L\ket{k} = -r_0\ket{k}$, so this allows us to obtain a bound

\begin{align*}
\vert{\frac{1}{N_j}\sum_{k,l\in A_j} \bra{k}L^\dagger L \ket{l}\vert} \geq \vert\sum_{i = 1}^{N_j}\frac{r_0^2}{N_j}\vert =r_0^2
\end{align*}

which is a lower bound that is independent of $N_j$. Doing this for each of these $d$-terms, we find a lower bound for the norm of these one-local terms as

\begin{equation}
\end{equation}

The rest of the terms are 2-local with $L_j$ acting on one site $L_k^\dagger$ acting on another site. If both $L_k$ and $L_j^\dagger$ are acting on each of their respective boundary registers, of which there are also $d$ of these terms $\mathcal{D}_k\ket{0}$ and $\mathcal{D}_j\ket{0}$ then the contribution will be of the same order given above, i.e.

\begin{equation}
\end{equation}

If only one of $L_j$ or $L_k^\dagger$ is acting on it's respective boundary ket, and the other is operating on another boundary register $\mathcal{D}_i\ket{0}$ then we expect to see a contribution of the order $\frac{1}{N_i}$, of which there are $d$ of. This gives an additional contribution that is

\begin{equation}
\end{equation}

Finally, if neither $L_j$ nor $L_k^\dagger$ is acting on their respective registers which there are $2(d^2-d)$ terms like this, we can expect an additional contribution that is $O(1/(N_l N_i))$, which gives a contribution that is

\begin{equation}
O\left(2d(d-1)\sum_{i,j=0}^{d-1}\frac{1}{N_i N_j}\right).
\end{equation}

Therefore the total success probability will be the sum of each of these terms

\begin{equation}
P_{\text{succ}} = \frac{1}{d^2\alpha^2}\left(2dr_0^2 + d\sum_{j=0}^{d-1}\frac{1}{N_j} + 2d(d-1)\sum_{i,j=0}^{d-1}\frac{1}{N_i N_j}\right).
\end{equation}

We can therefore obtain a lower bound to the success probability as

\begin{equation}
P_{\text{succ}} > \frac{2dr_0^2}{d^2 \alpha} > \frac{6r_0^2}{2d\pi^2} > \frac{6 \left(\frac{2\pi ^2}{3}\right)^2}{2d\pi^2} = \frac{4 \pi ^2}{3 d}
\end{equation}

where the second to last and last inequalities are direct applications of the bound given in (<ref>). Which is constant in $N$ and inverse linear in $d$. Knowing the particular boundary geometry allows one to obtain tighter bounds.

d = 1
The success probability can be estimated by computing $||A\ket{b}||$. Let's consider an example in one spatial dimension for one particle, where the boundary data specifies that $u(0) = 0$ and $u(1) = 1$. Here, the system of equations we form has the structure

\begin{equation}
\begin{bmatrix}
-2 & 0 & 0&0& \cdots & 0\\
1 & -2 & 1 &0&\cdots & 0\\
0 & 1 & -2 &1& \cdots & 0\\
\vdots &\vdots &\ddots & \ddots &\ddots & \vdots\\
0 & 0 & \ldots& 1 & -2 & 1\\
0 & 0 & \ldots& 0 & 0 & -2
\end{bmatrix}
\begin{bmatrix}
\vdots \\
\end{bmatrix}
\begin{bmatrix}
\vdots \\
\end{bmatrix}.
\end{equation}

Presented using the shift operators, this equation becomes

\begin{equation}
\left[-2I + \frac{1}{2}\left(S^1 + S^1_{N-2,N-1}\right) + \frac{1}{2}\left(S^{-1} + S^{-1}_{0,1}\right)\right]\ket{u(x)} = -2 \ket{N-1}
\label{eq:1d alg equation}
\end{equation}

We divide the matrix on the left hand side by the LCU subnormalization factor

\begin{equation}
\alpha = \sum_{i=1}^5 |c_i| = 4.
\end{equation}

Furthermore, we require the right hand side also be normalized, so we replace $-2\ket{N-1}$ with $\ket{N-1}$.

For the 1D case, we compute

\begin{equation}
||L\ket{N-1}||^2 = |\bra{N-1}L^\dagger L \ket{N-1}| = 5
\end{equation}

so that the success probability

\begin{equation}
P_{\text{succ}} = \frac{1}{\alpha^2}||L\ket{N-1}||^2
\end{equation}

\begin{equation}
\frac{5}{16} \sim .3125
\label{eq:1d success prob}
\end{equation}

which is a constant.

In $2$ dimensions, the boundary condition becomes

\begin{equation}
u(x=L,y) = 1
\end{equation}

which, upon discretization corresponds to the right hand side vector becoming

\begin{equation}
\ket{b} = \ket{N-1}\frac{1}{\sqrt{N}}\sum_{i=0}^{N-1}\ket{i}
\end{equation}

and the matrix representation of the operator also becomes

\begin{equation}
\frac{1}{2\alpha}\left(L\otimes I + I \otimes L\right).
\end{equation}

As above, we wish to bound

\begin{equation}
||\left(L\otimes I + I \otimes L\right)\ket{b}||^2 = |\bra{b}\left(L\otimes I + I \otimes L\right)^2\ket{b}|,
\end{equation}

after which we divide by $10$ in order to obtain the success probability with an LCU block-encoding. An equivalent expression for the norm is given by

\begin{equation}
|\bra{b}\left(L^\dagger L \otimes I + I \otimes L^\dagger L + L^\dagger \otimes L + L\otimes L^\dagger\right )\ket{b}|
\end{equation}

Taking this term by term, the first term

\begin{align*}
&\bra{b}L^\dagger L \otimes I\ket{b} = \bra{N-1}\frac{1}{N}\sum_{i=0}^{N-1}\bra{i}(L^\dagger L\otimes I)\ket{N-1}\sum_{j=0}^{N-1}\ket{j}\\
&=\frac{1}{N}\left(\bra{N-1}L^\dagger L\ket{N-1}\times \sum_{i,j}\bra{i}\ket{j}\right)\\
&=\frac{1}{N}\left(5\times \sum_{i,j}\bra{i}\ket{j}\right)\\
&=\frac{5}{N}\left(\times N\right)\\
\end{align*}

the second term

\begin{align*}
&\bra{b}I \otimes L^\dagger L \ket{b} = \bra{N-1}\frac{1}{N}\sum_{i=0}^{N-1}\bra{i}(I\otimes L^2)\ket{N-1}\sum_{j=0}^{N-1}\ket{j}\\
&=\frac{1}{N}\sum_{i,j=0}^{N-1}(\bra{N-1}\ket{N-1}\bra{i}L^\dagger L\ket{j})\\
&=\frac{1}{N}\sum_{i,j=0}^{N-1}(\bra{i}L^\dagger L \ket{j})\\
&= \frac{1}{N}\times8\\
\end{align*}

and the third term

\begin{align*}
&\bra{b}L^\dagger\otimes L\ket{b} = \frac{1}{N}\sum_{i,j=0}^{N-1}(\bra{N-1}L^\dagger \ket{N-1} \bra{i}L\ket{j})\\
\end{align*}

Here we used the facts that $\bra{N-1}L^\dagger \ket{N-1} = -2$ and $\sum_{ij}L_{ij} = -4$.

The fourth term is similar,

\begin{align*}
&\bra{b}L\otimes L^\dagger\ket{b} = \frac{1}{N}\sum_{i,j=0}^{N-1}(\bra{N-1}L \ket{N-1} \bra{i}L^\dagger\ket{j})\\
\end{align*}

Adding up each of the individual contributions, we obtain our estimate to the success probability as

\begin{equation}
5 + \frac{16}{N} + \frac{8}{N} = 5+\frac{24}{N}
\end{equation}

which, is again a constant. The LCU subnormalization  contributes a factor of $\frac{1}{(2\alpha)^2}$, to the success probability giving

\begin{equation}
P_{\text{succ}} = \frac{5 + \frac{33}{2N}}{64} \sim .078 \text{ as } N\rightarrow \infty
\end{equation}

We learned that the overlaps for each of these terms respectively satisfy

\begin{align*}
L^\dagger L \otimes I &\rightarrow 5\\
I\otimes L^\dagger L &\rightarrow \frac{8}{N}\\
L^\dagger\otimes L \text{ or h.c. } &\rightarrow \frac{8}{N}
\end{align*}

We work through the case where $d=3$, this time we have

\begin{equation}
\mathcal{L} = L\otimes I \otimes I + I\otimes L \otimes I + I\otimes I \otimes L
\end{equation}

\begin{equation}
\ket{b} = \ket{N-1}\otimes \frac{1}{N}\sum_{i,j}\ket{i}\otimes\ket{j}.
\end{equation}

We have terms of the form $L^\dagger L\otimes I \otimes I$ and $I\otimes L^\dagger L \otimes I$, in addition, we have terms of the form
$L^\dagger \otimes L \otimes I$ + h.c. and $I\otimes L^\dagger \otimes L$ + h.c.

We can use the results from the previous table, to get

\begin{equation}
\bra{b} L^\dagger L \otimes I \otimes I \ket{b} = 5
\end{equation}

\begin{equation}
\bra{b}(I\otimes L^\dagger L \otimes I + I\otimes I \otimes L^\dagger L)\ket{b} = \frac{16}{N}
\end{equation}

\begin{align}
&\bra{b}(L^\dagger \otimes I \otimes L + L^\dagger\otimes L \otimes I + h.c.)\ket{b} \nonumber\\
&=\frac{8}{N} + \frac{8}{N} +\frac{8}{N} + \frac{8}{N} \nonumber\\
&= \frac{32}{N}.
\end{align}

The final terms are the ones where we have two Laplacian operators acting on the uniform superposition sites, and none acting on the register containing $\ket{N-1}$

\begin{align}
&\bra{b}(I\otimes L \otimes L^\dagger + I\otimes L^\dagger \otimes L)\ket{b}.
\end{align}

It suffices to consider

\begin{align*}
&\frac{1}{N^2}\sum_{ijkl}\bra{i}\bra{j}(L^\dagger \otimes L)\ket{k}\ket{l}\\
&=\frac{1}{N^2}\sum_{ijkl}\bra{i}L^\dagger \ket{k} \bra{j}L\ket{l}\\
&=\frac{1}{N^2}\sum_{ik}\bra{i}L^\dagger \ket{k}\sum_{jl} \bra{j}L\ket{l}\\
\end{align*}

of which there are two of them, so we obtain

\begin{equation}
\bra{b}(I\otimes L \otimes L^\dagger + I\otimes L^\dagger \otimes L)\ket{b} = \frac{32}{N^2}
\end{equation}

Now, we add up all the contributions, and we obtain

\begin{equation}
||\mathcal{L}\ket{b}||^2 =  5 + \frac{48}{N} + \frac{32}{N^2}.
\end{equation}

We then divide through by the subnormalization factor $\frac{1}{(3\alpha)^2} = \frac{1}{144}$ to obtain

\begin{equation}
P_{\text{succ}} = \frac{1}{144}( 5 + \frac{33}{N} + \frac{32}{N^2})
\end{equation}

which in the large $N$ limit gives a success probability that is

\begin{equation}
P_{\text{succ}} = \frac{5}{144} \sim .0347
\end{equation}

§.§ Arbitrary dimension

Here we consider the general case of computing

\begin{equation}
\end{equation}

\begin{equation}
\mathcal{L} = \underset{d}{\underbrace{L\otimes \underset{d-1}{\underbrace{I \otimes I \otimes \cdots \otimes I}} + \cdots +I\otimes I \otimes \cdots \otimes L}}
\end{equation}

\begin{equation}
\ket{b} = \frac{1}{\sqrt{N^{d-1}}} \ket{N-1}\otimes \sum_{i_1,i_2,\ldots,i_{d-1}}\ket{i_1}\ket{i_2}\cdots\ket{i_{d-1}}.
\end{equation}

We will have one term of the form $L^\dagger L$ acting on $\ket{N-1}$, we will have (d-1) terms as $L^\dagger L$ acting on a uniform superposition state, and $2\binom{d}{2}$ terms with action an operator $L$ or $L^\dagger$ acting on two sites. Of the “weight two" terms, we have $2d$ terms where either $L$ or $L^\dagger$ is acting on the $\ket{N-1}$ state with an $L$ or $L^\dagger$ acting anywhere else and we have $2d^2 - 2d$ terms where both sites are acting on the uniform superposition registers.

The first term gives

\begin{equation}
\bra{N-1}L^\dagger L\ket{N-1} = 5
\end{equation}

the second kinds of terms, i.e. the weight 1, no support on 0th-direction register, of which there are $d-1$ of, each contribute $\frac{8}{N}$ in total contribute:

\begin{equation}
\frac{8(d-1)}{N}
\end{equation}

the third kinds of term where there is only $L$ or $L^\dagger$ acting on the 0th direction register and a random $L$ elsewhere, there are $2d$ of those, each contribute $\frac{8}{N}$ in total:

\begin{equation}
\frac{16d}{N}
\end{equation}

the fourth kinds of terms where neither $L$ nor $L^\dagger$ has support on the $\ket{N-1}$ state, each contribute $16/N^2$ of which there are $2(d^2-d)$ of, gives

\begin{equation}
\frac{32(d^2-d)}{N^2}
\end{equation}

Adding these up, we find

\begin{equation}
|\bra{b}\mathcal{L}\ket{b}|^2 = 5 +  \frac{8(d-1) + 16d}{N} + \frac{32(d^2-d)}{N^2},
\end{equation}

so that the success probability, which is multiplied by $\frac{1}{(4d)^2}$ becomes

\begin{equation}
P_{\text{succ}} = |\bra{b}\mathcal{L}\ket{b}|^2 =\frac{5 +  \frac{8(d-1) + 16d}{N} + \frac{32(d^2-d)}{N^2}}{d^2}
\end{equation}

whenever $d>2$.

In the large $d$ limit, we obtain

\begin{equation}
P_{\text{large dim}} = \frac{33}{16N^2}
\end{equation}

and in the large $N$ limit, we obtain

\begin{equation}
P_{\text{large grid}} = \frac{5}{16d^2}
\end{equation}

and in the $d = O(N)$ limit, we obtain

\begin{equation}
p_{d\sim N} = \frac{62N - 41}{16N^3} = O(1/N^2)
\end{equation}

For a large grid where $N \gg d$, we see that the expected scaling for the success probability is $O(1/d^2)$ with the preconstant $5/16 = 5/\alpha^2$.