The quantum beating and its numerical simulation
Abstract.
We examine the suppression of quantum beating in a one dimensional nonlinear double well potential, made up of two focusing nonlinear point interactions. The investigation of the Schrödinger dynamics is reduced to the study of a system of coupled nonlinear Volterra integral equations. For various values of the geometric and dynamical parameters of the model we give analytical and numerical results on the way states, which are initially confined in one well, evolve. We show that already for a nonlinearity exponent well below the critical value there is complete suppression of the typical beating behavior characterizing the linear quantum case.
Key words and phrases:
Keywords: nonlinear Schrödinger equation, weakly singular Volterra integral equations, numerical computation of highly oscillatory integrals, quantum beating effect.1. Introduction
Quantum beating may nowadays refer to many, often quite different, phenomena studied in various domains of quantum physics, ranging from quantum electrodynamics to particle physics, from solid state physics to molecular structure and dynamics.
A paradigmatic example in the latter field is the inversion in the ammonia molecule observed experimentally in 1935. The ammonia molecule is pyramidally shaped. Three hydrogen atoms form the base and the nitrogen atom is located in one of the two distinguishable states (enantiomers) on one side or the other with respect to the base (chirality). Experimentally it was tested that microwave radiation could induce a periodic transition from one state to the other (quantum beating). It was also observed that in several circumstances the pyramidal inversion was suppressed. In particular in an ammonia gas the transition frequency was recognized to decrease for increasing pressure and the beating process proved to be finally suppressed for pressures above a critical value.
A theoretical explanation of the quantum beating phenomenon was obtained modeling the nitrogen atom (better, the two “nonbonding“ electrons of nitrogen) as a quantum particle in a double well potential.
The double well potential is of ubiquitous use in theoretical physics. In our present context, its importance consists in the fact that, for particular values of the parameters, the ground state and the first excited state have very close energies, forming an almost single, degenerate, energy level. A superposition of these two states is shown to evolve concentrating periodically inside one well or the other, with a frequency proportional to the energy difference (see section 2.1 below).
According to the mathematical quantum theory of molecular structure developed in the second half of the last century (see [1, 2] and references therein; see also [3, 4] for studies of the pressure dependent transition mechanism) the effect of the ammonia molecule quantum environment can be modeled as a nonlinear perturbation term added to the double well potential. A detailed quantitative analysis of the physical mechanism giving rise to the nonlinear reaction of the environment, in the case of pyramidal molecules, can be found in [5].
Following this suggestion, Grecchi, Martinez and Sacchetti ([6, 7, 8, 9]) investigated the semiclassical limit of solutions to the nonlinear Schrödinger initial value problem obtained perturbing the double well potential with a nonlinear term breaking the rotational symmetry. They were able to prove that there is a critical value of the coupling constant, measuring the strength of the symmetry breaking nonlinear term, above which the beating period goes to infinity meaning that the beating phenomenon is suppressed.
In this paper we present a different model of a similar physical situation. We consider a hamiltonian with two concentrated nonlinear attractive point potentials and we investigate the corresponding Cauchy problem. The study of the evolution problem is reduced to the analysis
of a system of two Volterra integral equations whose solutions we examine via numerical computation. It is worth stressing that the nonlinear model we consider is governed by symmetric dynamical equations. Asymmetry will appear only as a consequence of the nonlinearity and of specific initial conditions. We will come back to this point in the conclusions.
In Sections 2.1 and 2.2 we recall the properties of the corresponding symmetric and asymmetric linear case in order to clarify the origin of the beating phenomenon and its destruction. Due to the great degree of solvability of point interaction hamiltonians, the characterization of the beating states as functions of the dynamical and geometrical parameters of the model, will be carried through at a high level of detail.
In Section 3 we investigate, via numerical studies, the evolution problem in the linear symmetric, linear asymmetric and nonlinear cases. We show that the asymmetry resulting from the nonlinearity causes beating suppression and a rapid localization in one of the wells as soon as the nonlinearity becomes relevant.
In the conclusions we compare our results with what is known in literature and we list some open problems and some possible extensions of our results.
2. The mathematical model  Concentrated nonlinearities
In order to introduce the double well potential model we shall investigate in the present paper, we first briefly recall the definition of point interaction hamiltonians in (see ([10]) for further details).
For two point scatterers placed in of strength , the formal hamiltonian reads
(2.1) 
where the (reduced) Planck constant has been taken equal to one and the particle mass equal to .
A rigorous definition of in dimension has been given in the early days of quantum mechanics, when such kind of hamiltonians were extensively used to investigate the dynamics of a quantum particle in various kinds of short range scatterer arrays. A complete characterization of point interaction hamiltonians in was only made available in the second half of last century (see [10] for details and for an exhaustive bibliography).
Restricted to the case of our interest, definition and main results in are shortly detailed below. Assume that the two points are placed symmetrically with respect to the origin and that . Then
(2.2) 
(2.3) 
are domain and action of a selfadjoint operator in which acts as the free laplacian on functions supported outside the two points . In (2.2) is the Green function for the free laplacian, given by
(2.4) 
and the matrix is defined as
(2.5) 
where the positive real number is chosen large enough to make the matrix invertible.
It is immediate to check that the derivative of has a jump in the origin, equal to . This in turn implies that every function in the domain satisfies the boundary conditions
(2.6) 
The dynamics generated by is then characterized as the free dynamics outside the two scatterers, satisfying at any time the boundary conditions (2.6).
Our aim is to investigate the behaviour of the solutions to the nonautonomous evolution problem
(2.7) 
where the time dependence of is nonlinearly determined by the values in of the solution itself.
An alternative way to examine the Cauchy problem (2.7) is to write down Duhamel’s formula corresponding to the formal Hamiltonian (2.1) with the coupling constants given in (2.7), then prove that the boundary conditions are satisfied at each time.
In detail, let be the integral kernel of the unitary group
Then from the ansatz
(2.8) 
one obtains for
(2.9) 
Explicitly
(2.10) 
It is easy to check that a function of the form (2.8) satisfies the nonlinear boundary conditions at all times (see
[11] for details). Following a standard use in higher dimensional cases, we will often employ in this paper the notation and refer to (2.10) as the “charge equations”.
The Cauchy problem (2.7) is then reduced to the computation of and the solutions of the system (2.10), corresponding to two coupled nonlinear Volterra integral equations. The whole wavefunction is then recovered via (2.8).
In the following we will show that the linear case is characterized by the presence of almost stationary states whose wave function evolves periodically between one well and the other (beating states). Along the lines traced by many authors in the past (see [1],[2] and [7]), we will then show that the nonlinearity destroys the beating phenomenon. The reduction in complexity we obtain, using linear and nonlinear point interactions, makes the investigation of the theoretical and computational aspects of the problem easier. In order to better understand how the beating effect occurs and the reasons why one expects suppression of beating by nonlinear perturbation, we develop in Sections 2.1 and
2.2 the symmetric and antisymmetric linear cases in some detail.
2.1. Linear point interactions  Symmetric double well
Let us consider the symmetric linear case, corresponding to and . We will show that the eigenstates relative to the lowest eigenvalues are explicitly computable for the hamiltonian we consider.
which implies that the integral kernel of the resolvent is
(2.11) 
As it is clear from (2.11), the resolvent of is a finite rank perturbation of the free laplacian resolvent operator. From the kernel representation of the resolvent, spectral and scattering properties of the operator are easily inquired in the case of interactions of equal strength (see [10], Theorem 2.1.3).
Only the second term appearing in the formula for the resolvent (2.11) can have polar singularities for those positive values of for which the matrix is not invertible . In particular, will be a negative eigenvalue of if and only if
(2.12) 
In the case of two point interactions of the same strength this condition reads
(2.13) 
For there are two solutions to the previous equation. The indices stand for “fundamental” resp. first “excited” state. The corresponding eigenfunctions are
(2.14) 
(2.15) 
where and are easily computable normalization factors.
In Fig. 1 we plotted the two eigenstates and , corresponding to the fundamental state (symmetric function) and the first excited state (antisymmetric function). Notice that the two eigenstates are relative to energies getting closer and closer as the value of increases (see the remark below). In the same limit the absolute values of the two eigenfunctions tend to coincide.
The stationary solutions corresponding to these eigenstates are given by
A superposition of the two stationary states,
evolves as
In particular the superposition
(2.16) 
is concentrated in the left well and will evolve in time as
(2.17) 
with a probability density given by
(2.18) 
It is clear that is an oscillating function with period concentrated successively on the left and on the right well, justifying the definition of (2.17) as a beating state.
The values assumed by the function in the centers of the two wells evolve as follows (see Figure 2)
(2.19) 
Remark 2.1.
Many authors analyzed the energy difference between the fundamental and the first exited state of a hamiltonian with double well potential in the semiclassical limit, roughly referred to as (see e.g. [5]). In the notes [12] a detailed computation of the energy difference for a point interaction hamiltonian with two attractive zero range potentials is performed, keeping all the standard dimensions of the physical constants. In terms of the dimensionless constant it is proved that in the limit
(2.20) 
The exponential decay of the energy difference when is then easily and rigorously obtained in the case of a zero range double well. Furthermore the result clarifies that the semiclassical limit is characterized by , which in our units reads .
2.2. Linear point interactions  Asymmetric double well
Let us now investigate the changes in the beating mechanism when the two zero range potentials have different strengths , . In this asymmetric case the equation permitting to compute the eigenvalues of the Hamiltonian reads:
(2.21) 
leading to
(2.22) 
All the relevant results we will need in the following are collected in the following lemma, concerning the resolution of this last equation.
Lemma 2.2.
Let , and let us define the ratio . Then one has:
Proof.
Defining equation (2.22) can be rewritten as
(2.24) 
The number of positive solutions to (2.24) depends on the values of the parameters and on the distance .
 a:

For both the right and the left side of (2.24) are convex functions of , taking the common value when . Furthermore, denoting by the left hand side, we have that for . We deduce that there are two solutions to (2.24) if and only if the derivative of for is larger than (see Figure 3).
It is easy to check that if at least one of the is positive there cannot be two positive solutions of (2.24).
 b:

(2.25)  c:
∎
The above results suggest that in the semiclassical limit with the fundamental eigenstate approaches the eigenstate of a single point interaction of strength placed in whereas the excited state approaches the eigenstate of a single point interaction of strength placed in .
In order to make this aspect clearer, we detail the steps needed to perform an exact computation of the eigenfunctions associated to the two eigenvalues.
The normalized eigenfunction relative to the lowest eigenvalue has the form (see again Theorem 2.1.3 in [10])
(2.27) 
where , and the coefficient are solutions of
(2.28) 
which, together with (2.22) gives
(2.29) 
The normalization condition finally gives for
(2.30) 
Under the assumptions we made on , there will be a second eigenvalue , whose corresponding normalized eigenfunction has the form
(2.31) 
where
(2.32) 
with
(2.33) 
The thus obtained functions (2.27) resp. (2.31) are the eigenfunctions corresponding to the eigenvalues . The initial condition we shall choose in the asymmetric case will be of the form
(2.34) 
where the exact timeevolution of this state is given by
(2.35) 
Let us also remark here that from (2.29) and (2.32) we deduce that both and become negligible in the semiclassical regime, if . Taking into account the normalization factors (2.30), (2.33) and part (c) of lemma (2.2) we finally obtain that the fundamental state tends to and the excited state to (see Figure 5).
As a consequence the product turns out to be small everywhere and any periodic cancellation in (2.18) becomes impossible. No beating phenomenon will occur in this cases, as will be shown by numerical computations.
It should be expected that the asymmetry due to the nonlinearity will produce a similar behavior on time scales depending on the initial condition and on the strength of the nonlinearity.
2.3. Nonlinear point interactions
A detailed analytical study of the nonlinear case (which is no longer explicitly solvable) can be found in ([13, 11]). The authors obtained general results about existence of solutions either local or global in time and proved existence of blow up solutions for . In this section we briefly review the results that will be relevant for our work.
In Section 3 we present the numerical simulation results for the evolution of a beating state, i.e., an initial state giving rise in the linear case to a beating motion of the particle, namely
(2.36) 
Our aim is to study how the nonlinearity influences the beating phenomenon. As we already mentioned we expect that even if the initial condition is almostsymmetric, the nonlinearity will have the effect of braking the symmetry.
Let us go back to the general problem (2.7) with initial conditions
(2.37) 
which we will investigate in the integral form (2.10).
From ([11], Theorem 6) we know that, if and one chooses an initial data , then the Cauchy problem has a unique solution which is global in time. Moreover in ([11], Theorem 23) it is proved that if and then there exist initial data such that the solutions of the Cauchy problem blowup in finite time.
For convenience of the reader we rewrite below the charge equation
3. The numerical discretization of the Volterrasystem
Let us come now to the numerical part of this work, namely the discretization and later on simulation of the Volterrasystem (2.38), in order to investigate the delicate phenomenon of beating. Linear (symmetric and asymmetric) as well as nonlinear cases will be treated, starting from an initial condition under one of the forms
(3.39) 
with some given constants and resp. defined in (2.14)(2.15) resp. (2.27)(2.31).
The discretization of the Volterrasystem (2.38) passes through the discretization of two different kind of integrals, an Abelintegral, which is of the form
(3.40) 
and a highlyoscillating integral of the form
(3.41) 
Besides, the free Schrödinger equation
(3.42) 
has to be solved to compute the right hand side of the Volterra system, i.e. , and one has also to take care of the nonlinearity, which will be treated iteratively by means of a linearization. The treatment of all these four steps shall be presented in the following subsections.
For the numerics we shall consider the truncated timespace domain and impose periodic boundary conditions in space. We shall furthermore fix a homogeneous discretization of this domain, defined as
3.1. The free Schrödinger evolution
We shall present now two different resolutions of the Schrödinger equation (3.42), a numerical resolution via the Fast Fourier Transform (fft,ifft) assuming periodic boundary conditions in space and an analytic, explicit resolution by means of the continuous Fourier Transform and based on the specific initial condition we choose.
The numerical resolution starts from the partial FourierTransform (in space) of (3.42)
where
and hence
(3.43) 
Remark that we supposed here periodic boundary conditions in the truncated spacedomain , where the appearance of the discrete Fouriervariable . Using the fft as well as ifftalgorithms permits hence to get from (3.43) a numerical approximation of the solution of the free Schrödinger equation (3.42).
Analytically, we shall situate us in the whole space and shall perform the same steps explicitly, taking advantage of the initial condition, which has the form
(3.44) 
where we recall that (see (2.4), (2.14), (2.15))
Thus one has with the definition of the Fouriertransform and its inverse
that
Let us now compute explicitly the Fourier transform of and and finally the inverse Fourier transform of . For this, remark that one has
leading to
Now, in the aim to resolve numerically the Volterrasystem (2.38), one needs only to compute the solution of (3.42) in the points , which means
To compute these two integrals, we shall take advantage of the following two formulae
where is the socalled errorfunction, defined by
After some straightforward computations, one gets
With the two expressions and one has now
which permits to have the righthand side of the Volterrasystem (2.38) analytically.
Let us observe that the same computations hold also for the asymmetric initial condition (2.34) with replaced by , as well as by .
3.2. The Abel integral
Let us now present a discretization of an Abelintegral of the form (3.40), based on a Gaussian quadrature. The time interval is discretized in a homogeneous manner, as proposed above, such that one can now approximate for as follows
Now, introducing the notation