Exact boundary conditions at finite distance for the timedependent Schrödinger equation
Abstract
Exact boundary conditions at finite distance for the solutions of the timedependent Schrodinger equation are derived. A numerical scheme based on CrankNicholson method is proposed to illustrate its applicability in several examples.
pacs:
: 03.65.Bz, 03.65.Nk, 03.80.+rI Introduction
The timedependent quantum mechanics formalism has been extensively used in the recent years in several problems of physics and chemistry. Since the work of McCullough and Wyatt [1] showing the interest of dealing directly with nonstationary states in calculating reactions probabilities, a consequent number of works has been devoted to the wavepacket dynamics. The reader could find the more recent developments in the reviews [2, 3, 4, 5] and references therein.
When solving the timedependent Schrödinger equation in a finite volume, one is faced with the problem of implementing the proper boundary conditions at finite distance. These boundary conditions are usually those of a free wave outgoing propagation. For the stationary solutions they are given by fixing the value of the logarithmic normal derivative on the boundaries of the interaction domain. However for the case of nonstationary wave packets a satisfactory solution is not known.
The usual way to overcome this difficulty is to impose the nullity of the wave on the boundaries and push them far enough to avoid the effect of parasite reflections in the region of interest [6]. Another way, widely used in other branches of physics and mechanics as well, is by means of the socalled absorbing boundary conditions [7, 8, 9]. They indeed minimize parasit reflections but perturbe the dynamics near the boundary and do not provide a solution of the initial equation in its neighborhood.
A third way is by the so called wave function splitting algorithm which consists in removing from the total wavefunction the part localized outside the interaction region before it reaches the boundary [10]. Other possibilities based on the interaction picture formalism have also been developed [11].
Although the above quoted methods provide a practical way to solve the problem they do not give a theoretical solution of it. The aim of this paper is to formulate an exact boundary condition (EBC) at finite distance for outgoing solutions of the timedependent Schrödinger equation. Consequently, it would allow a solution of timedependent quantum mechanical problems in a finite spatial domain without any approximation in the boundaries. This theoretical result is completed by implementing the derived boundary conditions in some practical calculations concerning the more usual problems of the wavepacket dynamics.
The paper is organized as follows: in Section 2 we derive the general expression for the EBC which has to be imposed at the boundary of an arbitrary integration domain containing the interaction. This conditions turns to be non local on time and allows a solution limited to the interaction domain of the timedependent Schrödinger equation. In Section 3 the boundary conditions are reformulated for a discretetime evolution problem in the frame of the CrankNicholson scheme. Section 4 is devoted to the one dimensional Schrödinger equation for which the use of these conditions is particularly simple. Several illustrative examples are given in Section 5: the spreading of a free wave packet, the scattering of a wavepacket on a static potential, the behavior of an initially localized wave packet under a timedependent perturbation. In the last example we derive explicit solutions of the Schrödinger equation in a timedependent delta potential. The applicability of the formulated conditions in a twodimensional problem is shown in Section 6. A discussion on the interest and limitations of the EBC is finally presented in the conclusion.
Ii Boundary conditions at finite distance for the timedependent Schrödinger equation
Let us consider the dimensionless timedependent Schrödinger equation in a ddimensional configuration space
(1) 
The potential , eventually timedependent, is assumed to vanish outside a finite domain named the interaction domain.
We wish to obtain a solution of (1) by solving the equation in a finite spatial region , called hereafter the integration domain. The domain has to be chosen such that it contains , i.e. such that its complementary with respect to the total configuration space is free of interaction. We have thus and ( see figure 1).
The solution of (1) in requires to formulate and implement the proper boundary conditions on its border, a surface . For this purpose let us first consider the solution of (1) in . The domain is, by definition, free of interaction and satisfies the free Schrödinger equation
(2) 
Let be the free retarded propagator, i.e. a solution of
(3) 
with and arbitrarily chosen.
We have in mind the outgoing solution at of a timedependent problem with initial value given by for . By multiplying equation (2) by , equation (3) by and substracting one gets the relation
(4)  
(5) 
with .
This relation is to be integrated over the dimensional spacetime volume . We will assume that the domain is connected. This assumption is justified for most of the cases except for , which will be discussed later. The integral of the right hand side is then transformed, by using the Green theorem, into a flux term through the boundary surface . We denote by the normal to pointing towards (see figure 1).
Taking into account that for t and that there is no contribution from the flux at the infinity for purely outgoing waves (for they have the same asymptotic behavior as the retarded free propagator) we obtain the following relations depending on the relative position of .
If one has
(6)  
The solution in the free region is thus obtained as a sum of a first term, representing the free propagation of the wave which was initially in the non interacting region , and a second term made only with the values of the wavefunction and its normal derivative at the boundary at the past times. Equation (6) can be consequently used to propagate outside a solution obtained in and can be viewed as a generalization of the Huyghens principle in a spacetime surface.
If one has
(7)  
This result can be obtained either as a limiting case of the preceding one (6) or by directly integrating equation (4) for . Equation (7) is the boundary condition we were looking for. It provides an integral relation between the value of the wavefunction on a point of at a given time and the values of the wavefunction and their derivatives in at the preceding times. Starting from a known initial state at time , equation (7) gives in its differential form the boundary condition on required for determining at time the solution of the timedependent Schrodinger equation in the finite domain .
Several comments are in order:

The relation (7) is non local in both space and time. It generalizes the outgoing boundary condition for the monochromatic free plane waves, which in the onedimensional case and at reads:
(8) 
The choice of ensures that only free outgoing waves will propagate across , in the same way as a solution obtained in an infinite domain. This result is exact and consequently will not generate any reflection on the boundary.

The first term of its right hand side corresponds to the free propagation of the initial state which was initially in the external region . Its numerical calculation needs some care due to the oscillatory character of the propagator and several ways to overcome it are proposed below:

It is always possible to eliminate this contribution by choosing the initial time and in such a way that the initial state is fully in .

In the case when the initial state is very extended, it is more suitable to include it entirely in . Indeed the volume integral term gives then the value at the boundary of a free evolution of the initial state. This can be easily calculated in the momentum space and is known analytically for some cases of physical interest like gaussian wavepackets.

In the case when the initial state is entirely in , there exists another possibility to remove the volume term by splitting the solution into its non perturbed and scattered parts, i.e.
in which is the solution of (2) which coincides at with the initial state and a solution of the inhomogeneous equation
(9) The unknown function vanishes at and satisfies in the free domain the same equations than . It will consequently satisfy on the boundary conditions (7) without the volume term.

The preceding relations, completed by the case , are summarized in the following equation
(10)  
with defined by
The above derivation has been done under the assumption that the potential vanishes outside a finite domain . This condition can be however weakened in some cases. For instance in the case of an interaction made of one short range part plus a static long range part (e.g. Coulomb, polarisation or centrifugal potentials) the same derivation holds. The interaction domain is then defined by whereas can contain all across. The only difference in equation (10) consists in replacing the free propagator by the corresponding propagator of the long range potential .
In the one dimensional case the interaction domain is and the free exterior domains has two connected components and . The boundary condition (10), based on the Green theorem, has to be applied to each of them. Let us take for instance . At , the normal derivative of the free propagator on the boundary vanishes, as it can be explicitely checked on its expression ():
(11) 
The EBC (7) becomes simply
(12) 
It is worth noticing that in case this result may be obtained in a quite straightforward way. Indeed let us consider the Fourier energy components of the wavefunction
(13) 
and let us impose to each of them the plane wave outgoing boundary condition (8) on the form
Inserting this relation in (13) one gets
(14) 
We recognize the Fourier transform of the free Green function for zero argument
(15) 
and we obtain the following relation:
(16) 
in agreement with the second term in the right hand side of (12). The first term does not contribute in the limit due to the fact that the free propagator tends to zero.
The boundary condition (10) generalizes the KKR [12] method which provides a relation between the solution and its normal derivative at the boundary of a free domain. This method has been derived and successfully applied in the case of static quantum billards. The formulation we have developed in this section appears as a natural extension of this method for spacetime billards.
The condition (10) is not directly useful in a numerical solution of the problem. Indeed the singularity of the propagator at would generate instabilities when the dependent Schrödinger equation is discretized. To overcome this difficulty one has to reformulate the boundary conditions in the particular approximate scheme chosen. The aim of the following section is to obtain such a formulation in the frame of the CrankNicholson scheme.
Iii Boundary conditions in the CrankNicholson scheme
The CranckNicholson [13] method for solving a discretetime evolution problem consists in approximating the infinitesimal evolution operator
by the expression:
(17) 
where and are the wavefunctions at two consecutive times separated by and the hamiltonian H is evaluated at the mean time .
By introducing the imaginary parameter , equation (17) can be rewritten in the form
(18) 
according to what the solution at timestep m is obtained by solving an inhomogeneous stationary Schrödinger equation with complex energy . The source term is provided by the solution at the preceding time.
By means of (18) the equivalent of equation (2) for the outgoing wavefunction in the free domain becomes
(19) 
and the discrete free retarded propagator is defined as the solution of
(20) 
normalised to
For one obtains the analogous to the continuum case by simply replacing the time integral by a sum over the discrete index p
In calculating the limit and like in the continuum case there is a singularity in the normal derivative of the free propagator. Its contribution to the surface term in the r.h.s of (III) is
for and the opposite for . The volume contribution can be in its turn rewritten by means of the regularized propagator defined as
(22) 
and one finally arrives to
This expression is the equivalent of (7) for the discrete time evolution.
To make use of the discrete EBC condition (III) it remains to derive an expression for the discretetime free propagator . In the momentum space equation (20) reads as
which leads to:
(24) 
One can see in equation (24) that in the limit the Fourier transform of the propagator contains the singular delta function manifested in (22). Nevertheless this term cancels in the sum of two successive propagators. This sum is given by the non singular integral:
(25) 
This integral can in principle be calculated by the contour Cauchy method. However, its value can also be obtained by means of the discrete Fourier tranform. We note indeed that a N discretetime Fourier transfom with respect to the index p of equation (25) gives the sum of a geometric series. With the inverse Fourier transform we obtain:
(26) 
where and is the required convergence factor of the geometric series which has to be chosen such that . is the standard free Green function with energy .
Although the formalism has been developed in view of solving the Schrodinger equation in a finite domain , equation (III) with provides the wavefunction in all the configuration space. However one can obtain a closed form for some observables in terms of the internal solution only. Let us consider for example the overlapping between two outgoing solutions et of (19) in the exterior domain. Standard algebra gives for the exterior part of the integral the expression
(27) 
which can be calculated from the solution in . By adding the corresponding equations for consecutive time steps one obtains the desirate scalar product. This expression is used e.g. in calculating the autocorrelation function (see section IV.5). The case provides the rate of variation of the probability in the exterior domain by time step. The contribution from can be calculated by numerical integration and we thus have the possibility to verify the conservation of the total norm .
In summary the solution of the time dependent Schrödinger equation in the interaction domain can be obtained by solving, at each time step, a stationary complex and inhomogeneous Schrödinger equation (17) in with the boundary conditions given by (III). The local character of this equation is preserved by these conditions. We are thus let with the solution of a banded linear system, whatever the discretisation algorithm used for the spatial variables. For static hamiltonians this linear system is in addition the same at any timestep.
Iv Examples in the one dimensional case
In the one dimensional problems, simplifications arise. Like in the continuum case, the normal derivatives of the propagator at the boundaries vanish. The sum of two successive propagators for can be explicitly calculated and gives:
(28) 
This integral can be performed by standard methods: it vanishes for odd values of p and for it gives
(29) 
with
The boundary conditions (III) give at :
(30) 
At one obtains in a similar way
(31) 
For illustration we have detailed some simple examples concerning the time evolution of wavepackets. The initial state has been taken for simplicity of gaussian form
(32) 
The examples have been solved by scaling the integration domain to the interval . The discretetime Schrödinger equation (17) is there solved at each of the N timesteps by finite differences method with equally spaced grid points . The values of the wave function at are the unknowns. The required equations are obtained by validating (17) at each grid point. However to evaluate the second derivative at the end of the grid the values of the wave function at the nearest points exterior of this domain is needed. These are written in terms of the internal points by using EBC conditions (30,31) with a symmetric formula to evaluate the derivative term. Substituting this exterior values in (17) we end with a tridiagonal linear system of dimension .
iv.1 Free propagation of a gaussian wave packet
The first example concerns the free propagation of a gaussian wavepacket. Its main interest lies in the fact that an analytical solution is known which allows to check the validity of our boundary conditions. In this case the integration domain is not fixed by the interaction but chosen such that it fully contains the initial wavepacket.
The wave packet (32) is initially centered () and, in absence of interaction, the probability density at time is known to be given by :
with and a dimensionless time variable. If the total evolution time in this units is , one has . . We have taken for the initial width the value
In figure (2) are displayed the density probabilities each 4 time steps. The time is equal to , the number of time steps is and . The wave packet velocity is chosen such that during the time a classical particle would travel half the interval . The initial state is shown by dashed line. One can see by simple inspection the absence of any parasite reflection nor anomalous behavior. The curve corresponding to the time is compared to the analytical result. The difference is not visible by eyes except at the vicinity of in which it reaches its maximum value, less than 1%. This difference has been arbitrarily reduced by increasing the number of time and space grid points showing the validity of this approach. The parameters chosen in this simple example may serve to illustrate the efficiency and accuracy of this method.
In the example considered above and due to the value of , the probability flux leaves the interaction domain mostly at . However, due to the spreading of the wavepacket, the wave function at is also being populated although its value remains very small and is not visible in the figure. By taking v=0 one has a symmetric situation for which the same quality of results has been obtained.
A last check has been done concerning the total probability conservation. The contribution coming from the interior domain has been obtained by integrating (trapezoidal rule) the calculated wavefunction. The contribution from the external part has been evaluated by equation (27) which becomes in that case
(33) 
In this example the probability is found to be conserved better than .
iv.2 Scattering by a static potential
In the second example, the initial wave packet (32) scatters on an attractive potential:
with and . The wavepacket parameters are , , .
The probability densities at different times are displayed in figure (3). Its initial value is plotted in dashed line and the attractive well in bold. The reflected and transmitted wave packets are clearly separated. The integrated probabilites in and have been calculated respectively by (33) and trapezoidal rule from the solution in . They are shown as a function of time in figure (4). The contribution from the and regions are respectively plotted in curves (a) and (b). They correspond to the reflected and transmitted wavepackets. The difference between both curves corresponds to the integrated probability in the integration domain. It tends to zero very slowly.
iv.3 Localized state under a timedependent perturbation
The third example concerns the evolution of a wave packet initially localized on an attractive gaussian potential whose amplitude varies periodically:
(34) 
with , and .
The time of evolution is T=80, with 800 points in space and 800 timesteps. The initial state is a superposition of the bound states plus a small component () of continuum states both in the corresponding static potential at . Curves (a) and (b) on figure (5) represent respectively the probability density in the intervals and . The difference between these curves, i.e. the probability density in the integration domain, tends to zero showing that the bound state is pulled out of the potential by the time dependent perturbation. The evolution in the corresponding static potential, shown in dashed lines, tends to a constant: the projection of the initial wavepacket into bound states.
iv.4 Tunneling
The following example illustrates the tunnel effect through a double well. A gaussian wavepacket with , initially centered and at rest evolves in a double repulsive well of the form
(35) 
with and . We have displayed in figure (6) the probability density at different times. The initial state is drawn in dashed line and the potential (arbitrary units) in bold. The wavefunction spreads and oscillates inside the potential except for a small part that tunnels through. The time evolution of the integrated probability densities corresponding to figure (6) is shown in figure (7). Curve (a) corresponds to the domain and curve (b) to . The distance between both curves represents the integrated probability in the integration domain . The small oscillations correspond to the back and forth reflections of the wavepacket inside the well.
iv.5 Timedependent Delta potential
An interesting example is provided by the limiting case of a timedependent delta potential
In this case the interaction domain is reduced to a point and it is possible to obtain closed analytical solutions for the discretetime evolution problems.
We will consider the evolution of a state which initially coincides with the bound state of the potential with strength . Its energy, in current units (), is and the normalized wavefunction
By integrating the discretetime Schrödinger equation (18) on both sides of the singularity one gets the discontinuity of the wave function derivative at
(36) 
where is the mean intensity of the potential between the and timesteps.
The solution of (18) will be written by splitting its non perturbed and perturbed parts
(37) 
By definition, is a stationary solution of equation (18) with . In the continuum case the time evolution of would be simply given by
In the discretetime evolution it is easy to show that
with
The perturbated part is a solution of (19) on both sides of and satisfies the boundary conditions (30). The sum of these two conditions gives
(38) 
Combining this last result with equations (36) and (37) one obtains at x=0
(39) 
This is an inhomogeneous linear system for the perturbed wave function with the inhomogeneous term given by . By isolating the q=0 contribution equation (39) results into an explicit recurrence relation
(40) 
The perturbed wave function at the origin can be thus easily obtained. The values at can be calculated from equations (III) and (26) This example illustrates how we can extract all the physical quantities from the knowledge of the wave function and its derivative on the boundary, even in a limiting case when it is reduced to a point.
In figure (8) are displayed the results corresponding to a periodically driven delta potential with
and where is the Bohr frequency of the bound state for . Curve (a) shows the square modulus of the autocorrelation function and curve (b) the square modulus of the wave function at the origin normalized to one at the initial time. For the long times the wave function near the origin becomes proportional to the bound state and these two quantites are close to each other. There are 1000 time steps for 40 perturbative pulses. We note that after a transient period, the wave function is near of an eigenvector of the Floquet operator [14] with complex energy.
The periodically driven delta potential is used as a model of weakly bound atom (1D delta atom) when studying the tunnel effect in a static exterior field [15]. In this case the exact propagator is known explicitly. In our case there is also an explicit solution.
V A two dimensional example
We present in this section a straightforward application of the formulated boundary conditions for a two dimensional problem. We will illustrate it by considering the free evolution of a two dimensional gaussian wave packet crossing obliquely the axis. The wave packet has initial widths and is initially centered at the point of the XY plane with velocity . We push to the infinity the domain for the y variable and consider as integration domain the band
Let us apply the boundary conditions (III) to a point in the line . The integral term has only contributions coming from itself and is thus reduced to a one dimensional integral although with a two dimensional propagator of nonvanishing argument .
(41) 
with and .
An efficient way to solve (41) is by Fourier transforming the wavefunction with respect to the variable
and impose the desired boundary condition to each of its Fourier components . The problem is then formally equivalent to a one dimensional one with a propagator given by expression (28) where the variable is replaced by . The resulting integral is no longer analytic as in the one dimensional case but can be easily evaluated by means of a discrete Fourier tranform according to equation (26).
The numerical results are displayed in figure 9. They correspond to the solution of the discrete time Schrödinger equation in the region and with time steps. The number of points in x direction is and points were required for an accurate Fourier transform of amplitudes. The countour plots of the probability density are shown for six equidistant times.
The figure shows the spreading of the wave packet in both directions and the crossing of the boundary without any parasit reflections. The comparison of calculated values with the exact analytic expressions showed a relative accuracy at the maximum of the wave packet better than .
Vi Conclusion
Exact boundary conditions at finite distance for the timedependent Schrödinger equation have been derived. These conditions allow a resolution of timedependent quantum mechanical problems in a finite spatial domain containing the interaction without introducing any parasite reflections, nor changing the dynamics near the boundary. They are specially useful to observe long time evolution of waves packets in a finite space domain.
The value of the wavefunction outside the integration domain can be easily obtained from the values of the wave function and its normal derivative at the boundaries for all the anterior times.
These conditions have been explicitely reformulated in a timediscrete dynamics and implemented in a numerical algorithm providing a numerical solution of the more usual onedimensional problems. They result in solving at each discrete time step a stationary Schrödinger equation with complex energy in a finite spatial domain. The local character of this equation is preserved by the boundary conditions and we are thus let with a band matrix problem whatever the discretisation algorithm used.
The economy resulting when solving the problem in a reduced spatial domain allows to considerably increase the number of grid points and to use for a same accuracy less powerful but simpler algorithms. The numerical results obtained with the Cranck Nicholson scheme are found to be very satisfactory. Other time and space discretisation algorithms can however be used as well at the cost of the simplicity.
The derived boundary conditions can be applied to more general situations than those considered in the examples. The case of higher dimensionality and/or of coupled equations, for instance, can be treated without modifying the formalism. A numerical application of the free evolution of a twodimensional gaussian wave packet is given. Long range static potentials may be taken into account as far as an analytic solution is known for the corresponding propagators. The generalization to the timedependent treatment of the three body problem is straightforward and could be of some interest in some molecular physics problems [17] or a doorway to solve the still open question of Coulomb problem above threshold.
The interest in formulating exact boundary conditions at finite distance for timedependent problems recently appeared in the field of statistical mechanics and conclusions similar to those presented in this work have been reached in [18]. But this formulation is also relevant in other branches of physics than those governed by the Schrödinger equation. In particular for the diffusion equation, for which the same results hold with imaginary time. The case of hyperbolic equation of interest in many fields of physics can be treated in the same footing [16].
Vii Acknowledgments
We acknowlegde Dr. K. Protassov for his support on some mathematical aspects of this work and Prof. C. Leforestier for helpful discussions.
References
 [1] E.A. McCullough, R.E. Wyatt, J. Chem. Phys. 54, 3578 (1971)
 [2] TimeDependent Quantum Molecular Dynamics,NATO ASI Series 299, Ed. by J. Broeckhove et al., Plenum Press 1992
 [3] M. Kleber, Phys. Rep. 236, 331 (1994)
 [4] B.M. Garraway and K.A. Suominen, Rep. Prog. Phys 58, 365 (1995)
 [5] N. Balakrishnan, C. Kalyanaraman, N. Sathyamurthy, Phys. Rep. 280, 79 (1997)
 [6] A. Goldberger, H.M. Schey, J.L. Scwwartz, Am. J. of Phys. 35 177 (1967)
 [7] C. Leforestier and R.E. Wyatt , J. Chem. Phys. 78, 2334 (1983)
 [8] R. Kosloff, D. Kosloff, J. of Comput. Phys. 63, 363 (1986)
 [9] A. Vibok, G.G. BalinKurti, J. Chem. Phys. 96 7615 (1992)
 [10] R. Heather and H. Metiu, J. Chem. Phys. 86, 5009 (1987); 88, 5496 (1988); R. Heather, Comput. Phys. Com. 63, 446 (1991)
 [11] J.Z.H. Zhang, Chem. Phys. Lett. 160, 417 (1989)

[12]
R. Balian and C. Bloch, Ann. Phys. 60, 401 (1970)
M. Berry, Ann. Phys. 131, 163 (1981)  [13] Numerical Recipes W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery Cambridge Univ. Press New York 1986
 [14] S. H. Shirley, Phys. Rev. 138, 3979 (1965)
 [15] S. Geltman J. Phys. B 11, 3323 (1978)
 [16] O. Meplan and C.Gignoux Phys. Rev. Lett. 76 408 (1996)
 [17] C. Leforestier, F. LeQuéré, K. Yamashita, K. Morokuma, J. Chem. Phys. 101, 3806 (1994)
 [18] J. R. Hellums and W. R. Frensley, Phys. Rev. B 49, 2904 (1994)