Absorbing boundary conditions for dynamical many-body quantum systems
In numerical studies of the dynamics of unbound quantum mechanical systems, absorbing boundary conditions are frequently applied. Although this certainly provides a useful tool in facilitating the description of the system, its applications to systems consisting of more than one particle is problematic. This is due to the fact that all information about the system is lost upon absorption of one particle; a formalism based solely on the Scrhödinger equation is not able to describe the remainder of the system as particles are lost. Here we demonstrate how the dynamics of a quantum system with a given number of identical fermions may be described in a manner which allows for particle loss. A consistent formalism which incorporates the evolution of sub-systems with a reduced number of particles is constructed through the Lindblad equation. Specifically, the transition from an -particle system to an -particle system due to a complex absorbing potential is achieved by relating the Lindblad operators to annihilation operators. The method allows for a straight forward interpretation of how many constituent particles have left the system after interaction. We illustrate the formalism using one-dimensional two-particle model problems.
1 Introduction and basic theory
Large efforts have been invested in understanding the quantum mechanical behaviour of dynamical unbound systems involving several particles. Experimental advances allow us to study processes in which internal interactions between the constituent particles play a crucial role – in addition to dynamics induced by external perturbations. Examples of such processes may be photo ionization of helium  and cascaded Auger processes following ionization of inner core electrons in atoms . Of course, in order to understand such processes, theoretical and numerical studies are required. Furthermore, the theoretical interest of such systems, be it within the context of solid state, molecular, atomic or nuclear physics, is spurred by the fact that they represent demanding tasks. One obvious challenge in describing unbound systems is that their spatial extension may become arbitrarily large. Moreover, even if one is able to represent the whole system within a finite space, extracting relevant information from the final wave function may be far from trivial.
The unbound quantum systems under study may often be thought of as having an interaction region of finite spatial extension and an asymptotic region where the unbound part of the system travels outwards. As this asymptotic behaviour often is well known, it may be desirable to describe only the dynamics of the part of the wave function belonging to the interaction region. In a numerical implementation this cannot be done by simply resorting to a representation of space that is smaller than the extension of the wave function as this leads to unphysical reflections at the boundary. However, it can be achieved by imposing absorbing boundary conditions, i.e. by demanding that the wave function vanishes as the particle approaches the edge of the numerical grid – preferably with as little reflection as possible, see references [3, 4] or the recent review of reference . Such absorbers are frequently referred to as perfectly matched layers in the context of general wave equations .
When propagating a wave packet on a numerical grid, a common way to impose absorbing boundary conditions is by adding a complex absorbing potential (CAP) to the Hamiltonian of the system . Complex absorbing potentials are widely used e.g. within molecular dynamics [8, 9, 10] and atomic physics [11, 12]. Alternatively, this way of absorbing particles may be formulated as multiplying the wave function with a masking function at each time step .
In any case the resulting effective Hamiltonian acquires an anti-Hermitian part,
where both and are Hermitian and is positive semi-definite ().
The wave function of the system obeys the Schrödinger equation
A density operator correspondingly obeys the von Neumann equation
It is easy to see that the evolution is non-unitary, viz,
so that probability is “lost” if the density matrix overlaps with anti-Hermitian part of the effective Hamiltonian.
Although methods involving non-Hermitian effective Hamiltonians may come in very handy in reducing the complexity in describing potentially unbound systems, there is one major problem when the system contains more than one particle: As one particle leaves the rest of the system and is subsequently absorbed, the entire wave function is lost. This is obvious since the wave function is normalized to the probability of finding all particles within the space defined by the numerical implementation. In other words: If an initial -particle system gradually “loses” one particle through some non-Hermitian “interaction”, the corresponding wave function gradually goes to zero – not to some wave-function corresponding to particles.
As it is desirable to be able to describe dynamics where several particles are lost one by one, one may try and construct a formalism where an -particle wave function is created as one particle is lost, and where this wave function may be propagated using the corresponding Schrödinger equation including some source term. Once being able to do this, the extension to particles etc. follows by induction. However, as it turns out, such a construction is not possible due to the fact that the process of losing particles in this way is irreversible: As a particle is absorbed, information is irretrievably lost. Hence, a Markovian master equation should be a more suitable starting point than a pure state approach. Lindblad  was able to show that in order for the evolution of a system to be Markovian, trace-conserving and positive, the density operator has to obey an equation of the form
with the so-called Lindbladian given by the generic expression 
Here, is the Hamiltonian governing the unitary part of the evolution. The operators are referred to as Lindblad operators. In a diagonal representation, the Lindbladian simplifies to
Equation (5) is usually used to describe energy dissipation to the environment . However, it has recently been demonstrated that the spontaneous decay of unstable particles may also be described within this formalism [16, 17]. In reference  a master equation of Lindblad from is obtained for a Hilbert space consisting of unstable particle states and vacuum, however without any dynamical degrees of freedom. The decay rates of the unstable particles are introduced as parameters. Similarly, in reference  decay from one system to another one consisting of its decay products is described by introducing a single Lindblad operator. In these works it is demonstrated that one may describe decoherence along with decay in this context.
In the present work these ideas have been used to generalize the standard one-particle absorbing boundary technique to an -body setting, and it is demonstrated that this is indeed the natural way to do this. Alternatively, in a more general context, this formalism allows the study of a system which gradually loses particles to some environment due to the non-Hermitian part of .
We comment that the absorption process does not affect the dynamics on its interior. However, removing a particle may still affect the other particles via their interaction. Hence, in order to obtain physically meaningful results, any dependence of the absorber must be eliminated by placing the absorber sufficiently far away from the interaction region – as is the case for any implementation featuring absorbing boundary conditions.
In the following section, section 2, the formalism will be formulated for identical fermions exposed to a complex absorbing potential. In section 3 two examples featuring two identical fermions in one dimension are given. Finally, in section 4, conclusions are drawn and a few future perspectives are outlined.
2 Fock space description
The natural setting for describing a system of a variable number of fermions is Fock space, the direct sum of all -fermion spaces . As we wish to describe at most particles, it suffices to consider
An arbitrary -fermion state can be written
where the field operator creates a particle at position . The field operators obey the usual anti-commutator relations
Here, may denote all degrees of freedom associated with a particle, including eigenspin. Moreover, is the anti-symmetric local wave function of the -particle system, and is the vacuum state, containing zero particles.
A system of identical fermions interacting with at most two-body forces has the Hamiltonian
Here, (resp. ) is a one-body (two-body) operator acting on the degrees-of-freedom associated with particle (and ). Using field operators, the Hamiltonian can be written
the point here being that given on this form is independent of the number of particles in the system. The exact expression is irrelevant at this point. The CAP, which is diagonal in , is conveniently introduced as
in order to reproduce the anti Hermitian “interaction” leading to absorption. Here we have allowed for an integral instead of a sum in (7). Furthermore, the last term of the Lindbladian, which is absent in the von Neumann equation, , should map an -particle system into an -particle system. Hence, should map into , which means that the Lindblad operator must be on the form
The simplest choice that satisfies (14) is the diagonal form
We will return to the justification of why this is the proper way to define the Lindblad operators shortly.
and the master equation (5) becomes
This is our fundamental dynamical formulation of particle loss due to a CAP.
Let us consider the master equation (18) in some detail. We may partition the density matrix into blocks, viz,
where , with being the orthogonal projector onto . Each block evolves according to
showing that the flow of depends on that of , but not the other way around. This is illustrated in figure 1. Also, notice that obeys the original von Neumann equation (3) as there are no couplings to other blocks in this case as particles do not enter the -particle system.
We now return to the question of whether the definition (16) of the Lindblad operators, which obviously is the simplest one, is the only adequate one. Indeed, during the transition from an -particle system to an -particle system, only the unabsorbed part of the original system should be reproduced within the -particle system, and this leads to the above choice. To see this, consider a somewhat idealized example consisting of two non-interacting particles. Since they do not interact, their wave function is given by a Slater determinant at all times, . The state does not overlap with the absorber at any time, i.e.
and corresponds to an unbound particle travelling outwards. We suppose that as approaches infinity, the particle in the state is completely absorbed, and our numerical representation of the final system should converge towards the pure one-particle state .
Indeed, the evolution dictated by (2) follows this pattern. The block of the density matrix corresponding to the two-particle system, , remains a pure state – albeit with decreasing norm as is absorbed – and the evolution of the one-particle block simplifies due to (21) to
Hence, the one-particle part of the system is always proportional to a pure -state, and, since the trace of the entire system is conserved, this simply integrates to , as it should.
On the other hand, with a non-diagonal form of the Lindblad operators, c.f. (15), we would have contributions to of the form and its Hermitian adjoint in addition to the pure state contribution. These contributions clearly cannot be allowed, as should be independent of the state of the outgoing particle. Hence, we find that, up to an arbitrary phase factor, the definition (16) is the proper way to define .
For a non-interacting -fermion system where particles are bound and particles are unbound, and with an initial state given by the Slater determinant , it may be verified by inspection that the source term for the -th block, , is always proportional to given that none of the -states overlap with . The intermediate density matrix blocks , where , will approach zero as , but transiently describe (mixed) states where particles have left.
Of course, in the more interesting context of interacting particles, the structure of the diagonal blocks of the density operator, , is more complex than in the special case of non-interacting particles.
2.1 Consequences and interpretation
Typically, our starting point is a pure -particle state, , in which case (2) reduces to the ordinary Schrödinger equation (2), which was our original formulation. Moreover, it is easy to see that will remain block diagonal for all in this case, i.e. if . But (2) shows that in general is a mixed state, unlike . This is due to the information loss when admitting ignorance of the whereabouts of the removed particle.
We stress that by construction, for all . Probability flows monotonically from into , and the particle absorbed from is not present in , but is erased, and is a proper description of an -fermion system. In this way, we see that the above construction is a natural generalization of the original non-Hermitian dynamics, which describes the classical removal of a particle, into one that also describes the remaining system in a consistent way. It is striking to notice that the CAP, , is already given as one of the terms in the Lindbladian, so that the Fock space formulation follows almost immediately.
It is worthwhile to mention that the traces of the blocks along the diagonal of have very simple interpretations as the probability of having particles in the system upon a measurement, i.e.
In particular, , and . Hence, within this formalism, distinguishing between single, double etc. ionization of atoms and molecules is straightforward.
For any observable the expectation value is given by . For example, the expected number of particles in the system is given by
It may also be useful to consider conditional expectation values:
i.e., the expectation value of given that the system is found in an -particle state.
Obviously, as particles are “removed” by the absorber, information is lost. This information loss may be quantified by the von Neumann entropy, , or the closely related notion of purity, defined by 
Unity minus this quantity, , is a measure of the amount of mixedness, and the purity is 1 for pure states only. Similar to conditional expectation values, c.f. (25), one may define conditional purity and von Neumann entropy as the corresponding quantity of the re-normalized block, i.e. .
Of course, for a full quantum mechanical description in terms of the complete (unabsorbed) -particle wave function, no information is lost. Indeed, the absorption of particles is a semi-classical concept. In the full -fermion quantum system no such separation is possible. On the other hand, as the Schrödinger equation for the particles is separable in the non-interacting and states discussed above, giving a Slater determinant of the time evolved one-particle states as the full solution, the Lindblad equation is seen to correctly construct the -particle Slater determinant resulting from the removal of the outgoing particles from this Slater determinant. Thus, the Lindblad equation exactly encapsulates the approximate separation of non-interacting quantum systems far apart.
3 Example: Two identical spin -particles in one dimension
In the following we consider two fermions in one dimension with the one-body Hamiltonian given by
where is some external, possibly time-dependent, potential. The fermions interact via a potential and the extension of the system is effectively reduced by imposing a CAP. With two interacting identical fermions, we are confined to the Hilbert space
We discretize this system using a uniform grid in the interval containing points , , with . The field operators can be replaced by a finite number of creation operators , creating a particle at grid point . These operators, which obey the usual fermion anti-commutator rules
map discrete anti-symmetrized -function bases for into bases for .
The Hamiltonian takes the form
where are the matrix elements of the one-body Hamiltonian, which depend on the chosen discretization of the second derivative. We choose a typical spectral approximation using the discrete Fourier transform, which also imposes periodic boundary conditions. The CAP is similarly discretized as
where is a non-negative function which increases as one approaches the boundary of the interval and is zero in most of the interior.
Concerning the density operator , we write , with . Initially the system is prepared in a two-particle pure state localized inside the absorption-free part of the grid. The master equation for is equivalent to the usual Schrödinger equation (2) for , while the master equation for acquires a source term, i.e.,
It is natural to view discrete the one-particle wave function as a vector of length equal to the number of discrete degrees of freedom. Similarly, it is natural to view two-particle wave functions as anti-symmetric matrices . Moreover, a one-particle density matrix becomes a Hermitian matrix . Using these notions, the master equation (3) can be compactly written as
where the last term contains matrix-matrix products. The extra factor 2 in originates from the fact that the matrix relates to a (redundant) “basis” of direct product states.
If the initial two-particle state is an eigenstate of the total spin and its projection, the source term is always proportional to a single one-particle spin state. Hence, the spin does not introduce any complication in the notation in this case. For parallel spins, the one particle spin has the same direction as the two-particle spins, and for spin projection the one-particle density operator has its spinor component given by , where is the total spin eigen value.
The equation for becomes
In principle, it is not necessary to include this equation in our calculations, as can be calculated from the constraint .
3.1 Collision in a Gaussian well
We will now focus on an example in which we set and place the electrons in a potential of Gaussian form,
The particles interact via a regularized Coulomb interaction
For this problem we chose a CAP of power form:
where and is the distance from the edges at which the CAP is “turned on”, see figure 2.
The system may for instance serve as a model for a quantum dot which couples to the conduction band and has narrow confinement in two dimensions.
Equation (33) is integrated using a scheme of second order in the time step based on a standard split-step operator scheme . It is instructive to consider a more general setting in order to introduce the time stepping scheme for the density matrix. Given a differential equation for an entity in a linear space on the form
where is a linear operator dependent on and is a source term independent of , we may integrate formally using standard time-ordering techniques to obtain
with . The source terms are not easily transformed using the time-ordering operator . Assuming that the case can be integrated to -th order in the time using , a -th order method for the general case can be obtained by keeping source terms and evaluating these with sufficiently high order quadrature. For example, using standard Strang splitting which gives a scheme of local error , or other schemes based on the Magnus expansion , we may approximate the term as (using trapezoidal quadrature) and the term as . Specifically, in our implementation we have used a second order scheme given by
Here is the kinetic energy. As this scheme is trace preserving to second order only, we have also solved (34) in order to check that our numerical time step is small enough to preserve the total trace.
Figure 3 shows the evolution of the particle density for a system in which a fermion collides with another one. In this case the potential depth and the width , the interaction strength , and the “softening” . The CAP has the power , the strength and . The starting point is a spatially anti-symmetrized state (spin triplet) consisting of a particle trapped in the well and an incoming wave packet of Gaussian shape. The trapped part corresponds to a superposition of the ground and the first excited one-particle states in the well, and the incoming wave function has mean momentum . It is seen that as absorption, both due to transmission and back-scattering, takes place, the two-particle density vanishes and a one-particle density emerges. It is also seen that, apart from in the absorption region, the total particle density, obtained by adding the two and one-particle densities, compares well to the “true” particle density obtained from solving the Schrödinger equation without absorber on a larger grid.
Due to the collision, there is a finite probability that both particles are absorbed. This is clearly seen in figure 4, which shows how the total trace is distributed between the sub-spaces , and as a function of time. In this particular case we have and . Also shown are the expectation value of the particle number and the purity of the system. Purity is reduced in two ways. Firstly, it is reduced as the trace becomes distributed between the three sub-systems, and secondly because is not a pure state within . This is seen from the fact that conditional purity, , converges towards 0.6, i.e. less than unity (not shown explicitly in the figure).
3.2 Laser ionization of a one-dimensional helium atom
In this next example we expose our system to an electric pulse of type
With this time-dependent perturbation, the one-particle Hamiltonian (27) acquires a time-dependent term, which in the length gauge representation reads
We have here set the charge of the particles (the electrons) to be . The static potential felt by the particles is chosen to be a regularized Coulomb potential,
whereas the interaction between them is still described by (36). By choosing a.u. and a.u., the ground state energy and the first ionization threshold coincide with those of a true three-dimensional atom, i.e. the ground state energy is a.u. and the ground state energy of ”He” is a.u.. By ”a.u.” we mean atomic units, defined by choosing the Bohr radius, the electron mass, the elementary charge and as units for their respective quantities. The ground state of the system, which is a spin singlet state, is easily obtained by propagation in imaginary time.
Figure 5 shows the evolution of the system exposed to a pulse of maximum field strength a.u., central frequency a.u. and a duration corresponding to three optical cycles. This central frequency corresponds to a photon energy which energetically allows for one photon double ionization. Rather than an absorber of power-form (37), we have here used a Manolopoulos-type absorber , which has the advantages that it is transmission free and dependent on one parameter only. Along with a figure showing how the partial traces, given by (23), evolve in time, and another one showing the time-dependence of the electric field, given by (41), we have included snapshots of the absolute values of the two-particle wave function and the one-particle density matrix . It is clearly seen that as absorption takes place in the two-particle sub-system, a one-particle density matrix emerges. Note that the lobes following the axes of the -coordinate system do not correspond to absorption of the second electron but rather loss of coherence within the one-particle sub-system. However, from the lower panel of figure 5, which clearly demonstrates how single ionization may be distinguished from double ionization, we see that there is a finite probability of absorbing both particles. In this case, specifically, the probability of ionizing only one electron is , and the probability of double ionization is .
In conclusion, we have demonstrated how the Lindblad equation in Fock space can be used to generalize the notion of absorbing boundary conditions for -particle systems. Specifically, the remainder of the system is preserved as a particle is absorbed. With this formalism it may be possible to describe the dynamics of unbound systems which otherwise would require an unrealistically large numerical grid.
As a consequence of this being a Markovian process, some coherence is lost, and the state after absorption is in general given by a mixed state rather than a wave function. Within Lindblad theory, the identification between the Lindblad operators and the creation and annihilation operators comes quite natural for complex absorbing potentials. We believe that the method outlined in this paper may be generalized to other kinds of non-Hermitian “interactions” as well.
Since it is considerably more involved to solve a master equation rather than a Schrödinger equation, c.f. , future work will aim to find lower rank methods for solving (2). Also, more sophisticated spatial approximations like sparse grids  may be utilized to simplify the treatment of more particles.
-  Jonathan Parker, K T Taylor, Charles W Clark, and Sayoko Blodgett-Ford. Intense-field multiphoton ionization of a two-electron atom. J. Phys. B, 29:L33, 1996.
-  M Uiberacker, Th. Uphues, M. Schultze, A. J. Verhoef, V. Yakovlev, M. F. Kling, J. Rauschenberger, N. M. Kabachnik, H. Schröder, M. Lezius, K. L. Kompa, H.-G. Muller, M. J. J. Vrakking, S Hendel, U. Kleineberg, U. Heinzmann, M. Drescher, and F. Krausz. Attosecond real-time observation of electron tunnelling in atoms. Nature, 446:627, 2007.
-  R. Kosloff and D. Kosloff. Absorbing boundaries for wave-propagation problems. J. Comput. Phys., 63(2):363–376, 1986.
-  Thomas Fevens and Hong Jiang. Absorbing boundary conditions for the schrödinger equation. SIAM J. Sci. Comput., 21(1):255–282, 1999.
-  Xavier Antoine, Anton Arnold, Christophe Besse, Matthias Ehrhardt, and Achim Schaedle. A review of transparent and artificial boundary conditions techniques for linear and nonlinear schrodinger equations. Comm. Comput. Phys., 4(4):729–796, 2008.
-  JP Berenger. A perfectly matched layer for the absorption of electromagnetic-waves. J. Comput. Phys., 114(2):185–200, OCT 1994.
-  David E. Manolopoulos. Derivation and reflection properties of a transmission-free absorbing potential. J. Chem. Phys., 117(21):9552–9559, 2002.
-  M Monnerville and JM Robbe. Optical potential discrete variable representation method in the adiabatic representation - application to the co(b-1 sigma(+)-d ‘(1)sigma(+), j=0) predissociation process. Eur. Phys. J. D, 5(3):381–387, 1999.
-  TP Grozdanov, L Andric, and R McCarroll. Calculations of partial cross sections for photofragmentation processes using complex absorbing potentials. J. Chem. Phys., 124(9), 2006.
-  Hans O. Karlsson. Stability of the complex symmetric lanczos algorithm for computing photodissociation cross sections using smooth exterior scaling or absorbing potentials. J Phys. B, 42:125205, 2009.
-  Kenneth C. Kulander. Multiphoton ionization of hydrogen: A time-dependent theory. Phys. Rev. A, 35(1):445–447, 1987.
-  M Protopapas, CH Keitel, and PL Knight. Relativistic mass shift effects in adiabatic intense laser field stabilization of atoms. J. Phys. B, 29(16):L591–L598, 1996.
-  G. Lindblad. On the generators of quantum dynamical semigroups. Comm. Math. Phys., 48:119, 1976.
-  M. Schlosshauer. Decoherence and the quantum-to-classical transition. Springer-Verlag, 2007.
-  A Sandulescu and H Scutaru. Open quantum-systems and the damping of collective modes in deep inelastic-collisions. Ann. Phys., 173(2):277–317, 1987.
-  Paweł Caban, Jakub Rembieliński, Kordian A. Smoliński, and Zbigniew Walczak. Unstable particles as open quantum systems. Phys. Rev. A, 72(3):032106, 2005.
-  Reinhod A. Bertlmann, Walter Grimus, and Beatrix C. Hiesmayr. Unstable particles as open quantum systems. Phys. Rev. A, 73(5):054101, 2006.
-  M.D Feit, J.A Fleck Jr., and A Steiger. Solution of the schrödinger equation by a spectral method. J.Comput. Phys., 47(3):412 – 433, 1982.
-  C. Lubich. From Quantum to Classical Molecular Dynamics: Reduced Models and Numerical Analysis. European Mathematical Society, 2008.
-  Jean Dalibard, Yvan Castin, and Klaus Mølmer. Wave-function approach to dissipative processes in quantum optics. Phys. Rev. Lett., 68(5):580–583, 1992.