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
###### Abstract.

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

 Hγ––,Yψ:=‘‘−d2dx2ψ+γ1δy1ψ+γ2δy2ψ‘‘, (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)
 (Hγ––,Y+λ)ψ=(−d2dx2+λ)ϕλ, (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

 Gλ(x):=e−√λ|x|2√λ, (2.4)

and the matrix is defined as

 (Γλγ––)ij :=1γiδij+Gλ(yi−yj), (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

 dψdx(y+j)−dψdx(y−j)=γjψ(yj),j=1,2. (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 non-autonomous evolution problem

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩ı∂ψ∂t=Hγ––(t),Yψ,∀(t,x)∈R+×R,ψ(0,x)=ψ0(x)∈D(Hγ––(0),Y)∀x∈R,γj(t):=γ|ψ(t,yj)|2σ,γ<0,σ≥0. (2.7)

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

 U(τ,y):=eı|y|24τ√4ıπτ,(U(t)ξ)(x)=∫∞−∞U(t;x−y)ξ(y)dy∀ξ∈L2(R).

Then from the ansatz

 ψ(t,x)=(U(t)ψ0)(x)−ıγ2∑j=1∫t0U(t−s;x−yj)|ψ(s,yj)|2σψ(s,yj)ds, (2.8)

one obtains for

 ψ(t,yi)=(U(t)ψ0)(yi)−ıγ2∑j=1∫t0U(t−s;yi−yj)|ψ(s,yj)|2σψ(s,yj)ds. (2.9)

Explicitly

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ψ(t,−a)+γ2√ıπ∫t0ψ(s,−a)|ψ(s,−a)|2σ√t−sds+γ2√ıπ∫t0ψ(s,a)|ψ(s,a)|2σ√t−seıa2(t−s)ds=(U(t)ψ0)(−a),ψ(t,a)+γ2√ıπ∫t0ψ(s,a)|ψ(s,a)|2σ√t−sds+γ2√ıπ∫t0ψ(s,−a)|ψ(s,−a)|2σ√t−seıa2(t−s)ds=(U(t)ψ0)(a). (2.10)

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

 (Hγ––,Y+λ)−1ξ=(−d2dx2+λ)−1ξ−2∑i,j=1(Γλγ––)−1ij[(−d2dx2+λ)−1ξ](yj)Gλ(⋅−yi)

which implies that the integral kernel of the resolvent is

 (Hγ––,Y+λ)−1(x,x′)=Gλ(x−x′)−2∑i,j=1(Γλγ––)−1ijGλ(x−yi)Gλ(x′−yj). (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

 detΓλγ––=0. (2.12)

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

 det⎛⎜⎝1γ+12√λGλ(2a)Gλ(2a)1γ+12√λ⎞⎟⎠=0. (2.13)

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

 ϕf(x)=Nf(Gλf(x+a)+Gλf(x−a)) (2.14)
 ϕe(x)=Ne(Gλe(x+a)−Gλe(x−a)), (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 (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.

The stationary solutions corresponding to these eigenstates are given by

 ψf(t,x)=eıλftϕf(x),ψe(t,x)=eıλetϕe(x).

A superposition of the two stationary states,

 ψ0(x)=αϕf(x)+βϕe(x),α,β∈R,

evolves as

 ψ(t,x)=αeıλftϕf(x)+βeıλetϕe(x).

In particular the superposition

 ψLbeat,0(x):=1√2(ϕf(x)+ϕe(x)), (2.16)

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

 ψLbeat(t,x)=1√2(eıλftϕf(x)+eıλetϕe(x)), (2.17)

with a probability density given by

 P(t,x)=12[|ϕf(x)|2+|ϕe(x)|2+2ϕf(x)ϕe(x)cos((λf−λe)t)]. (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)

 qL1(t)≡ψLbeat(t,−a)=1√2(eıλftϕf(−a)+eıλetϕe(−a))qL2(t)≡ψLbeat(t,a)=1√2(eıλftϕf(a)+eıλetϕe(a)). (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 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

 △E≃2mγ2ℏ2e−δ. (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:

 detΓλ(γ1,γ2)=det⎛⎜⎝1γ1+12√λGλ(2a)Gλ(2a)1γ2+12√λ⎞⎟⎠=0, (2.21)

 (1γ1+12√λ)(1γ2+12√λ)−(12√λ)2e−4√λa=0. (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:

a:

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

 1|γ1|+1|γ2|<2a. (2.23)
b:

For , , satisfying (2.23) and , one has

 Δλ:=λ0−λ1≥γ21(1−α2).

In particular as  .

c:

For , , satisfying (2.23) and , one has

 lim|γ1|→∞2√λ0/γ1=−1,lim|γ1|→∞2√λ1/γ2=−1.
###### Proof.

Defining equation (2.22) can be rewritten as

 ξ2γ1γ2+ξ(1γ1+1γ2)+1=e−2ξa. (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:

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

 4Δλ=ξ20−ξ21≥γ21−γ22=γ21(1−α2). (2.25)

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.

c:

Rewriting (2.24) in terms of and we obtain

 η2α−(1+αα)η+1=e−2η|γ1|a. (2.26)

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

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

 ϕ0(x)=c0Gλ0(x−y1)+c1Gλ0(x−y2), (2.27)

where , and the coefficient are solutions of

 ⎛⎜ ⎜⎝1γ1+12√λ012√λ0e−2√λ0a12√λ0e−2√λ0a1γ2+12√λ0⎞⎟ ⎟⎠(c0c1)=(00), (2.28)

which, together with (2.22) gives

 ∣∣∣c1c0∣∣∣= ⎷(2√λ0/γ1)+1(2√λ0/γ2)+1. (2.29)

The normalization condition finally gives for

 c0=2|γ1|λ3/40 ⎷γ1γ2(γ1+2√λ0)(γ2+2√λ0)+γ1(γ1+4√λ0+2√λ0a(γ1+2√λ0)). (2.30)

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

 ϕ1(x)=c2Gλ1(x−y1)+c3Gλ1(x−y2), (2.31)

where

 ∣∣∣c2c3∣∣∣= ⎷(2√λ1/γ2)+1(2√λ1/γ1)+1, (2.32)

with

 c3=2|γ2|λ3/41 ⎷γ1γ2(γ2+2√λ1)(γ1+2√λ1)−γ2(γ2+4√λ1+2√λ1a(γ2+2√λ1)). (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

 ψasy,0(x):=αϕ0(x)+βϕ1(x),α,β∈R, (2.34)

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

 ψasy(t,x):=αeiλ0tϕ0(x)+βeiλ1tϕ1(x). (2.35)

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

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

 ψ0(x):=αϕf(x)+βϕe(x),α,β∈R. (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 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

 ψ(0,x)=ψ0(x). (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 blow-up in finite time.

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

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩q1(t)+γ2√ıπ∫t0q1(s)|q1(s)|2σ√t−sds+γ2√ıπ∫t0q2(s)|q2(s)|2σ√t−seıa2t−sds=(U(t)ψ0)(−a),q2(t)+γ2√ıπ∫t0q2(s)|q2(s)|2σ√t−sds+γ2√ıπ∫t0q1(s)|q1(s)|2σ√t−seıa2t−sds=(U(t)ψ0)(a). (2.38)

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

 ψLbeat,0(x):=αϕf(x)+βϕe(x),ψasy,0(x):=αϕ0(x)+βϕ1(x), (3.39)

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

 Ab(t):=∫t0g(s)√t−sds, (3.40)

and a highly-oscillating integral of the form

 Ho(t):=∫t0g(s)√t−seia2t−sds. (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 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

 0=t1<⋯

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

 ⎧⎨⎩ddt^θk(t)=−ik2^θk(t),∀k∈Z,∀t∈R+,^θk(0)=^ψ0,k,∀k∈Z,

where

and hence

 ^θk(t)=e−ik2t^θk(0),∀(t,k)∈R+×Z. (3.43)

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

 (3.44)

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

 ϕf(x)=Nf[Gλf(x+a)+Gλf(x−a)],ϕe(x)=Ne[Gλe(x+a)−Gλe(x−a)].

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

 ^ϕ(ν):=1√2π∫∞−∞ϕ(x)e−ixνdx,ϕ(x)=1√2π∫∞−∞^ϕ(ν)eixνdν,

that

 ^ψ0(ν)=α^ϕf(ν)+β^ϕe(ν)⇒^ϑ(t,ν)=α^ϕf(ν)e−iν2t+β^ϕe(ν)e−iν2t(t,ν)∈R+×R.

Let us now compute explicitly the Fourier transform of and and finally the inverse Fourier transform of . For this, remark that one has

 ^Gλ(ν)=1√2π1λ+ν2,∀ν∈R,

 ^ϕf(ν)=2Nf√2πcos(νa)λf+ν2,^ϕe(ν)=−2iNe√2πsin(νa)λe+ν2.

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

 ϑ(t,−a)=αNf2π∫∞−∞1+e−2iaνλf+ν2e−iν2tdν+βNe2π∫∞−∞1−e−2iaνλe+ν2e−iν2tdν.
 ϑ(t,a)=αNf2π∫∞−∞1+e2iaνλf+ν2e−iν2tdν−βNe2π∫∞−∞1−e2iaνλe+ν2e−iν2tdν,

To compute these two integrals, we shall take advantage of the following two formulae

 IλB:=∫∞−∞e±2iaνλ+ν2e−iν2tdν=∫∞−∞cos(2aν)λ+ν2e−iν2tdν,

where is the so-called error-function, defined by

 erf(x):=2√π∫x0e−t2dt.

After some straightforward computations, one gets

 IλB=π2√λeiλt{e2√λa[1−erf(√iλt+a√it)]+e−2√λa[1−erf(√iλt−a√it)]}.

With the two expressions and one has now

 (U(t)ψ0)(−a)=ϑ(t,−a)=αNf2π[IλfA+IλfB]−βNe2π[IλeA+IλeB],
 (U(t)ψ0)(a)=ϑ(t,a)=αNf2π[IλfA+IλfB]+βNe2π[IλeA+IλeB],

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

 Ab(tl)=∫tl0g(s)√tl−sds=l−1∑k=1∫tk+1tkg(s)√tl−sds=l−1∑k=1√Δt∫10g(tk+ξΔt)√l−k−ξdξ.

Now, introducing the notation

 r(l)k:=l−k,φk(ξ):