The quantum beating

The quantum beating and its numerical simulation

Raffaele Carlone Università “Federico II” di Napoli, Dipartimento di Matematica e Applicazioni “R. Caccioppoli”, MSA, via Cinthia, I-80126, Napoli, Italy. Rodolfo Figari Università “Federico II” di Napoli, Dipartimento di Fisica e INFN Sezione di Napoli, MSA, I-80126, Napoli, via Cinthia, Italy.  and  Claudia Negulescu Université de Toulouse & CNRS, UPS, Institut de Mathématiques de Toulouse UMR 5219, F-31062 Toulouse, France

We examine the suppression of quantum beating in a one dimensional non-linear 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: non-linear 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 “non-bonding“ 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 non-linear perturbation term added to the double well potential. A detailed quantitative analysis of the physical mechanism giving rise to the non-linear 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 semi-classical limit of solutions to the non-linear Schrödinger initial value problem obtained perturbing the double well potential with a non-linear 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 non-linear 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 non-linear 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 non-linear model we consider is governed by symmetric dynamical equations. Asymmetry will appear only as a consequence of the non-linearity 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 non-linear cases. We show that the asymmetry resulting from the non-linearity causes beating suppression and a rapid localization in one of the wells as soon as the non-linearity 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


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


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


and the matrix is defined as


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


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 non-autonomous evolution problem


where the time dependence of is non-linearly 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


one obtains for




It is easy to check that a function of the form (2.8) satisfies the non-linear 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 wave-function 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 non-linear 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.

In fact, applying to both sides of (2.3) and using (2.2) we obtain that for all one has

which implies that the integral kernel of the resolvent is


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


In the case of two point interactions of the same strength this condition reads


For there are two solutions to the previous equation. The indices stand for “fundamental” resp. first “excited” state. The corresponding eigenfunctions are


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 (anti-symmetric 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.

Figure 1. Plot of the functions with a thicked blue line and with a dashed line.

The stationary solutions corresponding to these eigenstates are given by

A superposition of the two stationary states,

evolves as

In particular the superposition


is concentrated in the left well and will evolve in time as


with a probability density given by


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)

Figure 2. Plot of the time-evolution of the functions as a dashed line, as a dotted line and as a thick line.
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 semi-classical 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


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:


leading to


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:


There are two real solutions to equation (2.22) if and only if for and


For , , satisfying (2.23) and , one has

In particular as  .


For , , satisfying (2.23) and , one has


Defining equation (2.22) can be rewritten as


The number of positive solutions to (2.24) depends on the values of the parameters and on the distance .


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).

Figure 3. Plot of the functions in red, dashed for , and .

As it is clear from Figure 3, the two solutions to (2.24) are such that and . As a consequence


In the semi-classical regime , (see remark at the end of previous subsection), (2.25) implies that the energy difference becomes larger and larger whenever (see Figure 4). It is worth recalling that for (symmetric case) the energy difference goes to zero in the same limit.


Rewriting (2.24) in terms of and we obtain


Both solutions and are strictly larger than zero. In turn, this implies that in the semi-classical limit , the exponential term in (2.26) becomes negligible with respect to and , .

Figure 4. Contour plot of the functions solutions to the equation (2.22) as a function of and for , .

The above results suggest that in the semi-classical 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])


where , and the coefficient are solutions of


which, together with (2.22) gives


The normalization condition finally gives for


Under the assumptions we made on , there will be a second eigenvalue , whose corresponding normalized eigenfunction has the form






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


where the exact time-evolution of this state is given by


Let us also remark here that from (2.29) and (2.32) we deduce that both and become negligible in the semi-classical 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).

Figure 5. Plot of and for , , .

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 non-linearity 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 non-linear 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


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 almost-symmetric, the nonlinearity will have the effect of braking the symmetry.

Let us go back to the general problem (2.7) with initial conditions


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 blow-up in finite time.

For convenience of the reader we re-write below the charge equation


Next section is devoted to test the effectiveness of the integral form (2.38) of the evolution equations to find numerical solutions of the Cauchy problem (2.7).

3. The numerical discretization of the Volterra-system

Let us come now to the numerical part of this work, namely the discretization and later on simulation of the Volterra-system (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


with some given constants and resp. defined in (2.14)-(2.15) resp. (2.27)-(2.31).

The discretization of the Volterra-system (2.38) passes through the discretization of two different kind of integrals, an Abel-integral, which is of the form


and a highly-oscillating integral of the form


Besides, the free Schrödinger equation


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 non-linearity, 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 time-space 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 Fourier-Transform (in space) of (3.42)


and hence


Remark that we supposed here periodic boundary conditions in the truncated space-domain , where the appearance of the discrete Fourier-variable . Using the fft- as well as ifft-algorithms 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


where we recall that (see (2.4), (2.14), (2.15))

Thus one has with the definition of the Fourier-transform and its inverse


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 Volterra-system (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 so-called error-function, defined by

After some straightforward computations, one gets

With the two expressions and one has now

which permits to have the right-hand side of the Volterra-system (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 Abel-integral 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