Absorbing boundary conditions for dynamical many-body quantum systems

1.4 MB·English

by  Sølve Selstø

Absorbing boundary conditions for dynamical
many-body quantum systems
0
1
0
SølveSelstø
2
CentreofMathematicsforApplications,UniversityofOslo,N-0316Oslo,Norway
n
a E-mail:[email protected]
J

8
2
SimenKvaal

] CentreofMathematicsforApplications,UniversityofOslo,N-0316Oslo,Norway
h
p E-mail:[email protected]
-
t
n
a Abstract. In numerical studies of the dynamics of unbound quantum mechanical
u
systems, absorbing boundary conditions are frequently applied. Although this
q
[ 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
2
v is due to the fact that all information about the system is lost upon absorption of
6 one particle; a formalism based solely on the Scrho¨dinger equation is not able to
8 describetheremainderofthesystemasparticlesarelost. Herewedemonstratehow
0
the dynamics of a quantum system with a given number of identical fermions may
2
. be described in a manner which allows for particle loss. A consistent formalism
4
whichincorporatestheevolutionofsub-systemswith areducednumberofparticles
0
9 isconstructedthroughtheLindbladequation. Specifically,thetransitionfroman N-
0 particle system to an (N 1)-particle system due to a complex absorbing potential
: −
v isachievedbyrelatingtheLindbladoperatorstoannihilationoperators. Themethod
i
X allows for a straight forward interpretation of how many constituent particles have
leftthe system after interaction. We illustrate the formalism using one-dimensional
r
a two-particlemodelproblems.
PACSnumbers:31.15.-p,03.65.Ca,32.80.-t
Absorbingboundaryconditionsfordynamicalmany-bodyquantumsystems 2
1. Introduction and basictheory
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 bephoto ionization ofhelium [1] and
cascaded Auger processes following ionization of inner core electrons in atoms [2].
Ofcourse,inordertounderstandsuchprocesses,theoreticalandnumericalstudiesare
required. Furthermore,thetheoreticalinterestofsuchsystems,beitwithinthecontext
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
thattheirspatialextensionmaybecomearbitrarilylarge. Moreover,evenifoneisable
to represent the whole system within a finite space, extracting relevant information
from thefinal wavefunctionmaybefar fromtrivial.
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
wavefunctionbelongingto theinteraction region. In a numericalimplementationthis
cannotbedonebysimplyresortingtoarepresentationofspacethatissmallerthanthe
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 [5]. Such absorbers are frequently referred to as
perfectlymatchedlayers inthecontextofgeneral waveequations[6].
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
theHamiltonianof thesystem[7]. Complex absorbingpotentialsare widely used e.g.
within molecular dynamics [8, 9, 10] and atomic physics [11, 12]. Alternatively, this
wayofabsorbingparticlesmaybeformulatedasmultiplyingthewavefunctionwitha
maskingfunctionat each timestep[1].
In anycase theresultingeffectiveHamiltonianacquires an anti-Hermitianpart,
H H iG , (1)
eff
≡ −
whereboth H and G areHermitianand G ispositivesemi-definite(G 0).
≥
Thewavefunction Y (t) ofthesystemobeystheSchro¨dingerequation
| i
d
ih¯ Y (t) =H Y (t) . (2)
eff
dt| i | i
A densityoperatorρ(t)correspondinglyobeysthevonNeumannequation
d
ih¯ ρ(t)=H ρ(t) ρ(t)H† =[H,ρ(t)] i G ,ρ(t) . (3)
dt eff − eff − { }
Absorbingboundaryconditionsfordynamicalmany-bodyquantumsystems 3
It iseasy tosee thattheevolutionis non-unitary,viz,
d 2
Tr[ρ(t)]= Tr[G ρ(t)] 0, (4)
dt −h¯ ≤
so that probability is “lost” if the density matrix overlaps with anti-Hermitian part of
theeffectiveHamiltonian.
Althoughmethodsinvolvingnon-HermitianeffectiveHamiltoniansmay comein
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
functionislost. Thisisobvioussincethewavefunctionisnormalizedtotheprobability
of finding all particles within the space defined by the numerical implementation. In
otherwords: IfaninitialN-particlesystemgradually“loses”oneparticlethroughsome
non-Hermitian“interaction”,thecorrespondingwavefunction Y (t) graduallygoes
N
| i
tozero –not tosomewave-functioncorrespondingto N 1 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 (N 1)-particle wave
−
function Y (t) iscreatedasoneparticleislost,andwherethiswavefunctionmay
N 1
| − i
be propagated using the corresponding Schro¨dinger equation including some source
term. Once being able to do this, the extension to N 2 particles etc. follows by
−
induction. However, as it turns out, such a construction is not possibledue to the fact
thattheprocessoflosingparticlesinthiswayisirreversible: Asaparticleisabsorbed,
informationisirretrievablylost. Hence,aMarkovianmasterequationshouldbeamore
suitablestartingpointthan a purestateapproach. Lindblad[13]was able toshow that
in order for the evolutionof a system to be Markovian, trace-conserving and positive,
thedensityoperatorρ(t)hasto obeyan equationoftheform
d
ih¯ ρ=[H,ρ] D(ρ) (5)
dt −
withtheso-calledLindbladianD(ρ)givenby thegenericexpression[14]
D(ρ)=i(cid:229) γ A†A ρ+ρA†A 2A ρA† . (6)
m,n m n m n n m
−
m,n (cid:16) (cid:17)
Here, H is the Hamiltoniangoverningthe unitary part of the evolution. The operators
A arereferred toasLindbladoperators. Inadiagonalrepresentation,theLindbladian
n
simplifiesto
D(ρ)=i(cid:229) A†A ρ+ρA†A 2A ρA† . (7)
n n n n n n
−
n (cid:16) (cid:17)
Equation (5) is usually used to describe energy dissipation to the environment [15].
However, it has recently been demonstrated that the spontaneous decay of unstable
particles may also be described within this formalism [16, 17]. In reference [16] a
masterequationofLindbladfromisobtainedforaHilbertspaceconsistingofunstable
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[17]decayfromonesystemtoanotheroneconsistingofitsdecayproductsis
Absorbingboundaryconditionsfordynamicalmany-bodyquantumsystems 4
describedbyintroducingasingleLindbladoperator. Intheseworksitisdemonstrated
thatonemay describedecoherence alongwithdecay inthiscontext.
In the present work these ideas have been used to generalize the standard one-
particle absorbing boundary technique to an N-body setting, and it is demonstrated
that this is indeed the naturalway to do this. Alternatively,in a moregeneral context,
this formalism allows the study of a system which gradually loses particles to some
environmentdueto thenon-Hermitianpart of H .
eff
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,inordertoobtainphysicallymeaningfulresults,anydependenceof
theabsorbermustbeeliminatedbyplacingtheabsorbersufficientlyfarawayfromthe
interactionregion–asisthecaseforanyimplementationfeaturingabsorbingboundary
conditions.
In thefollowingsection,section2,theformalismwillbeformulatedforidentical
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,
conclusionsare drawnand afew futureperspectivesare outlined.
2. Fockspace description
The natural setting for describing a system of a variable number of fermions is Fock
space, the direct sum of all n-fermion spaces H . As we wish to describe at most N
n
particles,it sufficesto consider
N
H = H . (8)
n
Mn=0
An arbitrary n-fermionstate Y H can bewritten
n
| i∈
1
Y = dxnY (x x )ψ†(x ) ψ†(x ) , (9)
1 n 1 n
| i n!Z ··· ··· |−i
wherethefieldoperatorψ†(x)createsaparticleatposition x. Thefieldoperatorsobey
theusualanti-commutatorrelations
ψ(x),ψ(x) =0, ψ(x),ψ†(x) =δ(x x ). (10)
′ ′ ′
{ } { } −
Here, x may denote all degrees of freedom associated with a particle, including
eigenspin. Moreover, Y (x x ) is the anti-symmetric local wave function of the
1 n
···
n-particlesystem,and H is thevacuumstate, containingzero particles.
0
|−i∈
Asystemofnidenticalfermionsinteractingwithatmosttwo-bodyforceshasthe
Hamiltonian
n n
(cid:229) (cid:229)
H = h(x)+ u(x,x ). (11)
i i j
i=1 j<i
Here, h(x) (resp. u(x,x )) is a one-body (two-body) operator acting on the degrees-
i i j
of-freedom associated with particle i (and j). Using field operators, the Hamiltonian
Absorbingboundaryconditionsfordynamicalmany-bodyquantumsystems 5
can bewritten
H = dxψ†(x)h(x)ψ(x)+ (12)
Z
1
dxdx ψ†(x)ψ†(x)u(x,x)ψ(x)ψ(x),
′ ′ ′ ′
2Z
thepointherebeingthatH givenonthisformisindependentofthenumberofparticles
in the system. The exact expression is irrelevant at this point. The CAP, which is
diagonalin x, isconvenientlyintroducedas
iG = i dxψ†(x)G (x)ψ(x). (13)
− − Z
By comparing (5), (7) and (13) with the von Neumann equation (3), we see that
theLindbladoperators, foradiagonalrepresentationoftheLindbladian,mustfulfil
G = dxA†(x)A(x) (14)
Z
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, 2A ρA†
n n
− →
2A(x)ρA†(x), should map an N-particle system into an (N 1)-particle system.
− −
Hence, A(x) should map H into H , which means that the Lindblad operator
N N 1
−
A(x) mustbeontheform
A(x)= dxa(x,x )ψ(x). (15)
′ ′
Z
Thesimplestchoicethatsatisfies (14)is thediagonalform
A(x) G (x)ψ(x). (16)
≡
p
We will return to the justificationof why thisis theproper way to define theLindblad
operators shortly.
With(16), theLindbladian(7), may immediatelybewritten
D(ρ)=iG ρ+iρG 2i dxG (x)ψ(x)ρψ†(x), (17)
− Z
and themasterequation(5)becomes
d
ih¯ ρ=[H,ρ] i G ,ρ +2i dxG (x)ψ(x)ρψ†(x). (18)
dt − { } Z
ThisisourfundamentaldynamicalformulationofparticlelossduetoaCAP.
Let us consider the master equation (18) in some detail. We may partition the
densitymatrixρintoblocks,viz,
N N
(cid:229) (cid:229)
ρ= ρ , (19)
n,m
n=0m=0
where ρ = P ρP , with P being the orthogonal projector onto H . Each block
n,m n m n n
evolvesaccording to
d
ih¯ ρ =[H,ρ ] i G ,ρ + (20)
m,n m,n m,n
dt − { }
2i dxG (x)ψ(x)ρ ψ†(x),
m+1,n+1
Z
Absorbingboundaryconditionsfordynamicalmany-bodyquantumsystems 6
showingthattheflowofρ dependsonthatofρ ,butnottheotherwayaround.
n,m n+1,m+1
This is illustrated in figure 1. Also, notice that ρ obeys the original von Neumann
N,N
equation (3) as there are no couplings to other blocks in this case as particles do not
enterthe N-particlesystem.
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 N-particle system to an (N 1)-particle system, only
−
the unabsorbed part of the original system should be reproduced within the (N 1)-
−
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, Y =
[α(x ,t)β(x ,t) β(x ,t)α(x ,t)]/√2. The state α does not overlap with the absorber
1 2 1 2
−
at anytime,i.e.
G (x)α(x,t)=0, forall xand t, (21)
and β corresponds to an unbound particle travelling outwards. We suppose that as
t approaches infinity, the particle in the state β is completely absorbed, and our
numerical representation of the final system should converge towards the pure one-
particlestateα.
Indeed, the evolution dictated by (20) follows this pattern. The block of the
density matrix corresponding to the two-particle system, ρ , remains a pure state
2,2
– albeit with decreasing norm as β is absorbed – and the evolution of the one-particle
blockρ simplifiesdueto(21)to
1,1
ih¯ρ˙ =[h,ρ ]+2i G (x) β(x;t) 2dx α α . (22)
1,1 1,1
Z | | | ih |
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
ρ (t ¥ )= α α ,as itshould.
1,1
→ | ih |
On the otherhand, with anon-diagonal form ofthe Lindbladoperators, c.f. (15),
we would have contributions to ρ˙ of the form β α and its Hermitian adjoint in
1,1
| ih |
addition to the pure state contribution. These contributionsclearly cannot be allowed,
as ρ should be independent of the state β of the outgoing particle. Hence, we find
1,1
that,uptoanarbitraryphasefactor,thedefinition(16)istheproperwaytodefineA(x).
For a non-interacting N-fermion system where N particles are bound and N
α β
particles are unbound, and with an initial state given by the Slater determinant Y =
| i
α ,...,α ,β ,...,β , it may be verified by inspection that the source term for the
|{ 1 Nα 1 Nβ}i
N -th block, ρ , is always proportional to α ,...,α α ,...,α given that
α Nα,Nα |{ 1 Nα}ih{ 1 Nα}|
noneoftheα-statesoverlapwithG . Theintermediatedensitymatrixblocksρ ,
i N n,N n
where0<n<N ,willapproach zero ast ¥ , buttransientlydescribe(mixed−)sta−tes
β
→
wheren particleshaveleft.
Ofcourse,inthemoreinterestingcontextofinteractingparticles,thestructureof
the diagonal blocks of the density operator, ρ , is more complex than in the special
n,n
caseofnon-interactingparticles.
Absorbingboundaryconditionsfordynamicalmany-bodyquantumsystems 7
2.1. Consequences and interpretation
Typically,ourstartingpointisapureN-particlestate, Y (0) Y (0) ,inwhichcase(20)
| ih |
reduces totheordinarySchro¨dingerequation(2),which wasouroriginalformulation.
Moreover, it is easy to see that ρ(t) will remain block diagonal for all t in this case,
i.e. ρ = 0 if n = m. But (20) shows that ρ (t) in general is a mixed state, unlike
n,m n,n
6
ρ (t)= Y (t) Y (t) . This is due to the information loss when admittingignorance
N,N
| ih |
ofthewhereaboutsoftheremovedparticle.
We stress that by construction, Tr[ρ(t)] = 1 for all t. Probability flows
monotonically from H into H , and the particle absorbed from ρ is not present
N n N,N
in ρ , but is erased, and ρ is a proper description of an n-fermion system. In this
n,n n,n
way,weseethattheaboveconstructionisanaturalgeneralizationoftheoriginalnon-
Hermitian dynamics, which describes the classical removal of a particle, into one that
alsodescribestheremainingsysteminaconsistentway. Itisstrikingtonoticethatthe
CAP, iG (x), isalreadygivenas oneofthetermsintheLindbladian,so thattheFock
−
spaceformulationfollowsalmostimmediately.
Itisworthwhiletomentionthatthetracesoftheblocksρ alongthediagonalof
n,n
ρ have very simple interpretations as the probability P(n) of having n particles in the
systemuponameasurement,i.e.
P(n) Tr [ρ(t)] Tr[ρ (t)]. (23)
n n,n
≡ ≡
In particular, P(N) = Y Y 1, and P(0) = 1 (cid:229) N P(n). Hence, within
h | i ≤ − n=1
this formalism, distinguishing between single, double etc. ionization of atoms and
moleculesisstraightforward.
ForanyobservableAtheexpectationvalueisgivenby A Tr[Aρ]=(cid:229) Tr[Aρ ].
n n
h i≡
Forexample,theexpected numberofparticlesin thesystemisgivenby
N
N = dxTr[ψ†(x)ψ(x)ρ]= (cid:229) nP(n). (24)
h i Z
n=1
It mayalso beusefultoconsiderconditionalexpectationvalues:
Tr (Aρ)
n
A , (25)
n
h i ≡ Tr (ρ)
n
i.e., theexpectationvalueof A given thatthesystemisfoundin an n-particlestate.
Obviously, as particles are “removed” by the absorber, information is lost. This
informationlossmaybequantifiedby thevon Neumannentropy,S Tr[ρlogρ],or
≡−
thecloselyrelated notionofpurity,defined by [14]
ς Tr(ρ2) 1. (26)
≡ ≤
Unity minus this quantity, 1 ς, 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
quantityofthere-normalized block,i.e. ρ /Tr[ρ ].
n,n n,n
Of course, for a full quantum mechanical description in terms of the complete
(unabsorbed) N-particle wave function, no information is lost. Indeed, the absorption
Absorbingboundaryconditionsfordynamicalmany-bodyquantumsystems 8
of particles is a semi-classical concept. In the full N-fermion quantum system no
such separation is possible. On the other hand, as the Schro¨dinger equation for the
N particles is separable in thenon-interacting α and β states discussed above, giving
i i
a Slater determinant of the time evolved one-particle states as the full solution, the
Lindblad equation is seen to correctly construct the N -particle Slater determinant
α
resulting from the removal of the N outgoing particles from this Slater determinant.
β
Thus, the Lindblad equation exactly encapsulates the approximate separation of non-
interactingquantumsystemsfar apart.
3. Example: Two identical spin 1-particles inonedimension
2
In the following we consider two fermions in one dimension with the one-body
Hamiltonianh givenby
h¯2 ∂2
h= +V(x,t), (27)
−2m∂x2
where V(x,t) is some external, possibly time-dependent, potential. The fermions
interactviaapotentialU(x x )andtheextensionofthesystemiseffectivelyreduced
1 2
−
by imposing a CAP. With two interacting identical fermions, we are confined to the
Hilbertspace
H =H H H . (28)
2 1 0
⊕ ⊕
Wediscretizethissystemusingauniformgridintheinterval[0,x )containing
max
N points x N 1, x = jh,withh=x /N. Thefieldoperatorsψ†(x)canbereplaced
{ j}j=−0 j max
†
by a finite number of creation operators c , creating a particle at grid point x . These
j j
operators, whichobeytheusualfermionanti-commutatorrules
c ,c =0, c ,c† =δ , (29)
{ j k} { j k} j,k
mapdiscreteanti-symmetrizedδ-function bases forH intobases forH .
n n 1
±
TheHamiltoniantakes theform
H =(cid:229) h c†c +1(cid:229) U(x x )c†c†c c, (30)
ij i j 2 i− j i j j i
i,j i,j
where h are the matrix elements of the one-body Hamiltonian, which depend on
ij
the chosen discretization of the second derivative. We choose a typical spectral
approximation using the discrete Fourier transform, which also imposes periodic
boundaryconditions. TheCAP issimilarlydiscretized as
iG = i(cid:229) G (x )c†c , (31)
− − j j j
j
whereG (x)isanon-negativefunctionwhichincreasesasoneapproachestheboundary
oftheinterval[0,x )and iszero inmostoftheinterior.
max
Concerning the density operator ρ, we write ρ = Y Y +ρ +ρ , with ρ =
1 0 n
| ih |
P ρP . Initially the system is prepared in a two-particle pure state localized inside
n n
the absorption-free part of the grid. The master equation for ρ (t) is equivalent to the
2
Absorbingboundaryconditionsfordynamicalmany-bodyquantumsystems 9
usualSchro¨dingerequation(2)for Y (t) ,whilethemasterequationforρ (t)acquires
1
| i
asourceterm,i.e.,
d
ih¯ ρ (t)=[H,ρ (t)] G ,ρ (t) + (32)
1 1 1
dt −{ }
2i(cid:229) G (x )c Y (t) Y (t) c†.
j j| ih | j
j
It is natural to view discrete the one-particle wave function as a vector ψ (x ) of
1 j
length equal to the number of discrete degrees of freedom. Similarly, it is natural to
view two-particle wave functions as anti-symmetric matrices ψ (x ,x ). Moreover, a
2 j k
one-particledensitymatrixbecomesaHermitianmatrixρ (x,x). Usingthesenotions,
1 ′
themasterequation(32)can becompactlywrittenas
d
ih¯ ρ (t) =[H,ρ (t)] G ,ρ (t) +2iρ (33)
1 1 1 S
dt −{ }
with ρ 2hψ†G ψ ,
S ≡ 2 2
wherethelasttermcontainsmatrix-matrixproducts. Theextrafactor2inρ originates
S
fromthefactthatthematrixψ relatestoa(redundant)“basis”ofdirectproductstates.
2
If the initialtwo-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 M = 0 the one-particle density operator has its spinor component
S
givenby ( +( 1)S )/√2, whereS isthetotalspineigenvalue.
|↑i − |↓i
Theequationforρ = p (t) becomes
0 0
|−ih−|
d p (t)= 2(cid:229) G (x )c ρ c†. (34)
dt 0 h¯ j j 1 j
j
In principle, it is not necessary to include this equation in our calculations, as p (t)
0
can becalculated from theconstraintTr[ρ(t)]=1.
3.1. Collisionina Gaussianwell
We will now focus on an example in which we set h¯ =m=1 and place the electrons
ina potentialofGaussianform,
(x x )2
0
V(x)= V exp − . (35)
− 0 (cid:18)− 2σ2 (cid:19)
Theparticlesinteract viaaregularizedCoulombinteraction
λ
U( x x )= . (36)
1 2
| − | (x x )2+δ2
q 1− 2 C
Forthisproblemwechosea CAPofpowerform:
n
ξ
G (x)=C , (37)
(cid:18)x (cid:19)
T
ξ max 0,x x,x (x x ) ,
T max T
≡ { − − − }
Absorbingboundaryconditionsfordynamicalmany-bodyquantumsystems 10
where n 1 and x is the distance from the edges at which the CAP is “turned on”,
T
≥
seefigure2.
The system may for instance serve as a model for a quantum dot which couples
totheconductionbandand has narrowconfinement intwo dimensions.
Equation(33)isintegratedusingaschemeofsecondorderinthetimestepbased
on a standard split-step operator scheme [18]. It is instructive to consider a more
general setting in order to introduce the time stepping scheme for the density matrix.
Givena differentialequationforan entityy(t)in alinearspaceontheform
y˙(t)=L y(t)+ f(t), (38)
t
where L is a linear operator dependent on t and f(t) is a source term independent of
t
y(t),wemayintegrateformallyusingstandard time-orderingtechniquesto obtain
y(t)=T e 0tLsdsy(0)+F(t)+
R
¥
t s
(cid:229) n 1
− L L F(s )ds ds , (39)
Z0 ···Z0 s1··· sn n n··· 1
n=1
n-fold
t | {z }
with F(t) = f(s)ds. The source terms are not easily transformed using the time-
0
orderingoperRatorT . Assumingthatthecase f(t)=0canbeintegratedto p-thorderin
thetimetusingy(t)=U y(0),a p-thordermethodforthegeneralcasecanbeobtained
t
bykeeping psourcetermsandevaluatingthesewithsufficientlyhighorderquadrature.
Forexample,usingstandardStrangsplittingwhichgivesaschemeoflocalerrorO(t3),
or other schemes based on the Magnus expansion [19], we may approximatethe term
F(t) as F(t) t[f(0)+ f(t)]/2 (using trapezoidal quadrature) and the n = 1 term as
≈
tL F(s)ds t2L [f(0)+tf˙(0)/2]. Specifically, in our implementation we have
0 s t/2
≈
Rused asecondorderschemegivenby
ρ (t)=e i(V+U iG )t/2e iTte i(V+U iG )t/2ρ (0)
1 − − − − − 1
×
e+i(V+U+iG )t/2e iTte+i(V+U+iG )t/2+2tρ (0)+
− S
t2 ρ˙ (0) i(Hρ (0) ρ (0)H†) +O(t3) (40)
S S S
− −
h i
where ρ˙ =2h ψ˙ G ψ†+ψ G ψ˙† .
S 2 2 2 2
(cid:16) (cid:17)
Here T is the kinetic energy. As this scheme is trace preserving to second order only,
wehavealsosolved(34)inordertocheckthatournumericaltimestepissmallenough
topreservethetotaltrace.
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 V = 4 and the
0
width σ = 0.75, the interaction strength λ = 5, and the “softening” δ = 0.1. The
C
CAP has the power n = 3, the strength C = 4 and x = 5. The starting point is a
T
spatially anti-symmetrized state (spin triplet) consisting of a particle trapped in the
wellandanincomingwavepacketofGaussianshape. Thetrappedpartcorrespondsto
asuperpositionofthegroundandthefirstexcitedone-particlestatesinthewell,andthe
incomingwavefunctionhasmeanmomentumk =2. Itisseenthatasabsorption,both
0
duetotransmissionand back-scattering, takesplace, thetwo-particledensityvanishes