Quantum computing is an emerging technology that harnesses the laws of quantum mechanics to deliver unprecedented computational power [Hid19, NC02, Pre18, RP11]. Quantum algorithms operate on an $n$-qubit Hilbert space with dimension $2^{n}$, offering vast multidimensional spaces for computational models. Hence, it has a unique capability to handle large-scale scientific computing problems. A natural application is the PDEs from time-dependent Schrödinger equations (TDSE), which follow unitary evolutions and hence the wave functions can be coherently represented on quantum computers. Known as Hamiltonian simulations, a variety of efficient algorithms have been developed toward this end [BCC⁺15, BCK15, BCC⁺17, LC17, LC19, CGJ19, KSB19, AFL21, JLY22b, AFL22, FLT23].

While the TDSE is formulated in the entire physical space, in practice, the computation has to be done in a bounded domain, typically in locations where the electron density is concentrated. For situations when the electrons are being emitted outside the computational domain, such as the photoionization process, or when they are being drawn away, such as the ionization process, this approach, however, can result in an extraordinarily large computational domain. An artificial boundary condition (ABC) that may absorb outgoing wave packets and, more importantly, keep the size of the computational domain to a minimum, is frequently used to solve such problems. A correctly set up ABC will yield the same result as if the computation is done in the infinite domain, e.g., see the survey paper [AAB⁺08].

Due to the introduction of the ABCs, the computer simulation is now following a dynamics that is no longer unitary. As a result, existing Hamiltonian simulation techniques can not be directly applied. The goal of this paper is to bridge this emerging gap, by mapping the non-Hamiltonian dynamics back to a Schrödinger equation or Hamiltonian system, to enable the immediate Hamiltonian simulation capability. This is accomplished by using the recently introduced Schrödingerization technique [JLY22a, JLY22b]. It is a generic framework to convert any linear partial (and ordinary) differential equation or dynamical system into a system of Schrödinger equations or a Hamiltonian system by going to the space which is one dimensional higher, via the warped phase transformation. A key advantage of the Schrödingerization approach is its simplicity and generality, where the general Hamiltonian the system evolves under has a simple form. When the boundary value problem is considered, resulting in an inhomogeneous system that is not Hamiltonian, one just needs to introduce one auxiliary variable to make the system homogeneous, and then the Schrödingerization process still applies [JLY22a, JLY22b].

In this paper, we explore three popular types of ABCs, including the complex absorbing potential method, perfect-matched layers method, and Dirichlet-to-Neumann approach, to captivate the quantum implementation utilizing the Schrödingerization approach.

Another potential alternative to implementing the ABCs is by using quantum linear solvers (QLS) [HHL09, CKS17]. This can be done by first performing spatial discretization that reduces the models to a system of linear ODEs, followed by time-marching schemes. Runge-Kutta methods, spectral (collocation) methods [CL19], multi-step methods [BCOW17], and Dyson series [Kro22], are some of the important existing methods. In some earlier works [CL19], the analysis for this type of algorithm assumes that the matrix $A$ in the linear ODE system is diagonalizable, and the final complexity involves a condition number of the eigenvector matrix. This is later improved by Krovi [Kro22], where the condition number is replaced by the following bound,

In the context of ABCs, due to the dissipative nature, the matrix $A$ is typically stable, in that the real parts of the eigenvalues are non-positive. Indeed, such stability property is a key focus in the development and analysis of ABCs [EM79]. The stability implies that $C(A)=1$, thus completely removing the dependence of the complexity on the condition number.

However, QLS-based methods still require more involved steps like a truncation of Taylor series expansions or the construction of a unitary operation associated with the inverse of a matrix. The latter is often realised through the block-encoding formalism [GSLW19]. While these results might give similar query and gate complexity to the Schrödingerization approach, it is not straightforward to actually show explicitly how the inverse is obtained in the block-encoding formalism, which makes it difficult to perform in practice. This is in contrast to the Schrödingerization method, which is a much simpler scheme where the Hamiltonian to be simulated is explicitly given, which places the problem directly in the realm of Hamiltonian simulation. The cost in the Schrödingerization method is also independent of the condition number.

We note that there are also very recent alternatives to QLS for simulating amplitude-encoded solutions of differential equations, for instance schemes based on block-encoding [ALWZ22] or linear combination of unitaries (LCU) [ALL23]. The application of these methods to ABCs could provide a complementary viewpoint to our methods here.

We first review several types of artificial boundary conditions that are commonly used in practice. The goal of these numerical techniques is to be able to simulate the dynamics in the entire space $\mathbb{R}^{d}$, i.e.,

where $\hat{H}$ is a Hamiltonian operator expressed in reduced units,

We assume that the initial wave function $\psi(\bm{x},0)$ is compactly supported in a subdomain, here denoted by $\Omega\subset\mathbb{R}^{d}$. On the other hand, due to the initial momentum, or the influence of the external potential, the wave function can spread well beyond $\Omega$, as shown in Figure 1. To avoid the need for continuously expanding the computational domain, it is necessary to develop suitable boundary conditions that can propagate wave packets beyond the domain’s boundaries. These conditions, known as artificial boundary conditions (ABCs), are intended to prevent boundary reflections, and hence the interference with the dynamics inside the domain.

From the proceeding section, we observe that the TDSE under the three types of ABCs can be expressed in the following general form,

Here $\ket{\phi(t)}$ is a quantum state that encodes the wave function, possibly the auxiliary functions, in the subdomain. $H_{0}$ and $H_{1}$ are both Hermitian, and $H_{0}$ is a slight modification of the original Hamiltonian.

The Schrödingerization approach was first proposed in [JLY22b]. The key step is the warped phase transformation $\bm{v}(t,p)=e^{-p}\bm{\phi}(t)$, defined for $p>0$ and symmetrically extends the initial data to $p<0$. This transformation converts (16) to a system of linear convection equations:

On the other hand, the solution $\bm{\phi}(t)$ can be restored by a simple integration,

A more intuitive view is by discretizing the $p$ dimension and concatenating the corresponding function for each $p$. Toward this end, we observe that $\bm{v}(t,p)$ has an exponential decay in the $p$ space. Hence we may choose a sufficiently large interval, such that,

Next, we choose sub-intervals with size $\Delta p=(b-a)/M$ with $M=2N$ being an even number and the grid points denoted by $a=p_{0}<p_{1}<\cdots<p_{M}=b$. Then the mapping from $\bm{v}$ to $\psi$ can be implemented using a quadrature formula,

To compute $\bm{v}(t,p),$ let the vector $\bm{w}$ be the collection of the function $\bm{v}$ at these grid points, defined more precisely as follows,

with “;” indicating the straightening of $\{\bm{w}_{i}\}_{i\geq 1}$ into a column vector. This can also be expressed as a superposition state using $\ket{k}$ as a new basis,

By applying the discrete Fourier transformation in the $p$ direction, one arrives at,

Here, $P_{\mu}$ is the matrix expression of the momentum operator $-i\partial_{p}$, given by

where $\mu_{l}=2\pi l/(b-a)$ are the Fourier modes and

is the matrix representation of the discrete Fourier transform.

At this point, we have successfully mapped the dynamics back to a Hamiltonian system. By a change of variables $\tilde{\bm{w}}=(I_{u}\otimes\Phi^{-1})\bm{w}$, one has

This is more amenable to an approximation by a quantum algorithm. In particular, if $H_{0}$ and $H_{1}$ are sparse, then (20) is a Schrödinger equation with the Hamiltonian $H=H_{0}\otimes I-H_{1}\otimes D_{\mu}$ that inherits the sparsity.

We now present results from some numerical experiments. As a concrete example, we consider the TDSE in two dimensions, and choose a square domain $\Omega:=[-3,3]\times[-3,3].$ Following the example in [Zhe07], we choose the initial condition as follows,

Here $r=\sqrt{x^{2}+y^{2}}.$ For the potential, we choose $V(x,y)=\sin(2\pi r)$ in the unit disc and extend it to zero outside.

In the numerical experiments, we solve (20) using the matrix exponential:

This is implemented on a classical computer using MATLAB’s built-in function expm.m. In the implementation, we also take $N_{x}=N_{p}=64$ and $p\in[L,R]=[-5,5]$. It is important, however, to note that the purpose of these experiments is to demonstrate that the equation (20) in the Schrödingerization captures the dynamics under the ABCs. A quantum implementation of (20) is much more feasible than classical computers, due to the less dependence on the dimension.

We first examine the CAP method. In these experiments, the imaginary potential $W(x,y)$ is chosen as the sum of $w_{1}(x,y)$ and $w_{2}(x,y)$, with

as shown in Figure 2. The results from quantum simulations of (3) and (20) at $t=0.3,0.6$ and $t=0.9$ are displayed in Figure 3(d-f), from which we observe very similar results to the ones given in Figure 3(a-c) for the original form (3). This verifies the correctness of the protocol of the Schrödingerization approach.

Now we turn to the PML method. We implemented (7) on the same test problem, and the results are depicted in Figure 4. Here the function $\sigma$ is chosen as

To implement the DtN approach, one can solve a discrete boundary element equation [WL20] to determine the matrix $R$ in (12), which then determines the self-energy in (27). Since our purpose is to test the Schrödingerization approach, we compute $R$ from (12) by choosing a relatively large domain $D\supset\Omega$ and simply set $\Omega_{\text{I\!I}}=D\backslash\Omega_{\text{I}}$. For the test example, it suffices to choose $D=[-6,6]^{2}$.

Since all the matrices ($H_{\text{I},\text{I}}$, $H_{\text{I},\text{I\!I}}$, and $H_{\text{I\!I},\text{I\!I}}$) can be directly extracted from the Hamiltonian matrix on $D$, we can compute the matrices in (27) by direct matrix inversion and multiplications. In (15), we take $s_{0}=1$. The direct simulation with the DtN ABC can be done by solving the ODEs (26). The results are displayed in the top row in Figure 5. Meanwhile, a quantum algorithm would solve (20) with the two matrices $H_{0}$ and $H_{1}$ from (29). We observe that the classical algorithm is able to maintain the propagating pattern of the wavefront without boundary reflections. The solution based on the Schrödingerization has a similar performance. The slight difference can be attributed by the discretization in the $p$ space.

Once the quantum dynamics with ABCs is turned into a Hamiltonian system (20), one can apply a Hamiltonian simulation algorithm to produce the wave function $\ket{\psi(T)}.$ Due to the possible time-dependence potential $V(\bm{x},t)$, a Hamiltonian simulation algorithm for time-dependent Hamiltonian $H(t)$ should be considered. To assess the algorithm complexity associated with the implementation of the ABCs, we use the recent results by Berry et al. [BCS⁺20, Theorem 10], although the other algorithms can also be used for the assessment. Here we simply highlight the query complexity,

Here the norm is defined as,

We now discuss the complexity with implementing the time-dependent Schrödinger equation (20) obtained from the Schrödingerization of the three types of ABCs. Since the discretization error for equation (17) is $\mathcal{O}(\Delta p)$ (since $e^{-|p|}$ is only continuous but not continuously differentiable in $p$) and the discretization in real space has error $\mathcal{O}(\Delta x^{2})$, we set

to make the discretization error comparable.

Thus we can conclude that

where we have used $\|D_{\mu}\|=\mathcal{O}(1/\Delta p)$ together with (32). Meanwhile, the sparsity of $H$ is $s(H)=\mathcal{O}(1)$ (For $d$ dimensions, it is $\mathcal{O}(d)$). Combining these estimates gives the total complexity, as shown in the first row of Table 1.

Inserting these estimates into (29), we have

The other factor that contributes to the complexity is the sparsity of the matrices in (27). Numerical observations of these matrices suggest that they are close to identity matrices. For example, for the matrix $\Sigma_{\Gamma,\Gamma}$ in (27), we plotted in Figure 6 its entries by removing entries that are much smaller than its maximum norm. We let $s_{\Sigma}$ be such an effective sparsity and the overall complexity is linear in $s_{\Sigma}$.

In this article, we use the Schrödingerization approch recently introduced in [JLY22a, JLY22b] for quantum dynamics with articifial boundary conditions (ABCs). While a quantum dynamics with ABC is no longer a Hermitian Hamiltonian system, which is most natural for quantum simulation, the Schrödingerization approach makes it so in a simple fashion. We give the implementation details and estimate the computational complexities for three representative ABCs, including the complex absorbing potential method, perfectly matched layer methods, and Dirichlet-to-Neumann approach. Our numerical experiments validate this apporach, demonstrating that the Schrödingerized systems yield the same results as the original dynamics with ABCs.

Our approach provides a simple quantum simulation algorithms for quantum dynamics in unbounded domains. As pointed out in [JLY22a, JLY22b], the Schrödingerization apporach is not only applicable to quantum dynamics with ABCs, it also applies to general linear partial differential equations with physical boundary conditions. It can also be applied to such equations with interface conditions. These will be the subject of our future research.

SJ was partially supported by the NSFC grant No. 12031013, the Shanghai Municipal Science and Technology Major Project (2021SHZDZX0102), and the Innovation Program of Shanghai Municipal Education Commission (No. 2021-01-07-00-02-E00087). NL acknowledges funding from the Science and Technology Program of Shanghai, China (21JC1402900). YY was partially supported by China Postdoctoral Science Foundation (no. 2022M712080). XL has supported by a Seed Grant from the Institute of Computational and Data Science (ICDS) at Penn State.