Efficient solution of the non-unitary time-dependent Schr¨odinger equation on a quantum computer with complex absorbing potential

We explore the possibility of adding complex absorbing potential at the boundaries when solving the one-dimensional real-time Schr¨odinger evolution on a grid using a quantum computer with a fully quantum algorithm described on a n qubit register. Due to the complex potential, the evolution mixes real-and imaginary-time propagation and the wave function can potentially be continuously absorbed during the time propagation. We use the dilation quantum algorithm to treat the imaginary-time evolution in parallel to the real-time propagation. This method has the advantage of using only one reservoir qubit at a time, that is measured with a certain success probability to implement the desired imaginary-time evolution. We propose a specific prescription for the dilation method where the success probability is directly linked to the physical norm of the continuously absorbed state evolving on the mesh. We expect that the proposed prescription will have the advantage of keeping a high probability of success in most physical situations. Applications of the method are made on one-dimensional wave functions evolving on a mesh. Results obtained on a quantum computer identify with those obtained on a classical computer. We finally give a detailed discussion on the complexity of implementing the dilation matrix. Due to the local nature of the potential, for n qubits, the dilation matrix only requires 2 n CNOT and 2 n unitary rotation for each time step, whereas it would require of the order of 4 n +1 C-NOT gates to implement it using the best known al-gorithm for general unitary matrices.


Introduction
The time evolution of quantum systems is a basic ingredient of simulating microscopic physics.It is a common subject of numerous studies, in domains ranging from many-body systems [1,2], electronphonon interactions [3], to quantum field theory [4], or even hydrodynamics [5].Most of these simulations deal with the real-time evolution of unitary systems, which fits naturally into the framework of quantum simulation.Quantum algorithms have also been proposed to perform imaginary-time evolution (ITE) or to mix real-and imaginary-time evolution on quantum computers.Indeed, many quantum systems can be described as or reduced to systems whose norm is not conserved over time (dissipative systems, open systems, tunneling effect, instantons, . . .), and imaginary time propagation methods are widespread in many domains.Extracting the ground state energy of microscopic systems by relying on the exponential decay of excited states is only one of the many applications of these powerful methods, which have been applied with great success not only in particle and nuclear physics but also in condensed matter and quantum chemistry.Developing algorithms on quantum computers able to efficiently encode evolution operators -unitary or not -is thus crucial to simulate a broad range of microscopic systems on quantum devices.Quantum algorithms aiming at including imaginary-time evolution enter either into the class of fully quantum algorithms or into the category of techniques referred to as hybrid quantum-classical methods [6].The Quantum Imaginary-Time Evolution (QITE) method [7][8][9][10][11], which is based on the use of variational principles [12,13] enters into the second class.While expected to be less resilient to noise, fully quantum algorithms have been the subject of extensive efforts in recent years.Among the methods that are now being explored, one can mention the dilation method [14][15][16][17][18][19][20] allowing to treat dissipative systems, the probabilistic ITE [21,22], the "forward and backward" real-time evolution [23] or the recently proposed Singular Value Decomposition (SVD)-based approach of Ref. [24].Other methods, like the "Linear Combination of Unitaries" [25,26] can be used, but they might require a rather large number of ancillary qubits to implement the imaginary potential.
In this work, our primary motivation is to explore the simulation of the Schrödinger equation in position space using quantum computers, and relying completely on a quantum computer algorithm, i.e., avoiding quantum-classical hybrid methods.Although solving Schrödinger-like equation is rather standard on classical computers, especially on low dimension, and despite the fact that the strategy to perform such an evolution on quantum computer is known from decades [27,28], the practical simulation of such problems remains, even today, rather difficult, at least due to the quantum resources required and due to the algorithm complexity (see, for instance, the discussion in [29]).Here we explore the possibility of adding complex absorbing potential (CAP) to a real-time evolution.This method is standard in classical computing to reduce the numerical resources when the evolution is solved on a grid.It consists of adding an imaginary potential that absorbs the wave function escaping from a certain region of interest.The corresponding time-dependent Schrödinger equation then becomes: that corresponds to a mixing of real-and imaginarytime propagation.Such potential reduces the boundary conditions effects and, when the absorption is properly made, allows for a reduction in the number of mesh points needed to simulate the evolution accurately.Specific discussions and technical aspects related to the implementation of CAP on classical computers can be found in for instance in Ref. [30][31][32][33][34][35][36].
A similar advantage can be envisaged when the same equation is solved on a grid on quantum computers, when the grid points are encoded on the qubit's register.Said differently, if we can include absorbing boundary conditions on a quantum computer, we might significantly reduce the number of qubits needed to accurately treat the system.This reduction is quite attractive nowadays because the number of qubits on quantum platforms is rather restricted.In addition, efficient treatment of boundary conditions through absorption opens the perspective to treat certain phenomena like particle decay in real-time or scattering processes.A challenge in adding absorption at the boundaries is that the evolution becomes non-unitary.
We note that the inclusion of CAP in the context of quantum computation has been briefly discussed in [10] using QITE and applied to a simple model case.In [37], although no application with CAP was shown, it was proposed to use the notion of quantum pixel for grid simulation in position.
Here, we explore the possibility of using dilation techniques to simulate real-time evolution given by Eq. ( 1), including CAP in one dimension on a grid.This technique requires only one additional reservoir qubit while using fully quantum algorithms -that is without resorting to any hybrid quantum-classical computation.The application of the imaginary evolution operator is subject to a certain success probability, which is directly related to the "loss of norm" of the system when measuring the ancillary qubit.Special attention is given to improving the quantum algorithm to optimize the success probability.For long time evolution of systems weakly coupled to their environment or short time evolution of systems strongly coupled to their environment, the success probability of our algorithm is expected to remain close to one in many physical situations of interest.
This article is organised as follows.In section 2, we detail our algorithm to simulate a non-unitary evolution given by Eq. ( 1) using the dilation method.Section 3 presents the practical implementation illustrated with results obtained on an IBM quantum simulator.The complexity associated with the implementation of the dilation method is discussed in section 4, followed by concluding remarks in section 5.

Non-unitary propagation algorithm 2.1 Single time-step Trotter-Suzuki decomposition with imaginary time
In order to perform the evolution given by Eq. (1), we first rewrite it as a generalized propagator (we take the convention ℏc = 1): where H is a non-unitary operator, which can be written without loss of generality as H = K +V is the usual Hamiltonian that contains the kinetic (K) and potential (V ) terms, while W is an hermitian operator associated with the CAP.|Ψ(0)⟩ is the initial state, assumed to be known and normalized to 1.To implement the evolution, we consider the first time step ∆t and make use of the first order Trotter-Suzuki approximation [38,39] and write As in usual Trotter propagation, we can estimate the error induced by replacing (2) by (4).For a small time step, we have: where, having in mind solving Eq. ( 1), we implicitly assumed that [W, V ] = 0, which is true if both V and W are diagonal in position space.This will be our case since we will consider V and W local in space.
The last two factors appearing in Eq. ( 4) corresponds to unitary propagators and can usually be treated using standard techniques on digital computers (see illustrations in section 3 and discussion in appendix A).Note that the kinetic-term spectral norm is explicitly given in Eq. ( 23), while ∥W ∥ and ∥V ∥ identify simply with the maximal absolute values these potentials take on the position mesh.
Once the propagator is approximated by Eq. ( 4), a major difficulty in implementing it on a quantum computer stems from the fact that the norm of the wave function immediately becomes lower than one.Indeed, for a single time step, we have where, in the last expression, we used the fact that W is Hermitian.Assuming further that W is positive definite will therefore lead to a decrease of the norm.Such a loss, which could be seen as a loss of information on the system, is impossible if the system is treated on a isolated quantum register together with a perfect quantum computers, where only unitary transformations are possible.A natural solution to this problem is to add one or several qubits that act as a reservoir and can absorb parts of the information on the system by interacting with it.

General description of the method
We briefly recall here the dilation method that is becoming today a common method for implementing non-unitary operators on a system encoded on a qubit register [14][15][16][17][18][19][20].Let us assume that we want to implement |Ψ⟩ → M |Ψ⟩ where M is not unitary, but is assumed to be hermitian and satisfying ∥M ∥ ≤ 1.We assume the system is encoded on n qubits, leading to a computational basis of size 2 n .The dilation method relies on doubling the size of the computational basis by adding a single reservoir qubit and applying the following unitary matrix to the n + 1 qubits: Specifically, we have Circuit for the implementation of the dilation method where the reservoir+system register is prepared in the state |0r⟩ ⊗ |Ψ⟩.After applying the U matrix and measuring 0 in the reservoir qubit, the system state collapses to the state given in Eq. (6).
where {|0 r ⟩, |1 r ⟩} denote the two states associated to the reservoir qubit while |Ψ⟩ is the initial system state.
From the above relation, we see that we can perform the desired operation by preparing the initial state |0 r ⟩ ⊗ |Ψ⟩, applying the U matrix and measuring 0 in the reservoir state.The procedure is schematically represented in the circuit displayed in Fig. 1.More precisely, when measuring the reservoir qubit in state |0⟩, according to the Born measurement's rules, the resulting reservoir+system state is given by: p s is nothing but the probability to measure the reservoir qubit in state |0⟩, called hereafter simply success probability and given by:

Application to imaginary-time propagation
We use the dilation method to apply the CAP, i.e., M ∝ e −W ∆t .A similar problem was addressed in Ref. [20] where the dilation method was implemented to perform imaginary-time propagation using directly the full Hamiltonian, i.e., assuming W = H, to obtain the ground state of the problem.In that case, it was proposed to use the prescription We first implemented this prescription but realized that it has the clear drawback that the success probability exponentially tends to 0 when performing several time steps.This could indeed easily be seen, considering that for arbitrary time-step ∆t, we always have p s ≤ 1/2 provided that M is positive definite.If the dilation process is iterated r times to simulate the evolution up to t = r∆t, the success probability will be lower than 1/2 r .This aspect is critical for practical implementation since, after several steps, most measurements of the ancillary qubit will be rejected, and the method rapidly becomes inefficient.
This unwanted feature can easily be avoided by taking the simpler prescription M = e −W ∆t leading to The unitarity of the matrix U can be easily proven noting that [ √ I − e −2W ∆t , e −W ∆t ] = 0. Applying this prescription to a given initial state |Ψ ini ⟩, we see that the success probability given by Eq. ( 7) becomes: and provided the initial state is normalized to 1 at initial time, the success probability remains close to 1 if ∆t is small, which is much more favorable than the prescription of Ref. [20].
In terms of the specific steps taken for time evolution Eq. ( 1), one can take advantage of the Trotter decomposition(4) as follows: • For the first time step, the propagators e −iV ∆t e −iK∆t are directly implemented as a set of unitary gates on the system register.
• The CAP term is then implemented using the dilation method by measuring the reservoir qubit in state |0⟩.After the measurement, the wave function, denoted hereafter as |Ψ dil (t)⟩, is given by: where U approx is the evolution matrix from Eq.
(5).We denote here by p s (1) (resp.p s (r)) the success probability associated with the first (resp.the r th ) application of the dilation method.
• Iterating the above procedure, after the r th application, and denoting t = r∆t, we deduce that the system wave function is given on the qubit register as: where the global probability of success P s (t = r∆t) relates to r consecutive successful performances of the non-unitary evolution with CAP, i.e. the probability of obtaining r times zero when measuring the reservoir qubit.This probability is given by

Physical estimates of the total success probability evolution
A schematic representation of the iterative procedure is depicted in Fig. 2, where the norm of the state |Ψ dil (t)⟩ is shown as a function of time.Each time a zero is measured in the reservoir qubit, the norm is reinitialized to 1.
Denoting by |Ψ(t)⟩ the wave function obtained with good accuracy by solving numerically the Schrödinger equation ( 1) on a classical computer, according to Eq. ( 5), we see that we have: within the accumulated Trotter errors.Consequently, we also have: The last property makes the dilation method rather attractive for physical systems simulation.Indeed, we have shown that the success probability at a given time t tends, in the limit ∆t → 0, to the norm of the wave function that survives to the absorption when solving the problem on a classical computer.In many physical situations where CAP is useful, we are interested in systems that are rather localized in a certain region of space and where part of the wave function is emitted.This is, for instance, what happens when particles are emitted from a compact localized quantum object.In this case, most particles are emitted over certain time-scale t ≤ τ decay , implying that N (t) decreases and then reaches a stationary asymptotic value.Accordingly, we expect in such situation that the global probability of success P s (t > τ decay ) defined in equation (11) remains constant, which implies that single-step success probability p s (r) ≃ 1 for r ≫ τ decay /∆t.In brief, for such systems where most particles are emitted after a certain transient time, almost all events will be successful after this time.Note that, in section 3, we will actually consider the free wave-packet propagation, the physical situation where all particles escape the grid at infinite time.
This case corresponds to an extreme example of decay phenomena, where the importance of absorption is more pronounced since N (t) tends to 0 at infinite time.It is therefore a perfect example for testing the absorption of particles at the grid boundaries.

Expectation values of non-absorbed observables
In general, we are often interested in computing the expectation values of observables, generically denoted by O taken on a state |Ψ(t)⟩.Provided that we have |Ψ dil (t)⟩ encoded on the system register, we can compute approximately the quantity: by implementing either a Hadamard test [40] and performing a set of measurements, or, if the observable is written as a linear combination of Pauli chains, one can directly measure the system register after a proper change of the measurement basis to estimate each individual Pauli chain (see, for instance, [41], table 3).
In parallel, we can also estimate the total success probability, denoted by P ′ s (t) simply by counting the number of events where only zeros were measured in the reservoir qubit and divide this number by the total number of events.From these estimates, we can compute approximate expectation values on |Ψ(t)⟩ simply as P ′ s (t)⟨O⟩ dil .
3 Illustration: 1D Schrödinger evolution on a grid with complex absorbing potential We give here a proof of principle of applying the strategy discussed above.Specifically, we show examples of 1D Schrödinger equation on a grid with absorbing boundary conditions.We first concentrate on one critical situation where V (x) = 0.In this case, assuming that the particle is located inside the grid and will freely spread during the evolution, it should eventually completely be absorbed by the CAP in the infinite time limit.We then apply our method to the case of a wavepacket in a gaussian potential, where part of the wavefunction remains trapped in the potential over a long period of time.All calculations have been made below using the Qiskit software that emulates a perfect digital computer [42].

Free wave propagation with CAP
We consider a problem described on a restricted area of space x ∈ [x min , x max ].The wave function is assumed to be initially a Gaussian wave function located in the middle of the region with: , (16) with x 0 = (x max + x min )/2.Eventually, we will also consider the possibility that the wave function has an initial boost proportional to the parameter v.In the simulation, we use natural units with the convention We concentrate here on the free wave case, i.e., we assume V (x) = 0 in Eq. (1).Our goal is to remove the wave function while approaching the boundaries, i.e., to apply an absorbing potential within a certain distance D such that either |x − x min | < D or |x max − x| < D.
We consider here as absorbing potential W , the local amplitude operator suggested by Kosloff et al. in [43]: , for x − x min < D W (x) = 0, everywhere else.
The "amplitude reduction" technique described [43] consists in applying the usual evolution operator to propagate from |Ψ(t)⟩ to |Ψ(t + ∆t)⟩ and consecutively applying on |Ψ(t + ∆t)⟩ an "absorption" step according to This method is equivalent to implementing an absorbing potential at the boundaries using a first-order Trotter-Suzuki decomposition.This justifies the use of the amplitude operator W proposed in [43] in our simulation, even though we do not use amplitude reduction technique.
To solve the problem on a classical or quantum computer, we directly discretize the one-dimensional space on a grid of equally spaced mesh points.The wave function |Ψ⟩ is then written as: The space grid {x i = x min + i.∆x} i=0,Nx−1 ranges from x min to x max , with a mesh step given by ∆x = (x max − x min )/(N x − 1).On the discretized mesh, the CAP becomes where k = D/∆x is a fixed integer.As a reference calculation, we first solve the problem on a classical computer using the split-operator method [44,45], i.e., going back-and-forth from position to momentum space.Consistently with the quantum computation case (see below), we used the simplest first-order splitting.In the classical computer case, the absorbing potential is directly applied to the wave function by multiplying each component ψ(x i ) by e −∆tW (xi) .In this case, the norm of the wave function decreases in time, as shown in Fig. 2 directly probing the effect of the absorption at the two boundaries of the mesh.Such classical computing results are rather standard and will serve as reference calculations for the one obtained using the quantum computing algorithm.
We show in Fig. 3  The parameters U0 = 0.4 and α = 1.5 have been used for the absorption potential.In the quantum simulation case, for each time step, the amplitudes are reconstructed from the measurement of the register system qubit using 2 14 shots.For such a large number of shots and in the absence of device noise, the error bars are essentially zero.Note that the y-axis is zoomed in panels (c-f).
sults have been simulated using an initially localized state located in the middle of the mesh with a width σ = 0.4 and N x = 2 4 = 16 points.The CAP parameters are U 0 = 0.4 and α = 1.5.For figure 4, the boost parameter value is set to v = 4.In these figures, we display both the results of the calculations with and without absorption.In the absence of absorption, we see the accumulation of the wave function amplitude that is reflected by the boundaries and then starts to interfere with the incoming wave function.As it is well known from classical computing, such interference prevents from properly describing the wave function spreading.This interference is reduced when the CAP is included, although some interference is still visible, especially in the boosted case.Note that such interference can be reduced by fine-tuning the absorption potential, but this is not the aim of the present work.Our main objective is to give proof that the same evolution as obtained in a classical computer can be implemented on quantum computers using the dilation method discussed above.

1D free wave evolution with CAP on a quantum computer
To encode the problem onto the qubit register, we use the standard binary (SB) representation that consists of mapping each position |x i ⟩ into the register state |[i]⟩ where [i] denotes the binary representation of the integer i.
The first-order Trotter-Suzuki method is used to perform the time evolution.Note that at first order, Trotter-Suzuki decomposition is strictly equivalent for the evolution to the "amplitude reduction" technique implemented on the classical computer discussed previously.As in the classical computer simulation, we perform the evolution of the Eq. ( 1) by going back and forth from the position space where the operators (V, W ) are diagonal, and the associated propagator can be easily implemented to the momentum space, where the same holds for the kinetic term K. Changing from one basis to the other requires the use of the Quantum Fourier Transform (QFT) algorithm or its inverse (see appendix A).This gives the scheme depicted schematically in Fig. 5.
The absorbing potential is included using the dilation method introduced in section 2.2.The blue box shown in Fig. 5 is repeated N t times, each time the result of the ancilla measurement being stored in a classical register, and the total system is measured only at the end.We keep only runs for which all the Figure 5: Schematic illustration of one time-step evolution including the implementation of the CAP using the reservoir qubit as discussed in section 2.2 together with the first order Trotter-Suzuki method where the QFT is used to go back and forth from momentum to position space.After a certain number of time step iterations, the probability density is obtained by measuring the qubits used to describe the system.Note that in the specific schematic example shown here, there are only 3 such qubits.
measurements of the reservoir qubit lead to zero.The U matrix is implemented using the algorithm discussed in section 4. For a given time, the amplitude is reconstructed by measuring the system register.We used here 2 14 shots for each panels of Fig. 3 and 4. We see in both figures that wave function amplitudes obtained using the quantum algorithm perfectly match the ones obtained using the classical computer algorithm.Note that we have systematically compared the results obtained on classical and quantum processor units, including or not an initial boost and/or adding a local potential V (x), and always obtained a perfect matching in the results.

Success probability and absorption
The good agreement between the quantum and classical computer simulations is further confirmed in Fig. 6 where we compare the evolution of the wave function norm during the time evolution for the two cases shown in Fig. 3 and 4. In the classical computer, this norm can be directly computed by integrating the wave function over the grid.As discussed in detail in section 2.2.3 (see also Fig. 2), in the quantum calculation case, the system wave function is re-normalized to 1 after each measurement of the reservoir qubit.To compute the physical norm, i.e., the one corresponding to the non-absorbed particle, we can use its connection with the success probability given by Eq. ( 14).In practice, the success probability at a given time t = r∆t is given by the probability to only obtain 0s in the measurements of the reservoir qubit for time steps 1 to r.The norm shown in the quantum calculation is deduced, assuming a strict equality between the success probability and the norm of the wave function.We see that the two calculations are almost on top of each other.As a side remark, we also implemented the prescription of Ref. [20].As we discussed in the introduction, in this case, the success probability decays exponentially with the number of times  In each panel, the classical computer simulation (red filled circle) is compared to the one obtained on the quantum computer with the dilation method (blue cross -dashed line).In the latter case, we directly show the success probability (see Eq. ( 14)) that is equal to the wave function norm when the number of measurements tends to infinity.Here, 2 10 = 1024 measurements are used for each point obtained in the quantum simulation.The errorbars, displayed by the dashed area, correspond to the standard deviation obtained by repeating 20 runs, each with 2 10 events and by estimating the width around the mean value within the 20 runs.Note that using 2 14 events instead of 2 10 as is used in Fig. 3  the dilation is implemented.For comparison, without boost, the success probability is of the order of 1/2 r after the r th time step, which should be compared to the values reported in Fig. 6.

Wave propagation with a confining potential
In the previous section, we tested the methodology in the extreme situation where the wave function expands and will eventually disappear completely from the mesh.Accordingly, the probability of success in our case will rapidly tend to zero, rendering the approach rather difficult to use due to the high rejection rate.Most of the physical applications where the CAP is employed, however, will lead asymptotically to cases where part of the wave function remains trapped inside the numerical mesh.In contrast, the remaining part is emitted from the localized source.This is ,for instance, the case when a system is excited and emits particles.Then, the fraction of emitted particles will depend on the internal excitation and particle emission threshold.In these physical situations, an essential aspect of the method we process is that the total success probability will not asymptotically tend to zero but saturate to the fractions of particles remaining in the numerical mesh.This is illustrated here by adding a confining potential.
Specifically, in addition to the CAP, we now add a Gaussian fixed potential centered in the integration domain and given by with V 0 = −1 and σ V = 1.The evolution of the norm of the wave-function -which is the success probability of measuring the reservoir qubit in state |0⟩, see Eq. ( 12) -is displayed in Figure 7, over 100 time steps.We can clearly see on this plot the advantage of our method, which is to lead to asymptotically non-vanishing total probability of success since it is directly related to the probability that the system remains in the box.

Complexity analysis for the implementation of the dilation matrix
The dilation algorithm we use to implement nonunitary propagation is expected to be one of the most efficient algorithms in terms of ancillary qubits since it requires the addition of only one reservoir qubit.Still, it requires implementing a general matrix U which is quite demanding (see, for instance, Ref. [46] for a comprehensive overview).We restrict the discussion here on the complexity in terms of gates for implementing this matrix directly without adding more ancillary qubits besides the one needed in the dilation.If n qubits are used to simulate the physical system, it leads to a matrix U of dimension 2 n+1 ×2 n+1 = 2 d ×2 d , whose naive implementation, as a general unitary matrix, requires in principle of the order of d 2 ×4 d gates [40].In [47], it was shown that, as an alternative to the dilation technique, the general complexity of implementing a matrix treating the imaginary-time propagation can be reduced by a factor 2 when using the singular value decomposition (SVD) technique.We show here that, due to the specific diagonal nature of the absorbing potential, and using the Cosine-Sine decomposition, it is possible to implement the dilation matrix U with much lower number of gates than the SVD-based approach.
To decompose arbitrary unitary operators of dimension 2 d × 2 d in terms of C-NOT and one-qubit gates, one of the most efficient decomposition scheme is the Quantum Shannon Decomposition (QSD) [46,48].This algorithm is based on the Cosine-Sine Decomposition (CSD) [49] and requires (3/4) × 4 d − (3/2) × 2 d C-NOT gates, where d is the total qubits number.If we use the CSD algorithm alone, without the improvements brought by the Shannon decomposition, we get a total complexity of 4 d − 2 d+1 C-NOT gates and 4 d elementary one-qubit gates [50].QSD and CSD algorithms are particularly well-suited to our case since the dilation matrix can immediately be written in terms of cosine-sine decomposition: where C and S matrices verify C 2 + S 2 = 1.We note that, due to the particular structure of the dilation matrix, it can always be written directly in this form.The l.h.s.matrix in (19) is a trivial multiplexor that can be implemented with one Z gate acting on qubit d.However, since we are only interested in the diagonal part of the dilation matrix to implement the non-unitary evolution, we can even remove this Z gate and consider directly the non-trivial part of the decomposition, namely the matrix D = C S −S C1 .
Thanks to the particular form of the dilation matrix in the case of non-unitary propagation with absorbing potential (see Eq. ( 9)), we can even further reduce this complexity.Our matrices C and S are indeed diagonal, and we only need to resort to a small part of the CSD algorithm to implement U .These matrices can be written as with i = 0, . . . 2 n − 1.For a given time τ , the angles are related to the imaginary potential through: where we can restrict θ i angles to [0, π 2 ].The operator D is a uniformly controlled rotation about the y axis, also denoted as F n n+1 (R y ( ⃗ θ)) [50].It consists of n-fold controlled rotations of qubit n + 1 about the axis y, one R y rotation for each of the 2 n different classical values of the control qubits.The circuit representation of D = F n n+1 (R y ( ⃗ θ)) is displayed in Figure 8.In general, F n n+1 (R y ) is a product of 2 n two-level operators.In [50] is proposed an implementation of F n n+1 (R y ) in 2 n C-NOT and 2 n rotations acting on qubit n + 1.In the specific case we consider here, i.e., with diagonal CAP, the implementation of the dilation method thus requires 2 n C-NOT and 2 n one-qubit gates.For comparison, in the same situation considered here, the SVD-based matrix requires solely to implement the diagonal matrix Σ introduced in [24] of size d × d and will require 2 n+1 − 1 z-rotation unary gate and 2 n+1 − 2 CNOT gates.We compare in table 1 the number of required CNOT for some optimized methods as a function of the system size register n.We see that the dilation method outperforms the other techniques for the specific implementation of local complex potential.For instance, the numerical simulation considered in Fig. 3 and 4 have been made using n = 4 qubits, for which a direct implementation of the U matrix requires only 16 CNOT and 16 R y operations per time-step.

Conclusion
We analyze the possibility of performing wave function Schrödinger evolution on a grid where the evolution mixes both real-and imaginary-time evolution.The physical motivation behind adding an imaginary potential is the possibility to absorb particles that might be emitted from a localized source without the need to treat very extended space regions.This technique is standard in classical computing and might potentially lead to a significant reduction in the number of qubits when the grid itself is encoded on the qubit register.
For the real-time evolution, we use the standard first-order Trotter-Suzuki method, where the system is evolved alternatively in position and momentum The two first lines are adapted from table I of Ref. [46] and concerns general unitary matrices.The two last lines consider the specific case where the CAP using local complex potential is implemented using either the SVD-based method of [24], or the dilation method we propose to use here.
space.We show that absorbing boundary conditions can be efficiently implemented using the dilation method.In particular, we propose a specific prescription for the ingredient of the dilation method that ensures that the success probability remains significant in most situations due to its connection with the physical wave function norm.Proof of the principle of the technique is given in the case of a one-dimensional particle evolving freely on a mesh and being continuously absorbed in the two mesh sides.We show that results obtained on a quantum computer are identical to those obtained on a quantum computer and that, due to the local nature of the absorbing potential usually used in applications.Finally, it is demonstrated that the dilation matrix can be implemented efficiently, significantly reducing the required quantum resources.
The solution of Schrödinger equation with absorbing boundary conditions is an archetype of an example where real-time and imaginary time propagation are mixed.Developing a specific strategy to implement such a mixing serves us as a test bench to test efficient quantum algorithms.However, when the number of available qubits will increase the accessible dimension of meshes might be large enough so that imaginary potential might not be necessary.Nevertheless, we would like to mention that this exploratory study and the resulting algorithms might be exported to more complex problems in the near future, like configuration interaction techniques involving complex energy due to the presence of resonances [52] that are very hard to solve in large Hilbert spaces on classical computers.

Acknowledgments
This project has received financial support from the CNRS through the 80Prime program and the AIQI-IN2P3 project.J. Zhang is funded by the joint doctoral programme of Université Paris-Saclay and the Chinese Scholarship Council.This work is part of HQI initiative (www.hqi.fr) and is supported by France 2030 under the French National Research Agency award number "ANR-22-PNQC-0002".We acknowledge the use of IBM Q cloud as well as the use of the Qiskit software package [42] for perform-ing the quantum simulations.The Qiskit code used for the quantum simulation presented in this work is available on request from the authors.

A Practical aspects of the kinetic term implementation
For completeness, we give below some details regarding the implementation of the kinetic term in momentum space (see also [28]).Since K is diagonal in momentum representation, we use the Fourier transform to go from space to momentum representation, apply Once in Fourier space, we must implement the diagonal operator e −ip 2 ∆t/2m in the basis |p⟩.Our qubit basis is defined in terms of N = 2 n discretized space points, n being the number of qubits considered, such that x k = x min + k.∆x where x min ≤ x ≤ x max , ∆x = xmax−xmin 2 n −1 is the discretization step in position space.The discretized momentum p can be written as : p = ( 2π L ) n−1 j=0 2 j p j , with p j = {0, 1} and p ∈ [0; 2π L (2 n −1)].However, in this case, we have no negative values for the momentum.To adapt our Fourier transform to a wave packet with both negative and positive momenta, we shift the p range to p = 2π L n−1 j=0 2 j p j − 2 n−1 .Then, for the range of momentum, we have: Then, we deduce: We dropped the last term since it corresponds to a global phase on the qubits.To implement this evolution, we thus need only phase gates and controlledphase gates.An inverse Fourier transform brings us back to the direct space, where we can implement the part of the evolution operator containing the potential.
Finally, we mention that the above treatment of the momentum space leads to the kinetic term spectral norm (see the interval given in ( 22)):

Figure 2 :
Figure 2: Schematic view of the norm evolution when solving the Schrödinger equation given by Eq. (1) on a classical computer (blue solid line).The horizontal black line and vertical gray dotted lines are added as guidance.The red dashed line corresponds to the evolution of the norm of the system state encoded on the qubit register when the problem is solved using the dilation method with the prescription M = e −W ∆t .At each time step ∆t, when the reservoir qubit is measured in state |0⟩, the system state is renormalized to have a norm equal to 1.

Figure 3 :
Figure 3: Time evolution of a wave-packet with no boost.The initial wave packet is shown in panel (a).Its half-width is chosen to be σ = 0.4.Snapshots of the wave function amplitudes at time t = r∆t where r = 1, 2, • • • , 5 are shown in panel (b) to (f) respectively with a time step ∆t = 1.2.The parameters U0 = 0.4 and α = 1.5 have been used for the absorption potential.In the quantum simulation case, for each time step, the amplitudes are reconstructed from the measurement of the register system qubit using 214 shots.For such a large number of shots and in the absence of device noise, the error bars are essentially zero.Note that the y-axis is zoomed in panels (c-f).

Figure 4 :
Figure 4: Same as Fig. 3 but with a boost corresponding to v = 4.

Figure 6 :
Figure6: Evolution of the norm of the wave function decreasing in time due to the CAP for the case of a freely expanding wave function without (a) and with boost (b).In each panel, the classical computer simulation (red filled circle) is compared to the one obtained on the quantum computer with the dilation method (blue cross -dashed line).In the latter case, we directly show the success probability (see Eq. (14)) that is equal to the wave function norm when the number of measurements tends to infinity.Here, 2 10 = 1024 measurements are used for each point obtained in the quantum simulation.The errorbars, displayed by the dashed area, correspond to the standard deviation obtained by repeating 20 runs, each with 2 10 events and by estimating the width around the mean value within the 20 runs.Note that using 2 14 events instead of 2 10 as is used in Fig.3and 4 leads essentially to vanishing error bars.
Figure6: Evolution of the norm of the wave function decreasing in time due to the CAP for the case of a freely expanding wave function without (a) and with boost (b).In each panel, the classical computer simulation (red filled circle) is compared to the one obtained on the quantum computer with the dilation method (blue cross -dashed line).In the latter case, we directly show the success probability (see Eq. (14)) that is equal to the wave function norm when the number of measurements tends to infinity.Here, 2 10 = 1024 measurements are used for each point obtained in the quantum simulation.The errorbars, displayed by the dashed area, correspond to the standard deviation obtained by repeating 20 runs, each with 2 10 events and by estimating the width around the mean value within the 20 runs.Note that using 2 14 events instead of 2 10 as is used in Fig.3and 4 leads essentially to vanishing error bars.

Figure 7 :
Figure 7: Evolution of the norm of the wave function in the presence of a Gaussian potential.The shaded area represents the statistical error bars.The norm is constructed directly from the success probability and we used 2 1 4 shots.

Figure 8 :
Figure 8: Definition of an n-fold uniformly controlled rotation on qubit n + 1, F n n+1 (Ry( ⃗ θ) (see text).The black control bits stand for value 1 and the white for 0. This operator match the D matrix we use to implement the non-unitary propagation.

F 2 3 (
R y ( ⃗ θ)), as proposed in Ref. [50], is shown in Figure 10.The rotation angles α used in the efficient circuit are linked to the original rotation angles θ by a linear transformation.
p 2 2m ∆t , and then revert back to position representation.If F denotes the Fourier transform operator, we have:e −iK∆t = F † e −i p 2 2m ∆t F.So to apply the operator e −iK∆t , we first perform a standard QFT, whose general expression for a standard basis state |x⟩ is[40]

Table 1 :
Illustration of the number of CNOT required to implement a unitary matrix of dimension d × d with d = (n + 1).