Science Question

Hello 👋

How are you?

I have the task of writing an introduction to a scientific research so that the sentences must be rephrased in a scientific way that is far from cheating. Can you help me with this and write the references used in Leeds Harford format. I will send you a sample method for writing an introduction, which is a comparison of some of the methods used to solve some physical problems. The required references are three scientific papers in addition to the references used in the same papers.

I will send the 3 papers

But do not depend only on this three But also all references inside this three.

There are a number of styles in writing references, including a method called Leeds Harford style, and this is what I need.

I need to have coordination of ideas, coherence, and include the use of linking phrases between sentences.

I want to write an introduction, but I wanted to draw attention to the importance of clarifying the wording of the sentences, and not just copying them from the main reference

The idea of the introduction I will review several methods that have been used and compare them, then I will mention the method that I will use in my research

Please, I want to write an introduction in a professional manner, because it is a scientific paper for publication. Also, I want the citation at least three references.

Do not limit yourself to the number of words, the important thing is that the idea is complete.

Please, It important to gives a comprehensive idea of all the methods mentioned in the papers sent to you

As for the references, they must be written in the approved style that I sent to you.

Thank you

Multiconfigurational quantum propagation
with trajectory-guided generalized
coherent states
Cite as: J. Chem. Phys. 144, 094106 (2016); https://doi.org/10.1063/1.4942926
Submitted: 11 November 2015 . Accepted: 16 February 2016 . Published Online: 04 March 2016
Adriano Grigolo, Thiago F. Viscondi, and Marcus A. M. de Aguiar
ARTICLES YOU MAY BE INTERESTED IN
Quantum mechanics with the basis set guided by Ehrenfest trajectories: Theory and
application to spin-boson model
The Journal of Chemical Physics 130, 244101 (2009); https://doi.org/10.1063/1.3153302
Basis set sampling in the method of coupled coherent states: Coherent state swarms,
trains, and pancakes
The Journal of Chemical Physics 128, 054102 (2008); https://doi.org/10.1063/1.2828509
Zombie states for description of structure and dynamics of multi-electron systems
The Journal of Chemical Physics 148, 194109 (2018); https://doi.org/10.1063/1.5023209
J. Chem. Phys. 144, 094106 (2016); https://doi.org/10.1063/1.4942926
© 2016 AIP Publishing LLC.
144, 094106
THE JOURNAL OF CHEMICAL PHYSICS 144, 094106 (2016)
Multiconfigurational quantum propagation with trajectory-guided
generalized coherent states
Adriano Grigolo,1,a) Thiago F. Viscondi,2,b) and Marcus A. M. de Aguiar1,c)
1
Instituto de Física Gleb Wataghin, Universidade Estadual de Campinas, 777 Sérgio Buarque de Holanda,
13083-859 Campinas, Brazil
2
Instituto de Física, Universidade de São Paulo, 1371 Rua do Matão, 05508-090 São Paulo, Brazil
(Received 11 November 2015; accepted 16 February 2016; published online 4 March 2016)
A generalized version of the coupled coherent states method for coherent states of arbitrary Lie
groups is developed. In contrast to the original formulation, which is restricted to frozen-Gaussian
basis sets, the extended method is suitable for propagating quantum states of systems featuring
diversified physical properties, such as spin degrees of freedom or particle indistinguishability. The
approach is illustrated with simple models for interacting bosons trapped in double- and triplewell potentials, most adequately described in terms of SU(2) and SU(3) bosonic coherent states,
respectively. C 2016 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4942926]
I. INTRODUCTION
A vast number of physical systems exhibit the property
that some of their parts behave in a sort of classical way,
meaning that quantum effects play only a minor role in
the description of those parts. This distinctive classical
character of specific degrees of freedom is a much welcomed
attribute, for it makes possible the development of tractable
computational approaches capable of carrying out the timeevolution of complex quantum systems, being thus the
fundamental property upon which time-dependent trajectoryguided methods are based.
In this kind of technique, quantum states are represented
in terms of time-dependent basis functions or “configurations.”
Within a single configuration, those degrees of freedom
in which quantum effects are negligible are evolved
according to classical equations of motion. This classical
dynamics may be prescribed in a number of different ways
and different choices correspond to different propagation
schemes.
In spite of the fact that individual configurations have
some of their parts bound to obey classical laws, a
complete quantum solution is in principle attainable by
combining many configurations. The key idea behind such
“multiconfigurational” approaches is that trajectory-guided
basis functions are more likely to remain in the important
regions of the Hilbert space, thus being more efficient at
representing the quantum state in the sense that a reduced
number of basis elements are required in order to achieve an
accurate description. And it is precisely through a significant
reduction in the number of basis functions needed to
propagate the system that one hopes to escape the exponential
scaling of basis-set size with dimensionality typical of
standard static-basis formulations. This “mixed quantuma)Electronic mail: agrigolo@ifi.unicamp.br
b)Electronic mail: viscondi@if.usp.br
c)Electronic mail: aguiar@ifi.unicamp.br
0021-9606/2016/144(9)/094106/15/$30.00
classical” picture is adopted in many methods of quantum
chemistry.1
A recurrent theme in this field is the development of
techniques which, by means of equally simple recipes to
guide the basis functions, would be readily applicable to
systems presenting authentically non-classical qualities, such
as spin degrees of freedom or particle exchange symmetry.
Several works have been directed to that purpose, most often
aiming at a time-dependent description of the electronic
structure of molecules during non-adiabatic processes. One
particular example of such a recipe is the classical model
for electronic degrees of freedom proposed by Miller and
White2 in which a second-quantized fermionic Hamiltonian is
properly reduced to a classical function wherein number and
phase variables play the role of generalized coordinates. In
contrast, a more “mechanistic” approach to fermion dynamics
is found on the multiconfigurational formula proposed by
Kirrander and Shalashilin3 in which the basis functions consist
of antisymmetrized frozen Gaussians4 guided by fermionic
molecular dynamics.5
Yet if one seeks to describe non-classical degrees of
freedom by means of classical-like variables, then generalized
coherent states — defined in the group-theoretical sense —
are indisputably the most appropriate tools to be employed.
There are many reasons supporting this assertion.
First of all, coherent states are defined in terms of
non-redundant parameters and equations of motion for these
parameters can be readily obtained from the time-dependent
variational principle (TDVP).6 In this way, an optimized
time evolution can be assigned in an unambiguous manner.
Moreover, they are naturally able to capture the desired
symmetries of the system which are maintained during
propagation. Furthermore, the coherent-state parameters
evolve in a classical phase space in the strict sense of
the word, hence we automatically have at our disposal the
wealth of analytical techniques applicable to Hamiltonian
systems. At the same time, through this intimate connection
to classical dynamics, coherent states provide a compelling
144, 094106-1
© 2016 AIP Publishing LLC
094106-2
Grigolo, Viscondi, and de Aguiar
classical interpretation to quantum phenomena, in so far as
individual configurations are chosen to represent familiar
objects — i.e., in such a way that it is meaningful to discuss
the dynamics of the system in terms of their trajectories.
To this extent, coherent states — which are also minimum
uncertainty states (provided a proper meaning is assigned to
the term “uncertainty”)7,8 — are valuable tools in enhancing
our comprehension with respect to the semiclassical features
of the quantum system under investigation. In addition, the
group-theoretical formalism secures a well-defined integral
form for the coherent-state closure relation9 — a crucial
element to the developments presented in this work. This
list of virtues is not exhausted and other advantages of a
generalized coherent-state representation will be pointed out
throughout the paper.
Along these lines, Van Voorhis and Reichman10 have
considered a number of alternative representations of
electronic structure making use of different coherent-state
parametrizations and also examined their adequacy to a
variety of systems.12 Within the context of non-adiabatic
molecular dynamics, a particularly interesting coherent-state
representation13,14 is found in the simplest and most throughly
investigated version of the electron-nuclear dynamics theory,
developed by Deumens, Öhrn, and collaborators.15,16 The
same kind of coherent state has been discussed in detail, within
the field of nuclear physics, by Suzuki and Kuratsuji.17–19
Turning to bosonic dynamics, a semiclassical trajectory-based
formula in the SU(n) coherent-state representation has been
recently derived and successfully applied to a model of trapped
bosons.20,21
These methods are representative of the kind of
technique one has in mind when a description of
intrinsically quantum degrees of freedom in terms of
classical-like variables is desired. However, they either
constitute approximate single-configuration approaches15,16
or involve complicated trajectories that live in a duplicated
phase space,20,21 sometimes relying on sophisticated rootsearch techniques in order to determine them.10,22,23 It
seems that a multiconfigurational, generalized coherent-state
approach, based on simple — as opposed to duplicated —
phase-space trajectories would be more in the spirit
of the familiar time-dependent guided-basis methods of
quantum chemistry.24 This is precisely the direction we take
here.
In this work, a quantum initial value representation
method which employs a generalized coherent-state basis set
guided by classical trajectories is formulated. The resulting
propagation scheme is regarded as a natural extension of the
coupled coherent states (CCS) technique of Shalashilin and
Child25–27 in so far as (i) basis-set elements represent localized
quantum states; (ii) each element evolves independently in a
generalized classical phase space and carries an action phase;
and (iii) the quantum amplitudes associated with individual
elements obey fully coupled equations of motion which
present a number of attractive qualities.
We begin in Section II with an overview on generalized
coherent states, deliberately avoiding the underlying grouptheoretical formalism associated with their construction. In
particular, we demonstrate how the time-dependent variational
J. Chem. Phys. 144, 094106 (2016)
principle leads to classical equations of motion for the
coherent-state parameters in a curved phase space. This
preliminary discussion is illustrated with specific examples.
Next, in Section III, we set forth to derive the working
equations of the generalized coupled coherent states method,
first in continuum form and then in terms of a discrete basis
set. In the latter case, both unitary and non-unitary versions of
the formulas are devised. A primitive method, resulting from
propagation of a single configuration, is also discussed. In
Section IV, we use our approach to study model Hamiltonians
and compare our results against exact quantum data. We also
take the opportunity to expose certain particular aspects of the
generalized coherent-state formulation. Finally, conclusions
are drawn in Section V.
II. GENERALIZED COHERENT-STATE FORMALISM
Coherent states are most elegantly discussed within the
context of group theory and this is the point of view adopted
here. We shall not venture into the group-theoretical formalism
itself though — on that subject, the reader is referred to the
instructive review by Zhang, Feng, and Gilmore9 or to a recent
paper by the present authors.11 Here, we will rather follow
a more pragmatic approach according to which a coherent
state is given in the form of an expansion in a proper set
of orthonormal states and work its geometrical properties
thereon.
A. Preliminaries
Coherent states are Hilbert space vectors labeled by a
complex vector z = (z1, . . . , z d ).28 They can be understood as
the result of a Lie-group operator acting on a reference state,
which is recovered by setting all entries of the vector z to zero.
We shall denote a non-normalized coherent state by |z}. These
curly ket states are analytical in z, while the bra states {z| are
analytical in the complex conjugate variable, denoted by z ∗;
the normalized state |z⟩ is not analytical in z for it depends
1
on z ∗ through the normalization factor {z|z}− 2 . The length d
of z will be identified as the number of degrees of freedom of
the classical phase space associated with the dynamics of z.
Coherent states of different groups are characterized
by their distinct geometrical properties which, in turn, are
described in terms of a function f related to the scalar product
between two non-normalized coherent states,
f (z ∗, z ′) = log{z|z ′}.
(2.1)
The classical phase-space metric gα β is a Hermitian matrix
obtained by taking the cross derivatives of the real function
f (z ∗, z) with respect to its complex arguments, treating z and
z ∗ as independent variables,29
gα β (z ∗, z) =
∂ 2 f (z ∗, z)
.
∂zα ∂z ∗β
(2.2)
The (non-orthogonal) coherent states span an overcomplete
basis of the corresponding Hilbert space and a closure relation
094106-3
Grigolo, Viscondi, and de Aguiar
holds
J. Chem. Phys. 144, 094106 (2016)

1̂ =
dµ(z ∗, z)|z⟩⟨z|,
(2.3)
where the integration domain depends on the specific type of
coherent state being considered — for semisimple compact
Lie Groups or the Heisenberg-Weyl group, for example, it
extends over the entire d-dimensional complex plane. In
(2.3), the integration measure dµ(z ∗, z) is defined as
dµ(z ∗, z) = κ det[g(z ∗, z)]
d

d2 zα
π
α=1
,
(2.4)
where d2 zα = d(Rezα )d(Imzα ) and the constant κ is determined by normalization of (2.3) — e.g., by setting the
expectation value of (2.3) in the reference state to unity — and
therefore it depends on the quantum numbers that characterize
the particular Hilbert space that carries the coherent-state
representation (some examples are found in Section II C). In
order to shorten the notation, we shall write simply dµ(z), but
keeping in mind that the measure is a real function of both z
and z ∗.
B. Classical dynamics
The norm-invariant Lagrangian6 that gives rise to the
Schrödinger equation is


⟨ψ| Ĥ |ψ⟩
i~ ⟨ψ|ψ̇⟩ − ⟨ψ̇|ψ⟩
−
,
(2.5)
L(ψ) =
2
⟨ψ|ψ⟩
⟨ψ|ψ⟩
where |ψ⟩ = |ψ(t)⟩ is an arbitrary quantum state and Ĥ
represents the system’s Hamiltonian operator. The TDVP
states that, given a parametrization of |ψ⟩ in terms of some set
of variables, the Euler-Lagrange equations obtained from (2.5)
translate into equations of motion for those same variables.
If the parametrization is flexible enough, an exact quantum
solution is achieved. On the other hand, if the parametrization
contains less than the number of variables needed to span the
associated Hilbert space, the resulting dynamics will be only
approximate.
We shall define the classical equations of motion for z as
those equations obtained from the TDVP when a trial state
|ψ⟩ = |z} is substituted in (2.5) — that is, a trial state whose
dynamics is restricted to the nonlinear subspace consisting
only of coherent states. In this case, Equation (2.5) takes the
form

d 
i~  ∂ f (z ∗, z)
∂ f (z ∗, z) ∗
L(z) =
żα −
żα − H(z ∗, z), (2.6)
2 α=1
∂zα
∂zα∗
where the classical Hamiltonian is the diagonal element of
the operator Ĥ in the coherent-state representation,
{z| Ĥ |z}
= ⟨z| Ĥ |z⟩.
(2.7)
{z|z}
By means of the Euler-Lagrange equations, one
immediately finds that the dynamics of z obeys
H(z ∗, z) =
d

β=1
ż β gβα (z ∗, z) = −
i ∂H(z , z)
,
~ ∂zα∗
∗
(2.8)
and that an equivalent (complex conjugate) equation holds
for z ∗. Notice how the coherent-state geometry introduces a
curvature in phase space by means of the metric g. One now
can distinguish between two kinds of coupling between the
components of the vector z: a dynamical coupling via H and
a geometrical coupling induced by g.
The group-theoretical formalism assures us that Eq. (2.8)
describes a Hamiltonian system in the most strict sense:
the phase space exhibits a symplectic structure, since a
nondegenerate Poisson bracket can always be established.9
Furthermore, the measure (2.4) is invariant under the
classical dynamics given by Equation (2.8), that is, dµ(z(t 2)) =
dµ(z(t 1)), for any two instants t 1 and t 2 — a property that we
recognize as a generalized form of the Liouville theorem and
that remains valid even when the system’s Hamiltonian has
an explicit time dependence.30
Finally, we define the complex action A,
A(z ∗(t), z(0),t) = S(z) − i~2 [ f (z ∗(t), z(t)) + f (z ∗(0), z(0))] ,
(2.9)
where the first term is the time integral,
 t
S(z) =
L(z)dτ,
(2.10)
0
with the Lagrangian L(z), given by Equation (2.6), being
evaluated over a classical orbit, satisfying (2.8).
As can be seen in (2.9), the complex action A carries
imaginary surface terms, which ensure that this is a welldefined analytical function on both its complex arguments:
z ∗(t) and z(0), and also the time t. The derivatives with respect
to each of these variables are
i ∂ A(z ∗(t), z(0),t) ∂ f (z ∗(t), z(t))
=
,
~
∂zα∗ (t)
∂zα∗ (t)
i ∂ A(z ∗(t), z(0),t) ∂ f (z ∗(0), z(0))
=
,
~
∂zα (0)
∂zα (0)
∂ A(z ∗(t), z(0),t)
= −H(z ∗(t), z(t)).
∂t
(2.11a)
(2.11b)
(2.11c)
The above relations are recognized as the signature of a
properly defined classical action integral.11 Yet it is in terms
of the real quantity S of Eq. (2.10) that our results are most
conveniently expressed.33 Therefore we shall denominate S
the action.
C. Examples of coherent states
In order to illustrate the formalism presented above, we
consider simple examples of coherent states and evaluate some
of their geometrical elements, such as the metric matrix g and
integration measure dµ.
1. Canonical coherent states
Canonical coherent states have their functional definition
given in terms of a superposition of bosonic Fock states
with unrestricted occupation numbers; if the Hilbert space
comprises n modes, then the non-normalized canonical
094106-4
Grigolo, Viscondi, and de Aguiar
J. Chem. Phys. 144, 094106 (2016)
coherent state is defined by
form is
∞

∞ 

 n zαm α 
|z} =
…
√

 |m1, . . . , m n ⟩, (2.12)
mα ! 
m 1=0
m n =0  α=1
where the vacuum |0, . . . , 0⟩ is the reference state and the
length of the vector z equals the number of bosonic modes,
d = n.
Since the modes are assumed to be orthonormal, it follows
immediately from (2.12) that the overlap between canonical
coherent states is
n

{z|z ′} = exp *
zα∗ zα′ + ,
,α=1
–
(2.13)
n

zα∗ zα′ .
(2.14)
α=1
From (2.2), the phase-space metric matrix g is simply the
identity matrix,
gα β (z ∗, z) = δα β ,
(2.15)
which means that canonical coherent states give rise to a flat
phase space and therefore the degrees of freedom are not
geometrically coupled. The normalized measure, as defined
in (2.4), is then trivial,
dµ(z) =
n

d2 z α
α=1
π
.
(2.16)
In what concerns semiclassical trajectory-based methods,
canonical coherent states are undoubtedly the most widely
used type of coherent state. This is so because of the following
well-known homomorphism connecting the ladder operators
(aα† , aα ) of each bosonic mode in (2.12) with the position and
momentum operators (q̂α , p̂α ),
i~
γα
q̂α = √ (aα† + aα ), p̂α = √ (aα† − aα ),
2
γα 2
{m}
N!
m 1! . . . m n !
) 21 
 n−1 m α 
zα  |m1, . . . , m n ⟩,


 α=1
(2.17)
with [aα , a†β ] = δα β and where γα is an arbitrary constant that
sets the appropriate length scale in each mode. The position
representation of |z⟩ is then found to be a multidimensional
Gaussian wavepacket, whose mean position q and mean
momentum p are related to the real and imaginary parts of the
complex vector z, respectively. Furthermore, the dynamics
of (q, p), as obtained from the time-dependent variational
principle, is simply given by Hamilton’s classical equations
of motion in canonical form, the Hamiltonian being the
mean value given by Equation (2.7). This so-called “frozenGaussian representation” provides an obviously suitable
framework for semiclassical applications.31,32,34–40
2. SU(n) bosonic coherent states
The SU(n) bosonic coherent states are suitable for
describing systems in a Fock space comprising n modes
and a fixed total particle number N. Their non-normalized
(2.18)
where the reference state is |0, . . . , N⟩. In (2.18), the primed
summation symbol means that the set of occupation numbers
{m} must satisfy the condition m1 + m2 + · · · + m n = N.
Because of this constraint, the number of entries of the
vector z is one less than the number of modes, d = n − 1.
The overlap is easily evaluated from (2.18) with the help
of the multinomial theorem,

N
n−1

∗ ′

{z|z } = 1 +
zα zα  ;


α=1
′
thus, using (2.1), we identify the function f as
f (z ∗, z ′) =
|z} =
′ (
(2.19)
hence, by means of (2.1), all geometrical aspects of these
coherent states are codified in the function f given by


n−1

f (z ∗, z ′) = N log 1 +
zα∗ zα′  .


α=1
(2.20)
The metric matrix, according to (2.2), is
gα β (z ∗, z) = N
(1 + |z|2)δα β − zα∗ z β
,
(1 + |z|2)2
(2.21)
n−1 ∗
where |z|2 = γ=1
zγ zγ . Clearly, the fixed particle number
condition translates into a geometrical coupling among the
components of the vector z.
Despite the complications introduced by the curved
geometry, the metric’s determinant can be evaluated and
the integration measure, defined in (2.4), is found to be
dµ(z) =
n−1
(N + n − 1)!  d2 zα
.
N!(1 + |z|2)n α=1 π
(2.22)
Recently, semiclassical methods employing SU(n)
bosonic coherent states, including an initial value representation based on classical trajectories in a duplicated phase
space, have been developed and tested with an SU(3) model
Hamiltonian20,21 — in Section IV B, we shall have the
opportunity to revisit that same problem.
3. Spin coherent states
A particularly interesting SU(n) coherent state originates
when the bosonic Fock space has only n = 2 modes. The
(N + 1) states,
|N, 0⟩, |N − 1, 1⟩, |N − 2, 2⟩, . . . , |1, N − 1⟩, |0, N⟩,
can be put into one-to-one correspondence with the wellknown simultaneous eigenstates | J, M⟩ of the angular
momentum operators { Jˆ2, Jˆz } for a fixed total angular
momentum J = N/2. In this way, the non-normalized SU(2)
coherent states9,41 can be expressed as
(
) 21
J

2J
|z} =
z J +M | J, M⟩,
J+M
M =−J
(2.23)
094106-5
Grigolo, Viscondi, and de Aguiar
J. Chem. Phys. 144, 094106 (2016)
which are especially designated as atomic or spin coherent
states. Notice that, in this specific case, the complex vector
label z has dimensions d = 1 and | J, −J⟩ is the reference state.
The overlap {z|z ′} is easily evaluated using the
orthogonality of the | J, M⟩ states. Alternatively, we can simply
set n = 2 and N = 2J in Equation (2.19), thus obtaining
{z|z ′} = (1 + z ∗ z ′)2J ,
(2.24)
which leads to
f (z , z ) = 2J log(1 + z z ).
∗
′
∗ ′
{z|z ′} = det(I N + z † z ′),
(2.25)
According to Equation (2.2), it follows that the phasespace metric g (in this case, a scalar) is simply
2J
.
(1 + |z|2)2
The normalized measure is then found to be
g(z ∗, z) =
(2.26)
(2J + 1) d2 z
.
(2.27)
(1 + |z|2)2 π
The natural topology of the spin coherent state is that of
the surface of a sphere. In practical applications, one typically
writes z in terms of angles θ and φ,
dµ(z) =
z = tan(θ/2)e−iφ ,
determinant.13 Here, the label z is understood as an M × N
matrix, the number of degrees of freedom of the corresponding
phase space being d = M × N. The elements z µα describe the
mixing between occupied and virtual spin orbitals in such
a way that any single-determinantal state not orthogonal to
N
|Φ0⟩ = Â N α=1
|φ•α ⟩ can be represented by |z}.43
The overlap between two distinct Thouless-parametrized
Slater determinants is readily found to be
(2.28)
where θ ∈ [0, π] and φ ∈ [0, 2π). In these coordinates, the
integration measure (2.27) reads
dθdφ
.
(2.29)
4π
Spin coherent states are discussed in more detail in
Section IV A where we investigate a test model consisting of
an SU(2) Hamiltonian.
dµ(θ, φ) = (2J + 1) sin θ
4. SU(n) fermionic coherent states
Fermionic coherent states of the special unitary group
are suitable for describing a number-conserving system of N
fermions which are allowed to occupy a set of n orthonormal
single-particle states (n > N). While ultimately arbitrary,
these underlying single-particle states are often taken to be
eigenstates of the non-interacting part of the Hamiltonian or a
set of Hartree-Fock spin orbitals.9,14 A reference state |Φ0⟩ is
specified by constructing a Slater determinant out of N of such
spin orbitals (e.g., the ones having lowest energies). These are
denoted by |φ•α ⟩ with 1 ≤ α ≤ N, and are said to belong to
the occupied space. The remaining M ≡ n − N spin orbitals,
denoted by |φ◦µ ⟩ with 1 ≤ µ ≤ M, are said to belong to the
virtual space. Then, the non-normalized fermionic coherent
state can written as

N 
M



◦
•
|φα ⟩ +
|φ µ ⟩z µα  ,
(2.30)
|z} = Â N


α=1 
µ=1
where the symbol Â N instructs anti-symmetrization of
the direct product of the N single-particle states |φ•α ⟩
M
+ µ=1
|φ◦µ ⟩z µα , which are sometimes called dynamical
orbitals.16
In the context of quantum chemistry, the coherent state
(2.30) is designated as the Thouless parametrization of a Slater
(2.31)
where I N is the identity matrix of size N × N. Following
the definitions given in Section II, geometrical properties are
determined from the function
f (z ∗, z ′) = log[det(I N + z † z ′)].
(2.32)
From (2.2), one finds (through elementary determinant and
matrix identities) that the phase-space metric can be expressed
as
g µν,α β (z ∗, z) = [(I M + zz †)−1]ν µ [(I N + z † z)−1]α β ,
(2.33)
where I M is the M × M identity matrix and the derivatives
of f were taken with respect to z µα and zν∗ β . Despite the
complicated outlook of (2.33), the classical equations of
motion for the Thouless parameters z, as governed by a
standard many-body Hamiltonian consisting of one- and twobody terms, display a simple structure — see, for instance,
Section III A in Ref. 16.
Finally, the integration measure appearing in the closure
relation (2.3) is found to be
dµ(z ∗, z) = κ[det(I N + z † z)]−n
N 
M

d2 z µα
α=1 µ=1
π
. (2.34)
The normalization constant can be computed by evaluating
the required phase-space integral through a recurrence relation
 N (n+1−γ)!
method; the result is κ = γ=1
(N +1−γ)! .
As one might infer from the above, fermionic coherent
states of this kind are somewhat more intricate than the ones
discussed in the previous examples. And, although the content
of the present paper encompasses all basic ingredients needed
to implement the proposed method in terms of these states,
a proper treatment would nevertheless require introduction of
additional concepts. Therefore, we do not pursue applications
involving this particular set of coherent states in this paper.
III. GENERALIZED CCS METHOD
The CCS method, as originally developed by Shalashilin
and Child25–27 using canonical coherent states, belongs to
the family of multiconfigurational guided-basis methods. Its
characteristic attributes are the non-orthogonality of the basis
set and the use of simple classical mechanics42 to guide
the basis elements, as opposed to more complicated fullvariational approaches like the Gaussian-based version of
the multiconfigurational time-dependent Hartree (G-MCTDH)
method.44–46
Formulation of the method for Gaussian states is fairly
straightforward and the same is true in the generalized context.
In order to better appreciate the additional features that arise in
094106-6
Grigolo, Viscondi, and de Aguiar
J. Chem. Phys. 144, 094106 (2016)
the latter case, we will first present the method in its continuum
form; the discrete version is developed subsequently.
A. The continuum version
We begin by considering the coherent-state decomposition of an arbitrary quantum state,


|ψ⟩ = dµ(z)|z⟩⟨z|ψ⟩ = dµ(z0)|z⟩⟨z|ψ⟩,
(3.1)
which follows from (2.3). It is assumed that z = z(t) is bound
to obey classical equations of motion (2.8). By virtue of
phase-space volume conservation, we are allowed to transfer
the integration measure to the initial instant and conveniently
integrate over initial conditions z0 = z(0), as indicated in the
second equality in (3.1). The derivation of the CCS equations
amounts to finding a solution of the Schrödinger equation,
i~|ψ̇⟩ = Ĥ |ψ⟩,
(3.2)
for |ψ⟩ in the form given by (3.1) with the ansatz,
i
⟨z|ψ⟩ = C(z)e ~ S(z),
(3.3)
where S(z) is the action defined in (2.10). In other words, we
seek an equation of motion for the time-dependent amplitude
C(z) that solves (3.2). Let us make a few observations
regarding this particular choice of solution.
First, all quantities that specify |ψ⟩ — i.e., |z⟩, C(z),
and S(z) — are to be regarded as functions of the initial
conditions z0. Thus, the method is conceived as an initial
value representation from its onset.
Second, we note that it follows from (3.3) that C(z)
depends on the initial state |ψ0⟩ = |ψ(0)⟩ through the relation
C(z0) = ⟨z0|ψ0⟩. In numerical applications, the integral in
(3.1) has to be approximated somehow; the typical procedure
is to sample initial conditions z0 in phase space with the
overlap modulus |⟨z0|ψ0⟩| playing the role of a weight function.
Once the z0’s have been properly sampled, the values of the
corresponding C(z0)’s are uniquely defined.
Third, the motivation behind the factorization of ⟨z|ψ⟩ into
a complex amplitude times an action exponential comes from
a general result of semiclassical theory, according to which
the classical action provides a first-order approximation to the
phase of the quantum state. Since this phase accounts for most
of the wavefunction’s oscillatory behavior, C(z) is expected
to present a rather smooth time dependence, thus facilitating
numerical treatment.
We now proceed to look for a differential equation for
C(z). Taking the time derivative of (3.3) and making use of
the Schrödinger equation, we find (after rearranging terms)

i
i~Ċ(z) = i~⟨ ż|ψ⟩ + ⟨z| Ĥ |ψ⟩ + L(z)⟨z|ψ⟩ e− ~ S(z).
(3.4)
Next, we factor out |ψ⟩ by separating the scalar products
on the right-hand side
 of the equation with the help of the
closure relation 1̂ = dµ(z ′)|z ′⟩⟨z ′|, with z ′ = z ′(t), which
leads to

i
′
i~Ċ(z) = dµ(z0′ )⟨z|z ′⟩∆2 H(z ∗, z ′)C(z ′)e ~ [S(z )−S(z)]. (3.5)
Here we have already shifted the integration measure to the
initial instant [z0′ = z ′(0)] and replaced the ⟨z ′|ψ⟩ that appeared
i
′
under the integral sign for C(z ′)e ~ S(z ).
The coupling ∆2 H(z ∗, z ′) in (3.5) is given by
⟨ ż|z ′⟩
+ H(z ∗, z ′) + L(z),
⟨z|z ′⟩
where the non-diagonal matrix element,
∆2 H(z ∗, z ′) = i~
H(z ∗, z ′) =
{z| Ĥ |z ′} ⟨z| Ĥ |z ′⟩
=
,
{z|z ′}
⟨z|z ′⟩
(3.6)
(3.7)
is an analytical function of z ∗ and z ′ that can be directly
obtained by analytical continuation of classical Hamiltonian
(2.7).
As a last step, we express the first term in (3.6)
as a function of readily computable quantities. Since ⟨z|
1
∗
= e− 2 f (z , z){z|, we observe that
⟨ ż|z ′⟩ { ż|z ′} 1 d
=
−
f (z ∗, z).
⟨z|z ′⟩ {z|z ′} 2 dt
The total time derivative of f (z ∗, z) is simply

d 

∂ f (z ∗, z) ∗ ∂ f (z ∗, z)
d
ż
+
f (z ∗, z) =
ż
α ,
α
dt
∂zα∗
∂zα
α=1
while the remaining term involving { ż| can be written as
d
{ ż|z ′}  ∂ f (z ∗, z ′) ∗
żα ,
=
{z|z ′} α=1 ∂zα∗
owing to the analyticity of {z| on z ∗. Hence, collecting together
the above results and making the necessary substitutions in
(3.6), we find that the coupling may be expressed as
∆2 H(z ∗, z ′) = H(z ∗, z ′) − H(z ∗, z)

d 

∂ f (z ∗, z ′) ∂ f (z ∗, z) ∗
+ i~
−
żα ,
∂zα∗
∂zα∗
α=1
(3.8)
which is an analytic function on z ′.
By integrating equation of motion (3.5), the amplitudes
at time t can be determined from their initial values. Once the
amplitudes are known, we can reconstruct the quantum state
with (3.1), reproduced below in terms of C(z),

i
|ψ⟩ = dµ(z0)|z⟩C(z)e ~ S(z).
(3.9)
Integro-differential equation (3.5) — with ∆2 H(z ∗, z ′)
given by (3.8) — relates directly to the canonical coherent
states version of the CCS method and shares its attractive
characteristics, namely, (i) in the semiclassical regime,
according to the reasons mentioned earlier, the amplitude
C(z) is expected to have a smooth time dependence; (ii)
because of the coherent-state overlap ⟨z|z ′⟩, the z ′ integral is
localized around z;47 and (iii) the coupling between amplitudes
of different basis elements is not only sparse but also nondiagonal, since the integrand is identically zero when z ′ = z.
As a final remark, we should note that, if one performs a
series expansion of ∆2 H(z ∗, z ′) for small |z ′ − z|, one finds that
this series begins with a second-order term. In the generalized
coherent state case however, unlike the specific situation
for canonical coherent states, this does not coincide with
094106-7
Grigolo, Viscondi, and de Aguiar
J. Chem. Phys. 144, 094106 (2016)
the second- and higher-order terms in the Taylor series of
H(z ∗, z ′).
B. The discrete version
Here we re-derive the CCS equations using a discrete
set of coherent states as basis, as opposed to the continuous
set employed in Sec. III A. We shall find that the discrete
formulas do not differ from their analogue expressions in the
canonical coherent state case. This is due to the fact that all
information concerning distinct coherent-state geometries is
codified in a number of key quantities, namely, the overlap,
the phase-space metric, and the classical equations of motion.
Each of these quantities participates in the same way in
the working equations, regardless of the particular type of
coherent state being used — incidentally, a most desirable
feature for programming purposes, for it means that the core
subroutines of the method are essentially independent of
geometry. Nevertheless, for the sake of consistency, we must
review the discrete formulas, for they are the ones actually
used in practice.
1. Unitary propagation
The first step towards a discrete unitary formulation of
the generalized CCS method is the assumption that one is able
to write down a closure relation by employing a finite number
of basis elements as
1̂ =
M 
M

|z j ⟩Λ j k ⟨z k |,
(3.10)
j=1 k=1
where M is the size of the basis set.48 It suffices that this
closure relation represents the identity operator only on the
dynamically accessible part of the phase space and that it holds
only during the time interval upon which the propagation takes
place.
In order to properly represent the identity, the matrix Λ
in (3.10) must satisfy the relation,

δ jk =
Ω jl Λl k ,
(3.11)
l
where
notable weakness of time-dependent methods formulated
with non-orthogonal basis sets.49 In the event that Ω
becomes singular, one should take appropriate measures
before resuming the propagation. In this regard, a particularly
interesting protocol has been proposed by Habershon.50 The
“singularity problem,” though, did not occur in the simple
applications considered in this paper.
Upon these considerations, we introduce a discrete set
of M amplitudes C j = C(z j ) as well as their corresponding
action phases S j = S(z j ), writing
i
⟨z j |ψ⟩ = C j e ~ S j .
(3.14)
Next, we proceed exactly as in Section III A. The equation of
motion in the discrete unitary case is then readily found to be

i
i~Ċ j =
Ω j k ∆2 H j k Λk l Cl e ~ (Sl −S j ),
(3.15)
k,l
with coupling matrix given by
∆2 H j k = H(z ∗j , z k ) − H(z ∗j , z j )
d  ∂ f (z ∗, z )

∂ f (z ∗j , z j ) 

j k
 ż ∗ .
+ i~
−
(3.16)
 ∂z ∗
∗
∂z
 jα
jα
jα
α=1 
In practice, matrix Λ is never explicitly constructed;
rather, one introduces a set of auxiliary amplitudes D
= (D1, D2, . . . , D M ) which are related to the coefficients
C = (C1,C2, . . . ,CM ) according to

i
(3.17)
Ωl k Dk e ~ (S k −Sl ) = Cl .
k
Thus, at every time step, D is obtained from C by means of
the above intermediate equation — an operation that requires
solving a linear system of size M. Then, equation of motion
(3.15) can be recast as


i
i~Ċ j =
Ω j k ∆2 H j k e ~ (S k −S j ) Dk ,
(3.18)
k
while the quantum state is expressed in terms of amplitudes
D as

i
|ψ⟩ =
(3.19)
|z k ⟩Dk e ~ S k .
k
Ω j k = ⟨z j |z k ⟩
(3.12)
is the overlap matrix. This guarantees that we have a welldefined discrete coherent-state decomposition,

|ψ⟩ =
|z j ⟩Λ j k ⟨z k |ψ⟩.
(3.13)
j,k
In other words, Ω must be sufficiently well-conditioned so
that expressions involving its inverse Λ remain numerically
stable during the entire propagation. Therefore, the basis-set
initial conditions must be sampled in such a way as to ensure
that this requirement is fulfilled. One such procedure, that
results in a well-conditioned overlap matrix (at initial time),
is described in the Appendix.
Yet, nothing prevents that an initially well-conditioned
overlap matrix becomes singular at some later time — a
The propagation scheme comprised by Eqs. (2.8), (2.6), and
(2.10) together with Eqs. (3.17)–(3.19) represents the standard
form of the generalized coherent-state method proposed here.
It can be shown to preserve normalization — given by the

sum k Ck∗ Dk — as long as the overlap matrix remains
sufficiently well-conditioned. In addition, it preserves total
energy (for time-independent Hamiltonians) as long as the
identity operator can be resolved in terms of the basisset elements, though this situation is hardly achieved in
multidimensional problems. It has been pointed out that
energy conservation is closely related to the accuracy of
the CCS method.50 Thus, by monitoring total energy, one can
make an “on-the-fly” diagnosis as regards to the quality of
the CCS results; indeed we observe in our simulations that
deviations from the exact quantum solution are accompanied
by fluctuations in total energy.
094106-8
Grigolo, Viscondi, and de Aguiar
J. Chem. Phys. 144, 094106 (2016)
2. Non-unitary case
It may be convenient — particularly when the system
under study has only one or two degrees of freedom — to
attempt a more straightforward discretization of the closure
relation, as
1̂ ≈
M

solution. In other situations, it may provide a reasonable
approximation for very short times.
IV. APPLICATION EXAMPLES
A. Condensate in a double-well potential
|z k ⟩λ k ⟨z k |,
(3.20)
k=1
with λ k approximating the integration measure dµ(z k ) at each
phase-space point.
The equation of motion for C in this case can be obtained
at once from (3.15) by setting Λk l = λ k δ k l ,

 
i
(3.21)
i~Ċ j =
λ k Ω j k ∆2 H j k e ~ (S k −S j ) Ck .
k
Similarly, the quantum state in this case is given by

i
|ψ⟩ =
λ k |z k ⟩Ck e ~ S k .
(3.22)
k
This propagation scheme is computationally less demanding
than the one discussed in Sec. III B 1 — if the basis-set
size is kept the same — since there is no need to solve a
linear system at each time step. On the other hand, a larger
basis set (usually constructed as a grid in phase space) may
be necessary to converge the results if (3.21) is employed.
Moreover, as discussed by Shalashilin and Child,27 the direct
discretization of Equation (3.5) does not preserve the unitarity
of an exact quantum time evolution, that is, the norm of
the propagated quantum state is not automatically conserved,
meaning that results must be normalized on output.
C. Classical propagation
To end this section, we discuss the particular situation
whereupon a single coherent-state basis element is used to
describe the system,
i
|ψ⟩ = |z⟩e ~ S(z).
(3.23)
As can be seen from the above equation, this scheme only
applies if the quantum state to be propagated is a coherent state,
that is, |ψ0⟩ = |z0⟩. The form of the approximated |ψ⟩, with an
action phase, can be derived from the working equations of
Sec. III B by setting the basis-set size M = 1, in which case
we find that Ċ = 0 and hence C(t) = C(0) = 1.
We shall denominate the primitive method defined by
(3.23), together with (2.8), (2.6), and (2.10), as the classical
propagation scheme, in view of the fact that only classical
ingredients are present in its formulation. It serves as a reference method, against which more sophisticated approaches,
such as those described earlier, may be confronted —
which is useful, for example, in order to identify nonclassical behavior (defined in this sense), as in Section IV B.
It should be pointed out that if the Hamiltonian of the
system is such that application of the time-evolution operator
maps one coherent state onto another — or, more formally,
when the Hamiltonian is an element of the Lie algebra associated with the set of coherent states under consideration —
then the classical propagation scheme actually gives the exact
As a first application, we consider a simplified model
for the dynamics of an N-particle Bose-Einstein condensate
trapped in a double-well potential, where individual bosons
interact through contact forces, that is, with an interacting
potential V (x, x′) ∝ δ(x − x′). This model has been discussed
in detail in several works;51–53 here we briefly sketch its main
ideas.
In the two-mode approximation, it is assumed that
only the single-particle ground state and first excited state
of the double-well potential, as obtained from first-order
perturbation theory, have a significant occupation, so that the
dynamics of the system is restricted to these two levels. This
should be a good approximation if one seeks to describe the
low temperature regime, wherein most of the particles are
expected to be occupying the ground state — see Ref. 54 for
a discussion on the validity of the two-mode approximation
to the double-well problem.
Since the number of bosons is preserved, the particle
number operator N̂ is a constant of the motion. By means of
a well-known homomorphism between the algebra su(2) and
the bosonic creation and annihilation operators, the bosonic
dynamics can be described in terms of three independent
angular momentum operators — Schwinger’s pseudo-spin
operators55 — with total angular momentum J = N/2.
In terms of these operators the Hamiltonian of the model
is56
Ĥ = Ω Jˆx +
2 χ ˆ2
J ,
N −1 z
(4.1)
where the tunneling rate Ω equals the energy difference
between the two occupied single-particle states and the selfcollision parameter χ is proportional to the interaction strength
of bosons located in the same potential well.
Applying definitions (2.7) and (2.23) to (4.1), the classical
Hamiltonian becomes57
H(z ∗, z) =
NΩ z + z ∗
N χ (1 − z ∗ z)2
+
,
2 1 + z∗z
2 (1 + z ∗ z)2
(4.2)
in which we discarded, without any loss, a constant term.
From identity (2.8), the equation of motion for z(t) is
i ż =
Ω
z(z ∗ z − 1)
(1 − z 2) + 2 χ
.
2
1 + z∗z
(4.3)
Notice that the (N − 1)−1 scaling of the self-collision
parameter in (4.1) was specifically chosen so that equation of
motion (4.3) is independent of particle number N, therefore
representing a well-defined classical limit of the system when
N → ∞. This also means that, by employing the CCS method
with a fixed set of classical trajectories, we can obtain quantum
solutions for different particle number regimes.
094106-9
Grigolo, Viscondi, and de Aguiar
J. Chem. Phys. 144, 094106 (2016)
1. Non-unitary CCS with spin coherent states
We shall examine the quantum dynamics of Hamiltonian
(4.1) from the perspective of a non-unitary implementation
of the CCS method. Following the procedure outlined in
Section III B 2, a direct discretization of expressions (3.5) and
(3.9) is performed by replacing the integrals over phase-space
variables at t = 0 by finite summations over a discrete set of
initial conditions.
The accuracy of this straightforward discretization
procedure is strongly dependent on the selection of initial
values over the phase space. In the particular case of a
regular grid, results are found to be extremely sensitive to
characteristics such as spacing, placement, number of points,
and, most importantly, the choice of grid variables.
For the subsequent numerical calculations, we consider
a regular grid in the phase-space coordinates (η, ζ), which
are related to the SU(2) angular coordinates (θ, φ) defined in
Equation (2.28), according to
η = θ,
(4.4)
ζ = φ sin θ.
(4.5)
This choice aims at reducing the inherent difficulties imposed
by the spherical topology of spin coherent states at θ ≈ 0 or π.
In terms of this new set of coordinates, integration
measure (2.29) takes a simpler form,
2J + 1
dηdζ,
(4.6)
4π
which, unlike expressions (2.27) and (2.29), is independent
of the integration variables. Hence, the λ k ’s that appear in the
non-unitary evolution equation (3.21) and also in (3.22) are
all equal to the same time-independent λ given by
dµ =
2J + 1
∆η∆ζ,
(4.7)
4π
where ∆η and ∆ζ are the grid spacings in the η and ζ
directions, respectively.
The initial quantum state to be propagated is chosen to
be the spin coherent state |ψ0⟩ = |z ′⟩ with
λ=
z ′ = tan(π/8),
specified by angular coordinates (θ ′, φ ′) = ( π4 , 0). As it is
usually assumed in most guided-basis methods, the most
important dynamical contributions — at least for short-time
propagation — are expected to arise from classical trajectories
initially located at the same phase-space region occupied by
|ψ0⟩. Therefore, the grid of initial conditions is intentionally
centered at the point (θ ′, φ ′), as illustrated on the left panel of
Fig. 1 for the particular case of N = 100 particles — the right
panel portrays the classical dynamics in phase space. Notice
that the classical trajectory with initial value z ′ is located close
to a stable equilibrium point and within the boundaries of a
separatrix of motion.
With the purpose of quantifying the agreement between
the CCS-propagated quantum state |ψccs(t)⟩ and the exact
result |ψexact(t)⟩ — obtained via diagonalization of (4.1) in the
angular momentum basis — we shall compute the fidelity,
F (t) = |⟨ψexact(t)|ψccs(t)⟩|.
(4.8)
The fidelity reflects the physical similarity between two states:
its values are restricted to the interval 0 ≤ F ≤ 1, with
the maximum value corresponding to physically identical
quantum states.
We wish to analyze the behavior of F (t) as the number
of particles in the system is changed. There is a subtlety
involved, though: due to its fundamental property of minimal
uncertainty,7,8 the coherent state |z⟩ represents a quantum state
with maximal localization in phase space around the point z.
Moreover, for spin coherent states, the linear dimensions of
the phase-space region
√ effectively occupied by |z⟩ decrease
proportionally to 1/ N.58 Thus, in order to make a meaningful
comparison amongst trajectory-based propagations carried out
at different particle number regimes, one must account for the
“shrinking” of the initial state |ψ0⟩ = |z ′⟩ as N grows larger.
Therefore, in our simulations,
the grid spacing was reduced
√
proportionally to 1/ N whilst the number of initial conditions
in each run was roughly unchanged.
Fig. 2 shows F as a function of time for various values
of N. Notice that, even for long times, the non-unitary results
remain accurate (F > 0.96). On the other hand, it is clear
that the fidelity tends to decrease over time, evidencing the
FIG. 1. Left panel: regular grid of classical initial conditions with approximately three hundred initial values for the case of N = 100 particles. The red dot
indicates the center of the grid, whose angular coordinates are (θ ′, φ ′) = ( π4 , 0). This phase-space point represents the label of the initial coherent state |ψ 0⟩ = |z ′⟩.
Right panel: some examples of classical trajectories over the spherical phase space. The red orbit corresponds to the initial value z ′. The Cartesian coordinates
used to plot these trajectories are related to the coherent-state angular variables according to (x 1, x 2, x 3) = (sinθ cosφ, sinθ sinφ, −cosθ).
094106-10
Grigolo, Viscondi, and de Aguiar
FIG. 2. Fidelity as a function of the dimensionless time Ωt with parameters
Ω = 1.0, χ = 1.0 for different particle number regimes. Each run employed
roughly three hundred trajectories — the number of points varies slightly for
different runs due to the cropping of grid borders.
accumulation of numerical errors in the CCS dynamics. In
part, this inaccuracy stems from the non-unitary nature of the
discretization procedure. However, to a large extent, numerical
error arises due to dispersion of the classical trajectories over
phase space, resulting in an incomplete description of the
quantum system.
Also note that the fidelity is consistently higher for larger
particle numbers. Since the number of classical trajectories
was kept nearly constant for all the runs and the range of
the initial condition grid was adjusted so as to reflect the
width of the initial quantum state, we may conclude that
the CCS method is better suited for describing the dynamics
in the many-particle regime, which we recognize as the
semiclassical regime for this problem. Indeed, it is well
known that the limiting case N → ∞ (or J → ∞, as it can be
interpreted here) coincides with the classical limit of quantum
mechanics.59
In this sense, the CCS method, though in principle formulated as an exact method, when numerically implemented —
and therefore subjected to inevitable practical limitations —
should be regarded more as a semiclassical technique than
a genuine quantum approach, an assertion that can be made
with respect to both its non-unitary and unitary versions.
Finally, we point out that, for the system under
consideration, the computational cost of the method is
insensitive to particle number. Hence, for very large values
of N, the computational resources needed for evaluating the
exact quantum dynamics by diagonalization of Hamiltonian
(4.1) will certainly exceed the analogous requirements of the
non-unitary CCS method.
B. Condensate in a triple-well potential
Next, we consider a simplified model describing an Nparticle Bose-Einstein condensate trapped in a symmetric
triple-well potential where individual bosons are once more
J. Chem. Phys. 144, 094106 (2016)
assumed to interact by contact forces; the main ideas involved
are as follows: the triple-well trapping potential, under suitable
conditions, can be approximated by an harmonic expansion
around each of its three (symmetrically located) minima. The
three-fold degenerate fundamental states of this approximated
problem can be determined without difficulty. It is then
assumed that the dynamical regime is such that the energy
eigenspace spanned by these three local modes is sufficiently
isolated from the rest of the single-particle spectrum, so
that at low temperatures, they alone provide an adequate
description of the system. For more details on the derivation
and particularities of this model, see Refs. 60 and 61.
Let a1, a2, and a3 denote the annihilation operators associated with the aforementioned fundamental single-particle
modes related to the locally approximated wells. In terms of
these bosonic operators, the “three-mode approximation”61 to
the Hamiltonian is


χ
aγ† aγ† aγ aγ , (4.9)
Ĥ = Ω
aα† a β +
N
−
1
1≤γ ≤3
1≤α,β ≤3
where Ω is the tunneling rate, describing hops between
adjacent wells, and χ is the collision parameter, that controls
the strength of two-body interactions within the same well.62
Owing to particle number conservation, this system is
suitably described in terms of SU(3) bosonic coherent states
|z⟩ = |z1, z2⟩, which represent a particular case of the coherent
states discussed in Section II C 2. Using definition (2.18)
together with (4.9), we find from (2.7) that the classical
Hamiltonian is
(z ∗ z2 + z2∗ z1 + z1∗ + z1 + z2∗ + z2)
H(z ∗, z) = NΩ 1
1 + z1∗ z1 + z2∗ z2
+Nχ
(z1∗ z1)2 + (z2∗ z2)2 + 1
.
(1 + z1∗ z1 + z2∗ z2)2
(4.10)
From (2.8), it follows that the equations of motion are
i ż1 = Ω(1 + z1 + z2)(1 − z1) −
2 χz1(1 − |z1|2)
,
1 + |z1|2 + |z2|2
(4.11a)
2 χz2(1 − |z2|2)
.
(4.11b)
1 + |z1|2 + |z2|2
Similarly to the double-well model, we have deliberately
tuned the collision parameter χ with a (N − 1)−1 factor,
thereby making the classical dynamics independent of particle
number; in this way, the classical system is well-defined in
the limit N → ∞.
i ż2 = Ω(1 + z1 + z2)(1 − z2) −
1. Unitary CCS with SU(3) bosonic coherent states
The classical system defined in (4.11) has three
dynamically equivalent invariant subspaces, specified by the
constraints: z1 = z2, z1 = 1, and z2 = 1. Let us concentrate on
the first of these (z1 = z2) and refer to it as the Γ1 subspace.
Now, consider the set of operators b1, b2, and b3, defined
by the canonical transformation,
b1 = √12 (a1 + a2),
(4.12a)
b2 = a3,
b3 = √12 (a1 − a2).
(4.12b)
(4.12c)
094106-11
Grigolo, Viscondi, and de Aguiar
J. Chem. Phys. 144, 094106 (2016)
It can be demonstrated that Γ1 is an SU(2) subspace whose
two single-particle modes are associated with operators b1 and
b2.61 As a consequence, under the classical approximation
described in Section III C, any SU(3) coherent state initially
located in Γ1 will always display zero occupation of the b3
mode; as a matter of fact, the expectation value,
N (z1∗ − z2∗)(z1 − z2)
,
⟨z|b†3b3|z⟩ =
2 1 + z1∗ z1 + z2∗ z2
(4.13)
is identically null when z1 = z2.
This conclusion, however, does not apply to the actual
quantum problem: even if the initial state |ψ0⟩ has null
occupation in the b3 mode, this situation will not be preserved
as the system evolves in time. It is precisely this “non-classical
behavior” that we wish to describe using the SU(3) CCS
method.
Let us then consider the initial state |ψ0⟩ = |z1′ , z2′ ⟩, with
z1′ = z2′ = √12 tan(π/8),
and thus located on the classical invariant subspace Γ1. We
shall propagate this state employing the unitary method
developed in Section III B 1 and compute the occupation
Q(t) = ⟨ψ(t)|b†3b3|ψ(t)⟩ of the single-particle mode associated
with b3. Following (3.19), this function will be calculated
according to
Q=
 {z |b†b |z } 
j 3 3 k
 ei(S k −S j ),
D ∗j Dk Ω j k 
{z
|z
}


j
k
j,k


(4.14)
where the non-diagonal matrix elements between square
brackets can be obtained by analytic continuation of (4.13).
We shall also monitor the total energy of the system,

E=
D ∗j Dk Ω j k H(z ∗j , z k )ei(S k −S j ),
(4.15)
j,k
as a means to probe the quality of our results.
In order to construct the basis set, it is necessary to choose
adequate sampling variables. In the present case, we opt for
angular variables (θ 1, φ1, θ 2, φ2) defined by
z1 = tan(θ 1/2)e−iφ1,
(4.16)
.
(4.17)
z2 = tan(θ 2/2)e
−iφ 2
The initial conditions z(0) are then randomly sampled around
z ′ = (z1′ , z2′ ) from Gaussian probability distributions expressed
in terms of these angular variables, with each angle being
individually selected: for example, θ 1(0) and φ1(0) are selected
according to
∼ e−[θ1(0)−θ1] /2σ θ 1, ∼ e−[φ1(0)−φ1] /2σ φ 1.
′ 2
′ 2
(4.18)
Notice that the widths of these distributions are adjustable
parameters of the method. The actual sampling procedure —
which also comprises specific criteria for accepting and
neglecting candidate basis elements — is somewhat involved
and further details are reserved to the Appendix. Once the
basis-set initial conditions are known, the amplitudes C and
D of Eqs. (3.14) and (3.17) can be initialized and propagation
of (3.18) may be started.
We studied the case of N = 100 trapped bosons, with
tunneling rate Ω = −1.0 and collision parameter χ = −1.0.
The widths of the Gaussian distributions — exemplified in
π
π
for both θ 1 and θ 2 and to 10
for both
(4.18) — were set to 20
φ1 and φ2.
In order to visualize whether results converge as the
basis-set size M is increased, simulations have been performed
with different number of orbits while maintaining all other
parameters — including sampling widths and random number
generator seed63 — fixed.
Results obtained with the unitary SU(3) CCS method
and by a direct diagonalization of Hamiltonian (4.9) in the
bosonic Fock basis are compared in Fig. 3. On the left panel,
the occupation of the classically inaccessible b3 mode is
displayed for different basis-set sizes. On the right panel,
total energy expectation value (4.15) for each run is shown as
well as the corresponding behavior of basis-set conditioning
factors (see the Appendix).
It is found that CCS results improve as more basis
elements are included in the propagation, as expected —
although there is some anomalous behavior, with the M = 20
run performing better than the M = 30 case. This is most
likely due to the basis-set sizes involved, which are probably
too small for a uniform convergence to be observed.
We also note that the quality of CCS results is intimately
connected to conservation of total energy: the more accurate
runs are also those in which energy is better conserved.
FIG. 3. Left panel: occupation of the classically inaccessible b 3 mode as a function of dimensionless time |Ω|t for N = 100, Ω = −1.0, χ = −1.0, and different
basis-set sizes M ; inset: zoom on the time interval 9.0 ≤ |Ω|t ≤ 10.0. Right panel: energy fluctuation as a function of |Ω|t for the same runs shown on the left
panel; inset: basis-set conditioning factor ε(t) for the corresponding curves.
094106-12
Grigolo, Viscondi, and de Aguiar
The most computationally expensive run, with M = 60
basis elements, led to an initial basis-set conditioning
factor ε(0) ≈ 8.8 × 107, a value that diminished during the
propagation. This behavior was observed in all cases —
as shown in the inset of the right panel of Fig. 3 — indicating
spread of the classical trajectories. For this particular run
(M = 60), excellent agreement between the CCS and exact
solutions is found, except for slight discrepancies at |Ω|t & 8.0.
It should be noted that the propagation with M = 60
— which is quite a modest number of basis elements for a
problem with two degrees of freedom — not only proved to be
accurate64 but also considerably faster than solving the 5151dimensional quantum eigensystem. This gain in efficiency, at
least for the model discussed here, drastically increases as
more particles are added to the system, since the dimension of
the SU(3) bosonic Fock space, given by 12 (N + 1)(N + 2),
grows rapidly with N. At the same time, the system
becomes more classical as more particles participate in the
dynamics, therefore making a trajectory-based approach more
inviting.
Yet, despite these advantages, a reliable and positive way
to check for the convergence of the CCS method still remains
to be developed — the criterion of energy non-conservation
being only indicative of accuracy loss during propagation.
This represents a critical shortcoming when studying more
complex systems for which there are no exact numerical
results — or experimental data — available for comparison;
for in those cases no precise statements about the system could
be made based on CCS results alone.
V. SUMMARY AND CONCLUSIONS
In this work, we have formulated a multiconfigurational,
trajectory-guided quantum propagation method whose distinctive quality consists in employing generalized coherent states
as basis elements. In this sense, the technique is seen as a
natural extension of the coupled coherent states method of
Shalashilin and Child25–27 whereupon frozen Gaussians are
replaced by more general configurations; at the same time,
the main features of the original CCS are retained: quantum
amplitudes obey an integro-differential equation with sparse
coupling and present a smooth time dependence, owing to
their oscillatory behavior being partially compensated by the
classical motion of the basis elements and their action phases.
As pointed out in Section II, no deep understanding of
group-theory concepts is necessary, neither to derive the basic
equations of the method nor to implement it numerically —
in fact, we have seen that all geometrical quantities that
enter the formulas can be straightforwardly evaluated from
the coherent-state overlap function alone and that the discrete
versions of the working equations do not differ in overall
structure from their analogue expressions in the original CCS
approach.
In Section III, three versions of the method have been
devised: continuum, non-unitary, and unitary. The continuum
version most evidently displays the novel elements due to
the non-Euclidean geometry associated with the generalized
coherent-states and it serves primarily as a starting point
for a number of possible analytical approximations. The
J. Chem. Phys. 144, 094106 (2016)
non-unitary version, in turn, might be understood as a direct
attempt to reproduce the continuum formulas by reducing
phase-space integrals into finite sums. The unitary version,
meanwhile, is the standard form of the method, being the most
adequate for the majority of practical applications. We have
also briefly discussed one weakness of this last version of the
method which is the ill-conditioning of the overlap matrix. In
that respect, we note that the basis-set sampling procedure that
we propose — detailed in the Appendix — is more appealing
than the techniques used in previous CCS applications, for it
assures an initially well-conditioned overlap matrix and hence
stability of amplitude equation (3.18), at least for short time
propagation.
In Section IV, we have illustrated the general aspects of
the proposed approach with applications to simple models
of many-boson systems, described in terms of SU(2) and
SU(3) bosonic coherent states. One particular aspect, namely,
the choice of appropriate phase-space variables for either
sampling of initial conditions or grid construction, is found
to be of great importance — coordinates chosen for these
purposes must somehow reflect the intrinsic geometry of
the coherent state under consideration, otherwise only a
poor representation of the system is obtained. Moreover,
the accuracy and efficiency of the method were established
by comparing results against exact quantum calculations.
Excellent agreement was observed and, in the triple-well
system examined in Section IV B, small discrepancies in
the results were found to be correlated with energy nonconservation. These tests also allowed us to gain some insight
regarding the domain of applicability of our formulas. A
careful analysis of the fidelity in the double-well system
studied in Section IV A revealed that the CCS approach
is best suited for describing the regime of large particle
number, also identified as the semiclassical regime of that
problem. This conclusion can in fact be extended for SU(n)
bosonic systems in general. The same goes for the observations
concerning the computational advantages of the method (over
standard matrix diagonalization) when the number of particles
is large; assertions which were made within the context of the
triple-well problem.
Though we have exemplified the use of the method with
SU(2) and SU(3) bosonic Hamiltonians, we emphasize that the
formulation presented is by no means limited to these types of
systems — in addition to the SU(n) fermionic coherent states
discussed in Section II C 4, other possible coherent-state
parametrizations in terms of which the method could be
readily implemented were referenced throughout Section I.
Furthermore, the conclusion put forth in Section IV, namely,
that propagation of quantum states with the help of classically
guided basis sets is better suited for describing systems in
their semiclassical regime, is certainly expected to hold for
any coherent-state representation that might be employed —
provided one correctly interprets what “semiclassical regime”
means in each case.
A few words regarding the purely computational aspects
of the method are in place. In the numerical calculations
reported in this paper, we have adopted the simplest
possible strategy of implementation, which consists in
propagating simultaneously the coherent-state trajectories and
094106-13
Grigolo, Viscondi, and de Aguiar
their corresponding amplitudes. But since trajectories evolve
independently, an alternative approach would be to propagate
them separately, storing the coherent-state coordinates at predetermined time intervals, and use this information afterwards
in order to propagate the fully coupled amplitude equation —
perhaps employing interpolation algorithms if coordinates are
needed at intermediate instants.
This preliminary propagation of classical orbits could,
of course, take full advantage of parallelization, especially
if single basis elements are propagated with ease. It is
however in those cases in which the very propagation of
individual trajectories is a computationally demanding task —
either because the system has an extremely large number of
degrees of freedom or because no analytical expression for the
Hamiltonian is available — that this “two-stage” procedure
would most definitely be the strategy of choice. In addition,
it would allow for more sophisticated sampling techniques,
since it would then be possible to select trajectories based on
knowledge of their entire story and not relying just on their
initial proximity to the quantum state. Thus, for instance, one
could identify orbits which, though unimportant at early stages
of the propagation, give a significant contribution later on.
As a final and interesting remark, we note that the
multiconfigurational Ehrenfest method,66,67 as specifically
designed for “on-the-fly” non-adiabatic dynamics, can be
obtained at once from the presented formalism as the
particular case wherein each basis-set element is taken to
be a composite coherent state consisting of a canonical part
and an SU(n) part with the particle-number parameter N
set to unity.65 In this picture, different “diabatic” potential
surfaces would be represented by the n single-particle states
composing the SU(n) coherent state; further extension to BornOppenheimer “adiabatic” energy surfaces could be achieved
without difficulty.
ACKNOWLEDGMENTS
This work was supported by FAPESP under Project
Grant Nos. 2008/09491-9, 2011/20065-4, 2012/20452-0, and
2014/04036-2. In addition, M.A.M.A. acknowledges support
from CNPq. Numerical calculations made use of routines
from the GNU Scientific Library.68
APPENDIX: BASIS SET SAMPLING
The basis set sampling procedure that we propose for
the unitary CCS method (Section III B 1) is outlined here. It
applies to any type of coherent state once two geometrydependent ingredients are provided: adequate sampling
coordinates q = f (z) — with a known inverse z = f −1(q) —
and a weight distribution function w(q) according to which
these coordinates are to be randomly selected. The procedure
assumes that the initial state is a coherent state, i.e., |ψ0⟩ = |z ′⟩,
in which case the initial coherent-state sampling coordinates
q ′ = f (z ′) must also be supplied.
The sampling strategy follows a very simple “one-by-one”
protocol, which draws inspiration from previously developed
basis set conditioning techniques.50 The procedure amounts
to four steps.
J. Chem. Phys. 144, 094106 (2016)
(1) Take |z ′⟩ (the initial state itself) as the first basis element.
(2) Using the appropriate sampling coordinates q and weight
function w(q), randomly select a new basis element
z j = f −1(q j ) and temporarily add |z j ⟩ to the basis set.
(3) Compute the overlap matrix Ω and evaluate its conditioning factor ε = λ max/λ min, where λ max and λ min are the
largest and smallest eigenvalues of Ω, respectively.69
(4) If ε is less than some threshold value ε lim, accept |z j ⟩
permanently adding it to the basis set, whose size increases
by 1. Else, discard the selected basis element, in which
case the basis-set size does not change. In either case,
return to step (2).
The above procedure is then iterated until either a desired
basis-set size is achieved or saturation occurs, meaning that
the algorithm is unable to select a new |z j ⟩ that satisfy the ε
threshold condition (a certain maximum number of attempts
may be specified). How fast saturation takes place will depend
upon the system’s dimensionality, the threshold value ε lim,
the coherent-state parameters, and the details of the sampling
distribution w(q). Typically, we take ε lim ∼ 108–1012 and use a
predetermined basis-set size below saturation, thus ensuring a
reasonably well-conditioned overlap matrix at initial time and
hence the stability of the propagation (at least for sufficiently
short times).
This sampling protocol requires the eigenvalues of the
overlap matrix to be computed at every iteration. We point
out, however, that this does not compromise the method’s
efficiency since the sampling is performed only once, before
the actual propagation. Moreover, the overlap matrix typically
does not grow too large; this assertion holds even for
multidimensional systems, as long as the sampling distribution
is kept sufficiently localized around the initial-state coordinate
z ′ from where the most relevant contributions to the initial
value representation formula are expected to originate. Finally,
note also that the initial state |z ′⟩ is always included in the
basis set; this is crucial for accuracy of short-time results and
also secures that the initial norm is unity, regardless of how
the remaining basis elements are distributed in phase space.
1N. Makri, “Time-dependent quantum methods for large systems,” Annu.
Rev. Phys. Chem. 50(1), 167 (1999).
2W. H. Miller and K. A. White, “Classical models for electronic degrees
of freedom: The second-quantized many-electron Hamiltonian,” J. Chem.
Phys. 84(9), 5059 (1986).
3A. Kirrander and D. V. Shalashilin, “Quantum dynamics with fermion
coupled coherent states: Theory and application to electron dynamics in laser
fields,” Phys. Rev. A 84(3), 033406 (2011).
4More recently, Grossmann et al. [F. Grossmann, M. Buchholz, E. Pollak,
and M. Nest, “Spin effects and the Pauli principle in semiclassical electron
dynamics,” Phys. Rev. A 89(3), 032104 (2014)] have investigated, in a semiclassical context, whether propagation with antisymmetrized basis states is
essential for the description of electron scattering.
5H. Feldmeier and J. Schnack, “Molecular dynamics for fermions,” Rev. Mod.
Phys. 72(3), 655 (2000).
6P. Kramer and M. Saraceno, Geometry of the Time-Dependent Variational
Principle in Quantum Mechanics (Springer-Verlag, New York, 1981).
7R. Delbourgo, “Minimal uncertainty states for the rotation and allied
groups,” J. Phys. A: Math. Gen. 10(11), 1837 (1977).
8R. Delbourgo and J. R. Fox, “Maximum weight vectors possess minimal
uncertainty,” J. Phys. A: Math. Gen. 10(12), L233 (1977).
9W.-M. Zhang, D. H. Feng, and R. Gilmore, “Coherent states: Theory and
some applications,” Rev. Mod. Phys. 62(4), 867 (1990).
10T. Van Voorhis and D. R. Reichman, “Semiclassical representations of
electronic structure and dynamics,” J. Chem. Phys. 120(2), 579 (2004).
094106-14
Grigolo, Viscondi, and de Aguiar
J. Chem. Phys. 144, 094106 (2016)
11T. F. Viscondi, A. Grigolo, and M. A. M. de Aguiar, “Semiclassical propa-
37E. Kluk, M. F. Herman, and H. L. Davis, “Comparison of the propagation
gator in the generalized coherent-state representation,” e-print arXiv:1510.
05952/quant-ph (submitted).
12Their discussion is based on a rough extension of Solari’s semiclassical
propagators [see (a) H. G. Solari, “Glauber’s coherent states and the
semiclassical propagator,” J. Math. Phys. 27(5), 1351 (1986); (b) H. G.
Solari, “Semiclassical treatment of spin system by means of coherent states,”
J. Math. Phys. 28(5), 1097 (1987)]—a rigorous derivation of the generalized
coherent-state propagator can be found in a recent work.11
13D. J. Thouless, “Stability conditions and nuclear rotations in the HartreeFock theory,” Nucl. Phys. 21, 225 (1960).
14E. Deumens and Y. Öhrn, “Time-dependent dynamics of a determinantal
state,” J. Mol. Struct.: THEOCHEM 199, 23 (1989).
15E. Deumens, A. Diz, H. Taylor, and Y. Öhrn, “Time-dependent dynamics of
electrons and nuclei,” J. Chem. Phys. 96(9), 6820 (1992).
16E. Deumens, A. Diz, R. Longo, and Y. Öhrn, “Time-dependent theoretical
treatments of the dynamics of electrons and nuclei in molecular systems,”
Rev. Mod. Phys. 66(3), 917 (1994).
17T. Suzuki, “Classical and quantum mechanical aspects of time-dependent
Hartree-Fock trajectories,” Nucl. Phys. A 398, 557 (1983).
18H. Kuratsuji and T. Suzuki, “Path integral approach to many-nucleon systems and time-dependent Hartree-Fock,” Phys. Lett. B 92, 19 (1980).
19H. Kuratsuji and T. Suzuki, “Path integral approach to many-body systems and classical quantization of time-dependent mean field,” Suppl. Prog.
Theor. Phys. 74–75, 209 (1983).
20T. F. Viscondi and M. A. M. de Aguiar, “Semiclassical propagator for SU(n)
coherent states,” J. Math. Phys. 52(5), 052104 (2011).
21T. F. Viscondi and M. A. M. de Aguiar, “Initial value representation
for the SU(n) semiclassical propagator,” J. Chem. Phys. 134(23), 234105
(2011).
22T. Van Voorhis and E. J. Heller, “Nearly real trajectories in complex semiclassical dynamics,” Phys. Rev. A 66(5), 050501 (2002).
23T. Van Voorhis and E. J. Heller, “Similarity transformed semiclassical dynamics,” J. Chem. Phys. 119(23), 12153 (2003).
24We note that the approximations to the generalized coherent-state path
integral considered by Kuratsuji and Suzuki19—as well as specific
formulations for Slater determinants17,18—are very much akin to the
techniques developed in this paper.
25D. V. Shalashilin and M. S. Child, “Time dependent quantum propagation
in phase space,” J. Chem. Phys. 113(22), 10028 (2000).
26D. V. Shalashilin and M. S. Child, “Multidimensional quantum propagation
with the help of coupled coherent states,” J. Chem. Phys. 115(12), 5367
(2001).
27D. V. Shalashilin and M. S. Child, “The phase space CCS approach to quantum and semiclassical molecular dynamics for high-dimensional systems,”
Chem. Phys. 304, 103 (2004).
28The discussion is restricted to coherent states given in terms of analytical
complex parametrizations.
29Entries of the vector z will be referenced by Greek letters α, β, γ.
30T. F. Viscondi, Ph.D. thesis, Universidade Estadual de Campinas, Instituto
de Física Gleb Wataghin, Campinas, São Paulo, Brazil, 2013. Available in Portuguese, see http://www.bibliotecadigital.unicamp.br/document/
?code=000904028&opt=4&lg=en_US.
31M. Baranger, M. A. M. de Aguiar, F. Keck, H. J. Korsch, and B. Schellhaaß,
“Semiclassical approximations in phase space with coherent states,” J. Phys.
A: Math. Gen. 34, 7227 (2001).
32M. A. M. de Aguiar, S. A. Vitiello, and A. Grigolo, “An initial value
representation for the coherent state propagator with complex trajectories,”
Chem. Phys. 370, 42 (2010).
33In the majority of applications involving coherent states and classical
trajectories, the imaginary surface terms in A are canceled out by
normalization factors. However, when duplicated phase-space trajectories
are considered—i.e., when z and z ∗ are analytically continued in such a
way that they become completely independent variables (thereby doubling
the number of degrees of freedom in the system)—the complex action A
plays a prominent role, for its analytical properties are of the uttermost
importance in those cases.11,20,21,31,32
34E. J. Heller, “Frozen Gaussians: A very simple semiclassical approximation,” J. Chem. Phys. 75(6), 2923 (1981).
35E. J. Heller, “Cellular dynamics: A new semiclassical approach to timedependent quantum mechanics,” J. Chem. Phys. 94(4), 2723 (1991).
36M. F. Herman and E. Kluk, “A semiclassical justification for the use of
non-spreading wavepackets in dynamics calculations,” Chem. Phys. 91, 27
(1984).
of semiclassical frozen Gaussian wave functions with quantum propagation
for a highly excited anharmonic oscillator,” J. Chem. Phys. 84(1), 326
(1986).
38M. F. Herman, “Dynamics by semiclassical methods,” Annu. Rev. Phys.
Chem. 45, 83 (1994).
39K. G. Kay, “Integral expressions for the semiclassical time-dependent propagator,” J. Chem. Phys. 100(6), 4377 (1994).
40K. G. Kay, “The Herman-Kluk approximation: Derivation and semiclassical
corrections,” Chem. Phys. 322, 3 (2006).
41F. T. Arecchi, E. Courtens, R. Gilmore, and H. Thomas, “Atomic coherent
states in quantum optics,” Phys. Rev. A 6(6), 2211 (1972).
42Owing to the finite width of the Gaussian functions employed as basis, CCS
trajectories actually evolve under a smoothed-out classical Hamiltonian
whose width-dependent terms are usually referred to as “simple quantum
corrections.”
43SU(n) fermionic coherent states should be contrasted with the antisymmetrized Gaussian wavepackets mentioned in Sec. I. While the former are
defined in terms of a finite set of underlying molecular spin-orbitals, thus
providing a convenient description of the electronic structures of chemical
systems, the latter are more adequate for representing otherwise localized
fermions, such as protons and neutrons in the atomic nucleus5 or perhaps
outer shell electrons in highly excited atoms.3
44I. Burghardt, H.-D. Meyer, and L. S. Cederbaum, “Approaches to the approximate treatment of complex molecular systems by the multiconfiguration
time-dependent Hartree method,” J. Chem. Phys. 111(7), 2927 (1999).
45G. A. Worth and I. Burghardt, “Full quantum mechanical molecular dynamics using Gaussian wavepackets,” Chem. Phys. Lett. 368, 502 (2003).
46D. V. Shalashilin and I. Burghardt, “Gaussian-based techniques for quantum
propagation from the time-dependent variational principle: Formulation in
terms of trajectories of coupled classical and quantum variables,” J. Chem.
Phys. 129(8), 084104 (2008).
47For coherent states other than canonical, |⟨z|z ′⟩| is no longer a Gaussian
distribution, but is still localized in phase space.
48Latin letters j, k, l will be used for labeling basis-set elements and, for
convenience, we henceforth abbreviate basis-set summations by omitting
their ranges.
49K. G. Kay, “The matrix singularity problem in the time-dependent variational method,” Chem. Phys. 137, 165 (1989).
50S. Habershon, “Linear dependence and energy conservation in Gaussian
wavepacket basis sets,” J. Chem. Phys. 136(1), 014109 (2012).
51G. J. Milburn, J. Corney, E. M. Wright, and D. F. Walls, “Quantum dynamics
of an atomic Bose-Einstein condensate in a double-well potential,” Phys.
Rev. A 55(6), 4318 (1997).
52T. F. Viscondi, K. Furuya, and M. C. de Oliveira, “Generalized purity and
quantum phase transition for Bose-Einstein condensates in a symmetric
double well,” Phys. Rev. A 80, 013610 (2009).
53T. F. Viscondi, K. Furuya, and M. C. de Oliveira, “Coherent state approach
to the cross-collisional effects in the population dynamics of a two-mode
Bose-Einstein condensate,” Ann. Phys. 324(9), 1837 (2009).
54K. Sakmann, A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, “Exact
quantum dynamics of a bosonic Josephson junction,” Phys. Rev. Lett. 103,
220601 (2009).
55J. Schwinger, “On angular momentum,” in Quantum Theory of Angular
Momentum: A Collection of Reprints and Original Papers, edited by L. C.
Biedenharn and H. Van Dam (Academic Press, London, 1965), p. 229.
56Throughout this section units are such that ~ = 1.
57To compute Hamiltonian Eq. (4.2), one could either use decomposition
Eq. (2.23) to evaluate the required matrix elements of angular momentum
operators or employ SU(2) generating functions.9
58This can be verified by calculating the standard deviations of an appropriate
quantum phase-space distribution with respect to the angular coordinates
θ and φ.
59L. G. Yaffe, “Large N limits as classical mechanics,” Rev. Mod. Phys. 54(2),
407 (1982).
60T. F. Viscondi, K. Furuya, and M. C. de Oliveira, “Phase transition,
entanglement and squeezing in a triple-well condensate,” EPL 90, 10014
(2010).
61T. F. Viscondi and K. Furuya, “Dynamics of a Bose-Einstein condensate
in a symmetric triple-well trap,” J. Phys. A: Math. Theor. 44, 175301
(2011).
62We note that in the triple-well model, the energy difference between
the ground state and doubly degenerate excited states of the onebody Hamiltonian is 3~Ω—within the three-mode approximation, these
094106-15
Grigolo, Viscondi, and de Aguiar
stationary states span the same eigenspace as the local modes associated
with operators a 1, a 2, and a 3.61 Also, in Eq. (4.9), we have made an
additional simplification by excluding cross-collision terms, which arise
from the interaction between bosons in different wells.
63This means that the smaller basis sets are embedded in the larger ones.
64No significant improvement of the results was observed for larger basis
sets constructed with the same sampling parameters.
65As can be seen from Eq. (2.18), for the particular case of N = 1,
 n−1
the bosonic SU(n) coherent state reduces to |z} = |φ n ⟩ + α=1
z α |φ α ⟩,
where the occupation number representation has been replaced by “firstquantized” notation, with |φ α ⟩ being the single-particle orbital associated
with the α-th mode. This is nothing but a standard decomposition of a
quantum state in a finite basis of size n, only written without redundant
parameters.
J. Chem. Phys. 144, 094106 (2016)
66D. V. Shalashilin, “Nonadiabatic dynamics with the help of multiconfigura-
tional Ehrenfest method: Improved theory and fully quantum 24D simulation
of pyrazine,” J. Chem. Phys. 132(24), 244111 (2010).
67K. Saita and D. V. Shalashilin, “On-the-fly ab initio molecular dynamics with
multiconfigurational Ehrenfest method,” J. Chem. Phys. 137(22), 22A506
(2012).
68M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman, M. Booth, and F.
Rossi, GNU Scientific Library: Reference Manual, 3rd ed. (Network Theory
Ltd., 2003), http://www.gnu.org/software/gsl/.
69The overlap matrix is Hermitian and positive-definite, meaning that
its eigenvalues are real and positive, though numerical diagonalization
may produce null or very small negative eigenvalues. Alternatively, one
could employ a singular value decomposition and carry on the sampling
procedure using the singular values rather than the eigenvalues.
Article
pubs.acs.org/JPCA
Numerical Implementation and Test of the Modiﬁed Variational
Multiconﬁgurational Gaussian Method for High-Dimensional
Quantum Dynamics
Miklos Ronto and Dmitrii V. Shalashilin*
Downloaded via UNIV OF LEEDS on February 24, 2023 at 11:30:33 (UTC).
See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.
School of Chemistry, University of Leeds, Leeds LS2 9JT, U.K.
ABSTRACT: In this paper, a new numerical implementation and a test of the
modiﬁed variational multiconﬁgurational Gaussian (vMCG) equations are
presented. In vMCG, the wave function is represented as a superposition of
trajectory guided Gaussian coherent states, and the time derivatives of the wave
function parameters are found from a system of linear equations, which in turn
follows from the variational principle applied simultaneously to all wave
function parameters. In the original formulation of vMCG, the corresponding
matrix was not well-behaved and needed regularization, which required matrix
inversion. The new implementation of the modiﬁed vMCG equations seems to
have improved the method, which now enables straightforward solution of the
linear system without matrix inversion, thus achieving greater eﬃciency,
stability and robustness. Here the new version of the vMCG approach is tested
against a number of benchmarks, which previously have been studied by splitoperator, multiconﬁgurational time-dependent Hartree (MCTDH) and multilayer MCTDH (ML-MCTDH) techniques. The
accuracy and eﬃciency of the new implementation of vMCG is directly compared with the method of coupled coherent states
(CCS), another technique that uses trajectory guided grids. More generally we demonstrate that trajectory guided Gaussian
based methods are capable of simulating quantum systems with tens or even hundreds of degrees of freedom previously
accessible only for MCTDH and ML-MCTDH.
I. INTRODUCTION
Exact analytical solvability of the time-dependent Schrödinger
equation for systems with large number of degrees of freedom
(DOF) is limited to a few simple models. There are two
problems that make multidimensional quantum mechanics
diﬃcult to deal with. First, determining the potential energy
surface is a complicated problem, which is also present in
classical molecular dynamics. Recently, substantial progress has
been made1,2 with various forms of PES parametrizations and
ﬁts. The second problem is the scaling of quantum mechanics
with the number of degrees of freedom (DOF): the number of
quantum states increases exponentially with the size of
quantum system. Impressive progress has been made in
calculations of quantum states for the systems comprised of
large numbers of coupled vibrational modes3 some of which
can be very “ﬂoppy”.3,4 In dynamical calculationswhen
evolution of the wave packet is not restricted to a certain
areathe “exponential curse” of quantum mechanics is perhaps
the most severe. If a static grid of l states is used for a singlemode, a problem of M DOF requires
N = lM
of methods both semiclassical and formally exact have been
developed. The list of most important techniques can be
divided into two categories. Semiclassical methods include but
are not limited to the early application of frozen Gaussians5 and
semiclassical Herman−Kluk propagator.6,7 For more details and
new developments of the semiclassical Gaussian-based
methods, see review.8 Fully quantum techniques that at least
in principle can be converged to a fully quantum result include
multiple spawning,9 Gaussian multiconﬁgurational time-dependent Hartree (G-MCTDH),10 coupled coherent states
(CCS),11−13 variational Gaussian approach14 and variational
multiconﬁgurational Gaussians (vMCG).15 All these techniques
use grids (or basis sets) of trajectory-guided frozen Gaussian
coherent states (CSs), which follow the wave function, thus
economizing the basis set size. Another advantage of such
techniques is that a randomly sampled basis can be used, which
is advantageous in high-dimensional problems because Monte
Carlo techniques scale with dimensionality much better than eq
1.1. Importance sampling, which is the crucial part of all the
above methods, allows the basis to be built only around the
dynamically important phase-space region. With such random
(1.1)
grid points, making the challenge of exponential scaling almost
insurmountable.
In the past few decades a family of concepts based on
trajectory guided Gaussians, which rely on locally deﬁned
adaptable basis sets, gained considerable importance. A number
© 2013 American Chemical Society
Special Issue: Joel M. Bowman Festschrift
Received: November 6, 2012
Revised: April 8, 2013
Published: April 15, 2013
6948
dx.doi.org/10.1021/jp310976d | J. Phys. Chem. A 2013, 117, 6948−6959
The Journal of Physical Chemistry A
Article
1D problem allows one to visualize complicated variational
trajectories. Then Henon−Heiles (HH) model was investigated, and the accuracy and eﬃciency of vMCG was
compared with those of CCS, using the MCTDH benchmark
obtained previously for 2D, 6D, 10D and 18D HH systems.22
To demonstrate the limitations of vMCG and its scaling with
the number of degrees of freedom, we also performed CCS
calculation for 1458D model previously studied in ref 23 by
ML-MCTDH, and show that vMCG calculation for a system of
such size is not feasible, but CCS yields a good result. We also
use this example to make some general remarks about the
accuracy and convergence of vMCG, CCS, MCTDH and MLMCTDH methods24 of high-dimensional quantum mechanics.
Throughout this paper, natural units were used with ℏ = 1,
and the CSs were set to have ω = 1 and m = 1.
sampling, the Gaussian based methods can potentially scale
quadratically as N ∝ M2, although in reality the scaling is often
worse than that. However, even with the ever increasing
computational power, the methods based on trajectory-guided
Gaussians sometimes suﬀer from two diﬃculties: scaling and
robustness.
VMCG10,15 is potentially one of the most eﬃcient CS-based
methods of high-dimensional quantum dynamics. In vMCG,
the wave function is represented as a superposition of N
trajectory-guided Gaussian wave packets
∑ an(t )|pn(t ), q n(t ) = ∑ an(t )|zn(t )⟩
|Ψ(t )⟩ =
n = 1, N
n = 1, N
(1.2)
The time evolution of the wave function can be determined
from the variational principle. It yields time derivatives of its
parameters, which can be obtained from a system of linear
equations
Dα̇ = b
II. THEORY
II.1. Coherent States. In recent years, CS-based
approaches gained considerable importance in molecular
dynamics. As they are minimum uncertainty states and
Gaussian-functions, they are compatible with both semiclassical
and quantum mechanical descriptions. A single-mode CS z is
an element of the phase-space Γ = {(q,p)}, where qi and pi are
canonically conjugated variables. Using phase space coordinates, a 1D CS z ∈  can be written as
(1.3)
where α is the vector of wave function parameters, which
includes both the amplitudes an and all positions zn of N CSs
(M)
(each one of them M-dimensional zn = z(1)
n , …, zn ) and D is
the matrix that can be derived, for example, using the elegant
formalism developed by Kramer and Saraceno16 (See ref 17 for
more details). Matrix D in eq 1.3, as well as similar matrices in
other variational approaches,18−21 is often numerically nearly
singular and has to be regularized. This may be done by matrix
inversion, during which the lowest eigenvalues of D are
increased by a small and somewhat arbitrary regularization
parameter so that the inverse matrix does not have extremely
large elements. As a result, the equation is eventually written as
α̇ = D−1b
z=
γ
i
q+
ℏ
2
1
p
2γ
(2.1)
γ
i
q−
ℏ
2
1
p
2γ
(2.2)
and
z* =
(1.4)
with γ = mω/ℏ, where z and z* are eigenstates of the
annihilation and creation operators:
which is of course equivalent to eq 1.3 from the formal point of
view. Numerically, however, matrix inversion is much more
expensive than the solution of system of linear eqs 1.3.
In a recent paper,17 the equations of vMCG were modiﬁed
and written in a form similar to the equations of the CCS
technique,12 another method based on trajectory guided
Gaussians. In CCS, special eﬀorts were made to minimize
coupling between the amplitudes by making the coupling
matrix small, smooth and sparse. CCS introduces a
preexponential factor to smooth out the rapid oscillations of
quantum amplitudes. It appears that this simple trick works also
for vMCG and makes matrix D suﬃciently well behaved. In the
tests of the vMCG equations17 presented here, we were able to
solve eq 1.3 without inverting the matrix. Since our version of
vMCG is closely connected to the CCS theory, we also
compare the accuracy and eﬃciency of the two techniques.
With the same basis size, vMCG is more accurate simply
because it employs more variational parameters, as the time
dependence of all the parameters in eq 1.2 are determined
variationally. For the same number of variational parameters,
however, the accuracy of the two techniques is close, and both
methods show similar levels of robustness and stability.
In section II we brieﬂy sketch our version of the vMCG
theory, which, in the original paper,17 was presented for the
one-dimensional (1D) case only. Here we present the theory in
multidimensional form, which is conceptually simple but
involves long algebra, given in the Appendix A1. Chapter III
provides details of the numerical tests. First, our implementation of vMCG was tested on simple Morse oscillator because a
a|̂ z⟩ = z|z⟩
(2.3)
⟨z|a†̂ = ⟨z|z*
(2.4)
A multidimensional CS is a product of m single-mode CSs:
m
|zi(t )⟩ = ∏ |zi(k)(t )⟩
k=1
(2.5)
and represents a point in M-dimensional phase space ΓM =
{(qm,pm); m = 1,2,…,M}. The overlap of two M-dimensional
CSs
Ωij = ⟨zi|zj⟩
⎛
zjz*j ⎞
z z*
⎟
= exp⎜⎜z*i zj − i i −
2
2 ⎟⎠
⎝
M
= ∏ ⟨zi(k)|z(j k)⟩
k=1
⎛
|z(j k)|2 ⎞
|zi(k)|2
(k) (k)
⎜
⎟
*
= ∏ exp⎜zi zj −
−
2
2 ⎟⎠
⎝
k=1
M
(2.6)
is a product of M 1D overlaps. Although continuum basis of
CSs is overcomplete in numerical calculations, we always deal
with a ﬁnite set of CSs that simply represent a nonorthogonal
6949
dx.doi.org/10.1021/jp310976d | J. Phys. Chem. A 2013, 117, 6948−6959
The Journal of Physical Chemistry A
Article
basis with the overlap matrix 2.6. In this ﬁnite basis of CS the
identity becomes
be reparametrized with amplitudes and CS phase space
positions as17
 = ∑ |zi⟩Ω−ij 1⟨zj|
Λ(α , α*, α̇ , α̇*) = Λ(a , z , a*, z*, a ̇, z ̇, a*̇ , z*̇ )
⎡
⎛
i
= ∑ Ωij⎢ai*aj̇ − ai̇ *aj + ai*aj⎜z*i zj̇ − z*̇i zj
⎝
⎣
2
(2.7)
ij
where Ω−1
ij are the elements of the inverse overlap matrix of CS.
ij
The identity 2.7 is a discretization of that of the continuous
manifold covering all z-space where  is the product of M 1D
integral identities
M
⎛1
 = ∏⎜
⎝π
k=1
(2.8)
2
with the single-mode integral measure d z = (dq dp)/2ℏ. For
more details about CS notations, refer to ref 12. With the CS
representation, any time-dependent wave function and its
conjugate can be written up in CS basis in the general form 1.2.
In coordinate representation, a CS is simply a multidimensional
Gaussian wave packet with the following ansatz:
∂Λ
d ∂Λ
−
=0
∂α
dt ∂α̇
Dα̇ =
(2.10)
(2.11)
The index “ord” simply reminds about the terms that originated
from commuting â†,â. The Hamilton operator of the system is
the sum over the dimensions of the individual single-mode
Hamiltonians and their coupling terms:
⟨zi |: H(a†̂ (k) , a(̂ k)):| zj⟩ = ΩijHord(z*i , zj)
2
(2.13)
where the Lagrangian is
Λ = ⟨Ψ(t )|i ∂t⃡ − H |Ψ(t )⟩
(2.18)
but for smooth preexponential factor dj, where Sj is the action
along the trajectory. As a result, matrix D in eq 2.17 has many
small elements and is therefore sparse and smooth.
In vMCG, at each time step, one has to ﬁnd the derivatives of
the parameters from the system of coupled N × (M + 1) linear
eqs 2.17, where N is the basis set size and M is the number of
degrees of freedom. CCS is another method that utilizes exactly
the same parametrization of the wave function 1.2 as vMCG.
The diﬀerence is that in CCS the trajectories zn(t) are
predetermined and calculated from essentially classical
equations of motion. Only N amplitudes dn(t) are found
from quantum variational principle. As a result, the system of
linear equations for the derivatives of dn(t) is much smaller and
much simpler than in vMCG. Matrix D used in vMCG has the
size of [N × (M + 1)] × [N × (M + 1)] and includes that of
CCS as a small N × N block. The elements of this small N × N
matrix are simply those of the overlap matrix multiplied by the
exponentials of the classical actions. CCS trajectories are driven
by a classical Hamiltonian with quantum corrections Hord (eq
2.12), which is simply the expectation value of the classical
Hamiltonian with the Gaussian CS. The mathematical structure
of the two methods has been compared in ref 17. In this paper,
(2.12)
t1
∫t Λ(α , α*, α̇ , α̇*) dt
(2.17)
aj = dj exp(iSj)
2.2. Dynamics. Time-dependent variational principle
(TDVP) provides a generic way to derive various forms of
the time-dependent Schrödinger equation. Several formulations
of TDVP exist;16,18−21 our description here is based on the
approach16 where TDVP is presented in a form similar to the
principle of least action in classical mechanics and deﬁning
equations of motion through Euler−Lagrange equations.
According to the principle of least action, the equations of
motion can be obtained from the extremum of the functional
σ=
∂/
∂α*
where D is the matrix with the elements Dij = ∂παi/∂α*j and
⟨/⟩ is the eﬀective Hamiltonian obtained from the Lagrangian
(eq 2.14) in the usual way as ⟨/⟩ = παα − Λ . The operator
⟨/⟩ should not be confused with the actual physical
Hamiltonian of the system. The Lagrangian Λ and the eﬀective
Hamiltonian ⟨/⟩ are simply a tool to formally work out the
quantum equations of motion. We have used Greek letters π, σ,
and Λ to denote general momentum, action and the Lagrangian
associated with it to distinguish them from p, S and L, the
momentum, action and Lagrangian of the actual trajectory of a
physical coordinate q. More details about the equations and
their derivation can be found in ref 17 and in Appendix A1.
Here we only point out that equations17 have been written not
for the oscillating amplitude
(2.9)
This is also applicable to the single-mode Hamilton operators:
⟨z |: H(a†̂ , a):
̂ | z′⟩ = ⟨z|z′⟩Hord(z*, z)
(2.16)
∂Λ
Matrix elements of an arbitrary operator can be found via its
̂ in which the powers of the
normal ordered form : 6(a†̂ , a):
creation operator precede those of the annihilation operator if
â†,â are replaced by corresponding z or its complex conjugate
and the result is multiplied by the overlap:
⟨z |: 6(a†̂ , a):
̂ | z′⟩ = ⟨z|z′⟩6ord(z*, z)
(2.15)
which, after introduction of a conjugated momentum πα = 2 ∂α̇
can be presented in the form of Hamilton’s equations16
⎛ γ ⎞m /4 ⎛ γ
⟨x i|zj⟩ = ⎜ ⎟ exp⎜ − (x i − q i)T (x i − q i)
⎝π ⎠
⎝ 2
i pjq j ⎞
i
⎟
+ pj(x i − q i) +
ℏ
2ℏ ⎠
⎞⎤
1 *
(z j zj̇ + z*̇ j zj − z*i zi̇ − z*̇i zi) + 2i Hord(z*i , zj)⎟⎥
⎠⎦
2
where α = {a,z} is the vector of the wave function parameters
that includes all N × M components of the M-dimensional
complex vectors zj=1,N describing the phase space positions of
all basis CSs and N their amplitudes aj=1,N. Quantum equations
of motion can now be written as standard Lagrange equations
for the wave function parameters
(k) 2 (k)⎞
⎟
∫ |z ⟩⟨z |d z ⎠
(k)
−
(2.14)
with the diﬀerential-operator ∂t⃡ acting separately on bras and
kets. Any state vector |Ψ⟩ can be rewritten in CS form 1.2
|Ψ(α(t))⟩ = |Ψ(a(t),z(t))⟩, and therefore the Lagrangian can
6950
dx.doi.org/10.1021/jp310976d | J. Phys. Chem. A 2013, 117, 6948−6959
The Journal of Physical Chemist…

Turn in your highest-quality paper
Get a qualified writer to help you with

“ Science Question ”

Get high-quality paper

Guarantee! All work is written by expert writers!

Still stressed from student homework?

Get quality assistance from academic writers!

Order now