text:
# Discrete matrix Riccati equations with superposition formulas

## 1 Introduction

The well-known Riccati equation

$$\dot{w}=a(t)+b(t)w+c(t)w^2$$
(1)
arises in numerous contexts. For example, one can find particular cases of (1) as equations describing Darboux transformation for the Schrödinger operator , as equations describing particular solutions of Painlevé equations , or as Bäcklund transformations for soliton like equations .

The Riccati equation has many interesting properties, one of them is the fact that it is a non-linear equation with a superposition formula. A differential equation has a superposition formula if it is possible to express its general solution as a function of a finite number of particular solutions and arbitrary constants.

It is easy to prove that the Riccati equation (1) has a superposition formula. Let $`w_0,w_1,w_2`$ be three different particular solutions of Eq. (1). It is well-known that using one particular solution one can transform (1) into a linear equation. Indeed, let $`w=y+w_0,`$ then we have $`\dot{y}=(b+2cw_0)y+cy^2.`$ Now let $`y=\frac{1}{x},`$ then we obtain a linear equation $`\dot{x}=(b+2cw_0)xc.`$ We know two particular solutions $`x_1=\frac{1}{w_1w_0}`$ and $`x_2=\frac{1}{w_2w_0}`$ of this linear equation, hence we can find a general solution of this linear equation using $`x_1`$, $`x_2`$ and one arbitrary constant. This gives us a formula for the general solution of the Riccati equation in terms of three particular solutions $`w_0,w_1,w_2`$ and one arbitrary constant $`c:`$

$$w=w_0+\frac{(w_2w_0)(w_1w_0)}{w_2w_0+c(w_1w_2)}.$$
(2)
One can also consider the matrix Riccati equation

$$\dot{W}=A(t)+B(t)W+WC(t)+WD(t)W,$$
(3)
where $`W`$ is a $`n\times k`$-matrix and $`A,B,C,D`$ are matrices of appropriate sizes. This equation also has a superposition formula .

One can discretize the Riccati equations (1) and (3) in different ways, but we are interested in discretizations possessing superposition formulas. Our main result is the following. Let us consider a two-parameter class of operators $`U_{q,h}:`$

$$U_{q,h}f(t)=\{\begin{array}{cc}f^{}(t)\hfill & \text{if}q=1,h=0\hfill \\ \frac{f(qt+h)f(t)}{(q1)t+h}\hfill & \text{in other cases.}\hfill \end{array}$$
(4)
As we shall show below the equation

$$U_{q,h}w(t)=a(t)+b(t)w(t)+w(qt+h)c(t)+w(qt+h)d(t)w(t)$$
(5)
has a superposition formula. In the case where $`q=1,h=0`$ we obtain the Riccati equation (1). It can be easily seen that in the case $`q=1,h0`$ we obtain the $`h`$-Riccati equation

$$\frac{w(t+h)w(t)}{h}=a(t)+b(t)w(t)+w(t+h)c(t)+w(t+h)d(t)w(t),$$

and in the case $`h=0,q1`$ we obtain the $`q`$-Riccati equation

$$\frac{w(qt)w(t)}{(q1)t}=a(t)+b(t)w(t)+w(qt)c(t)+w(qt)d(t)w(t).$$

Analogous results hold also in the matrix case.

It is interesting to remark that particular cases of such discretizations appear in the same contexts as their classical analogue, for example as Darboux transformations for the discrete Schrödinger operator , or as equations describing particular solutions of discrete Painlevé equations . Thus the discretizations constructed in the present paper are very natural.

## 2 Differential and difference equations with superposition formulas

Let us consider a (vector) first order ordinary differential equation

$$\dot{y}^k=f^k(𝐲,t),k=1,\mathrm{},n.$$
(6)
We shall say that this equation has a superposition formula if its general solution $`𝐲(t)`$ can be expressed as a function of a finite number $`m`$ of particular solutions $`𝐲_1,\mathrm{},𝐲_m`$ and $`n`$ free constants

$$𝐲(t)=𝐅(𝐲_1,\mathrm{},𝐲_m,c_1,\mathrm{},c_n).$$
(7)
The formula (7) is called a superposition formula for Eq. (6).

The study of ordinary differential equations with superposition formulas has a long history that goes back to Sophus Lie. He proved the following fundamental theorem.

###### Theorem 1 (Lie )

The equation (6) has a superposition formula if and only if the function $`𝐟`$ has the form

$$f^k(𝐲,t)=\underset{l=1}{\overset{r}{}}Z_l(t)\xi _l^k(𝐲),$$

where the functions $`\xi _l(𝐲)`$ are such that the vector fields

$$X_l=\underset{k=1}{\overset{n}{}}\xi _l^k(𝐲)_{y^k}$$

generate a finite dimensional Lie subalgebra $`𝔥`$ of the algebra of vector fields of $`^N`$ or $`^N,`$ i.e.

$$[X_i,X_j]=\underset{k=1}{\overset{r}{}}C_{ijk}X_k.$$

Difference equations with superposition formulas can be defined in a similar way. In order to see what can be done in the difference case, let us reformulate the Lie theorem in different terms.

Let $`G`$ be a Lie group acting on a manifold $`M.`$ Let us consider a curve $`g(t)`$ on $`G`$ such that $`g(t_0)=e`$ and a point $`u_0M.`$ It is easy to find a differential equation for a curve $`u(t)=g(t)u_0`$ on $`M.`$ Indeed,

$$\dot{u}(t)=\dot{g}(t)u_0=\dot{g}(t)g(t)^1g(t)u_0=\dot{g}(t)g(t)^1u(t).$$

Since $`\dot{g}g^1𝔤,`$ we see that the following statement holds.

Statement 1. A curve $`u(t)=g(t)u_0,`$ such that $`g(t_0)=e,`$ is a solution of a differential equation

$$\dot{u}(t)=\xi (t)u(t),$$
(8)
where $`\xi =\dot{g}g^1𝔤.`$ The point $`u_0`$ plays the role of an initial condition, $`u(t_0)=u_0.`$

Let us now consider Eq. (8) with given $`\xi (t).`$ Let $`u(t)`$ be a solution. Let us choose a fixed value of the parameter $`t_0`$ and the initial value of the solution $`u_0=u(t_0).`$ There exists a unique curve $`g(t)G`$ such that $`\dot{g}g^1=\xi `$ and $`g(t_0)=e.`$ We see that the curve $`\stackrel{~}{u}(t)=g(t)u_0`$ is also a solution of Eq. (8), and it also satisfies the initial condition $`\stackrel{~}{u}(t_0)=u_0.`$ The following statement follows from the uniqueness theorem for solutions of ODEs.

Statement 2. Let $`\xi (t)`$ be a curve on $`𝔤.`$ Let $`t_0`$ be a fixed value of the parameter. Any solution $`u(t)`$ of equation (8) has the form $`u(t)=g(t)u_0,`$ where $`u_0M`$ is the initial condition $`u(t_0)=u_0`$ and $`g(t)`$ is a curve on $`G`$ such that $`g(t_0)=e.`$ The curve $`g(t)G`$ is uniquely determined by the conditions $`\dot{g}g^1=\xi `$ and $`g(t_0)=e,`$ hence $`g(t)`$ is the same for all solutions $`u(t)`$ of equation (8).

Statements 1 and 2 mean that Eq. (8) has a superposition formula. Indeed, let $`u_1(t),\mathrm{},u_m(t)`$ be particular solutions of (8). Let us choose an initial value of the parameter $`t_0.`$ It follows from Statement 2 that $`u_1(t),\mathrm{},u_m(t)`$ have the form $`u_i(t)=g(t)u_i(t_0),`$ where $`g(t)`$ is the same for all $`i.`$ Then we have the following system of equations for $`g(t):`$

$$\{\begin{array}{ccc}\hfill u_1(t)& =& g(t)u_1(t_0),\hfill \\ & \mathrm{}& \\ \hfill u_m(t)& =& g(t)u_m(t_0).\hfill \end{array}$$

If $`m`$ is sufficiently large, $`g(t)`$ can be expressed in terms of $`u_i(t)`$ and $`u_i(t_0):`$

$$g(t)=F(u_1(t),\mathrm{},u_m(t);u_1(t_0),\mathrm{},u_m(t_0)).$$

It was shown by Lie that $`m`$ should at least satisfy the inequality $`mnr,`$ where $`n=dimM`$ and $`r`$ is the same $`r`$ as in Theorem 1, i. e. the dimension of the image of $`𝔤`$ in the algebra of vector fields on $`M.`$

It follows that the general solution has the form

$$u(t)=g(t)u_0=F(u_1(t),\mathrm{},u_m(t);u_1(t_0),\mathrm{},u_m(t_0))u(t_0),$$
(9)
the initial condition $`u(t_0)`$ plays the role of the arbitrary constant in the superposition formula (9). It should be remarked that the function $`F`$ in the superposition formula (9) does not depend on the choice of $`t_0,`$ hence the superposition formula is essentially the same for any choice of $`t_0.`$ The superposition formula obtained by choosing another initial value of the parameter $`t_1`$ can be obtained form (9) by expressing the constants $`u_1(t_0),\mathrm{},u_m(t_0),u(t_0)`$ in terms of $`u_1(t_1),\mathrm{},u_m(t_1),u(t_1).`$

It should also be mentioned that generally we also have some independence conditions on the particular solutions in order to find $`g(t).`$ For example in the case of a homogeneous linear equation of second order any solution is a linear combination of two particular solutions under the condition that these two solutions are linearly independent.

Statements 1 and 2 are not complicated. However, the Lie theorem also contains the following statement which is more difficult to prove.

Statement 3. Every ordinary differential equation with a superposition formula must have the form (8) for some group $`G`$ acting on some manifold $`M`$ .

A modern exposition of the proof can be found in . Statements 1, 2 and 3 are exactly what Lie’s theorem states. There is a very important thing used implicitly. In the proof of Statement 2 a crucial step makes use of the uniqueness theorem for solutions of ODEs. We are working implicitly with some class of functions such that for solutions from this class this theorem holds. In the same implicit manner we assume that equation (8) is an equation for which the uniqueness theorem holds. For differential equations the assumptions underlying the uniqueness theorem are very natural. One can assume that Eq. (8) is defined by a continuous curve $`\xi (t)`$ and that we are considering solutions in the class $`C^1(),`$ i. e. the class of differentiable functions with continuous derivative. For difference equations the situation is different and the appropriate class of solutions must be carefully chosen and we need a more rigorous definition.

Definition. A differential (difference) equation is said to have a superposition formula in some class of solutions if the general solution belonging to this class can be expressed as a function of a finite number $`m`$ of particular solutions $`𝐲_1,\mathrm{},𝐲_m,`$ also belonging to this class, and free constants

$$𝐲=𝐅(𝐲_1,\mathrm{},𝐲_m,c_1,\mathrm{},c_n),$$
(10)
where $`n=dimM.`$ Eq. (10) is called a superposition formula.

The question of describing all difference equations with superposition principles is beyond the scope of this paper. It is a very difficult question, especially since different difference operators require the consideration of completely different classes of solutions. We will see examples later in the case of the discrete Riccati equation. Let us concentrate on Statements 1 and 2 and see how one can construct difference equations with superposition formulas.

We can start again from a Lie group $`G,`$ a manifold $`M,`$ a curve $`g(t)M`$ such that $`g(t_0)=e`$ and a point $`u_0M.`$ Let $`u(t)=g(t)u_0.`$ Let $`U`$ be a first order difference operator, for example of the form (4) defined in the Introduction. The difference derivative $`Uu(t)`$ can be defined only in a fixed coordinate system on $`M.`$ It is no longer an invariant object like the vector $`\dot{u}(t).`$ It is impossible, in general, to write something like $`Uu(t)=Ug(t)u_0`$ and reduce all to the Lie algebra, since there is no general chain rule for a difference operator and $`Ug(t)`$ is not, in general, an element of $`𝔤.`$ However, in some cases we can construct surrogates of statements 1 and 2 in the following way. Let us assume that for some $`G`$ and $`M`$ (with a fixed coordinate chart) the following statements are true.

Statement 1’. A curve $`u(t)=g(t)u_0`$ such that $`g(t_0)=e`$ is a solution of some difference equation

$$Uu(t)=H(u(t)).$$
(11)
Statement 2’. Let $`t_0`$ be a fixed value of the parameter. Any solution of equation (11), belonging to some class of solutions, has the form $`u(t)=g(t)u_0,`$ where $`u_0M`$ is the initial condition $`u(t_0)=u_0`$ and $`g(t)`$ is a curve on $`G`$ such that $`g(t_0)=e.`$ The curve $`g(t)G`$ is uniquely determined by the condition $`g(t_0)=e,`$ i. e. $`g(t)`$ is the same for all solutions $`u(t)`$ of the equation (11).

In this case we can easily prove in the same way as for ODEs that the difference equation (11) has a superposition formula in the class of solutions mentioned in the Statement 2’. Moreover, this formula will be the same superposition formula (9) as for the differential equation, since in both cases we are working with solutions of the same form $`u(t)=g(t)u_0.`$

Let us now show how this works in the case of the matrix Riccati equation.

## 3 Differential and difference Riccati equations

It is easier in this Section to consider the matrix Riccati equation from the beginning. We will omit the word “matrix” for simplification.

Let us now recall how the Riccati equation arises from the construction described in Section 2. Consider $`G=GL(N)`$ and let $`M`$ be a Grassman manifold $`G_{n+k,k}`$ of $`k`$-planes in an $`N=(n+k)`$-dimensional space. The homogeneous coordinates of a point $`pM`$ are given by the components of an $`(n+k)\times k`$-dimensional matrix

$$\left(\begin{array}{c}X\\ Y\end{array}\right),$$
(12)
where $`X`$ is a $`n\times k`$-matrix and $`Y`$ is a $`k\times k`$-matrix. The columns of (12) span a $`k`$-plane defining $`p.`$ The point $`p`$ is thus identified with the equivalence class $`\left[\left(\begin{array}{c}X\\ Y\end{array}\right)\right]`$ under the relation

$$\left(\begin{array}{c}X\\ Y\end{array}\right)\left(\begin{array}{c}Xh\\ Yh\end{array}\right),hGL(k),$$

identifying different bases for the same $`k`$-plane. The action of an element

$$g=\left(\begin{array}{cc}M& N\\ P& Q\end{array}\right)GL(N)$$

upon $`M`$ is obtained by the projection

$$\pi :\left(\begin{array}{c}X\\ Y\end{array}\right)\left[\left(\begin{array}{c}X\\ Y\end{array}\right)\right]$$

from a linear action

$$g\left(\begin{array}{c}X\\ Y\end{array}\right)=\left(\begin{array}{cc}M& N\\ P& Q\end{array}\right)\left(\begin{array}{c}X\\ Y\end{array}\right).$$

On the affine subspace defined by $`detY0`$ we may define inhomogeneous coordinates $`W=XY^1,`$ $`W`$ is an $`n\times k`$-matrix. In the inhomogeneous coordinates the group action is given by the formula

$$gW=(MW+N)(PW+Q)^1.$$
(13)
Let us now consider Eq. (8). In this case it reduces to Eq. (3), where $`A,B,C,D`$ are defined by $`\xi (t)`$:

$$\xi (t)=\left(\begin{array}{cc}B(t)& A(t)\\ D(t)& C(t)\end{array}\right)𝔤𝔩(N).$$
(14)
Let us now construct a difference version of the Riccati equation. We use the two-parameter class of operators $`U_{q,h}`$ defined by the formula (4) in the Introduction. This class of operators can be described in an axiomatic way by the following properties:

1. $`U_{q,h}\mathrm{\hspace{0.17em}1}=0.`$
2. $`U_{q,h}(af+bg)(t)=aU_{q,h}f(t)+bU_{q,h}g(t),`$ where $`a`$ and $`b`$ are constants.
3. $`[U_{q,h},t]f(t)=f(qt+h).`$

For us the important property of these operators is the modified Leibnitz rule:

$$U_{q,h}(fg)(t)=g(t)U_{q,h}f(t)+f(qt+h)U_{q,h}g(t).$$

Let us proceed as described in Section 2. Let $`G=GL(N)`$ and $`M`$ be the Grassman manifold $`G_{n+k,k}`$ as before, where $`N=n+k.`$ We need to choose a coordinate chart on $`M:`$ let us take the inhomogeneous coordinates described above.

Let $`g(t)`$ be a curve on $`GL(N)`$ such that $`g(t_0)=e`$ and $`W_0`$ be the inhomogeneous coordinates of a point on $`M.`$ We will represent $`g(t)`$ as

$$g(t)=\left(\begin{array}{cc}M(t)& N(t)\\ P(t)& Q(t)\end{array}\right).$$

Let us find a difference equation for $`W(t)=g(t)W_0.`$ The group action is described by the formula (13), so we have

$$W(t)=(M(t)W_0+N(t))(P(t)W_0+Q(t))^1.$$

We can find $`U_{q,h}W(t)`$ directly, but it is better to proceed in the following way. Let us choose homogeneous coordinates for $`W_0,`$ i. e. matrices $`X_0,`$ $`Y_0`$ such that $`X_0Y_0^1=W_0.`$ Then we can consider $`X(t)`$ and $`Y(t)`$ defined by

$$\left(\begin{array}{c}X(t)\\ Y(t)\end{array}\right)=g(t)\left(\begin{array}{c}X_0\\ Y_0\end{array}\right)=\left(\begin{array}{c}M(t)X_0+N(t)Y_0\\ P(t)X_0+Q(t)Y_0\end{array}\right).$$

It is easy to see that $`W(t)=X(t)Y(t)^1.`$ Since the group action in homogeneous coordinates is simply the matrix multiplication, it is easy to see that

$$U_{q,h}\left(\begin{array}{c}X(t)\\ Y(t)\end{array}\right)=\left(U_{q,h}g(t)\right)\left(\begin{array}{c}X_0\\ Y_0\end{array}\right)=\left(U_{q,h}g(t)\right)g(t)^1\left(\begin{array}{c}X(t)\\ Y(t)\end{array}\right).$$

Let us represent $`\left(U_{q,h}g(t)\right)g(t)^1`$ in the following way:

$$\left(U_{q,h}g(t)\right)g(t)^1=\left(\begin{array}{cc}B(t)& A(t)\\ D(t)& C(t)\end{array}\right),$$

in the case of $`U_{1,0}=\frac{d}{dx}`$ the matrix functions $`A(t),\mathrm{},D(t)`$ are the same as in (14). We obtain the following linear difference (or differential in the case $`U_{1,0}`$) equations:

$`U_{q,h}X(t)`$ $`=`$ $`B(t)X(t)+A(t)Y(t),`$ (15)
$`U_{q,h}Y(t)`$ $`=`$ $`D(t)X(t)C(t)Y(t).`$ (16)
It is easy to prove that for two matrices $`A(t)`$ and $`B(t)`$ we have the identity

$$U_{q,h}\left[AB^1\right](t)=\left[U_{q,h}A(t)\right]B(t)^1\left[AB^1\right](qt+h)\left[U_{q,h}B(t)\right]B(t)^1.$$

Using this identity and Eqs. (1516) we can finally find a difference equation for $`W(t):`$

$$U_{q,h}W(t)=U_{q,h}\left[XY^1\right](t)=$$

$$=\left[U_{q,h}X(t)\right]Y(t)^1\left[XY^1\right](qt+h)\left[U_{q,h}Y(t)\right]Y(t)^1=$$

$$=A(t)+B(t)W(t)+W(qt+h)C(t)+W(qt+h)D(t)W(t).$$

Thus we have found Statement 1’ from section 2 in our case:

###### Theorem 2

A function

$$W(t)=(M(t)W_0+N(t))(P(t)W_0+Q(t))^1,$$
(17)
such that

$$g(t)=\left(\begin{array}{cc}M(t)& N(t)\\ P(t)& Q(t)\end{array}\right)GL(N)$$
(18)
and $`g(t_0)=I,`$ is a solution of the difference matrix Riccati equation

$$U_{q,h}W(t)=A(t)+B(t)W(t)+W(qt+h)C(t)+W(qt+h)D(t)W(t)$$
(19)
with the initial condition $`W(t_0)=W_0.`$

The matrices $`A(t),\mathrm{},D(t)`$ are defined by the formula

$$\left(\begin{array}{cc}B(t)& A(t)\\ D(t)& C(t)\end{array}\right)=\left[U_{q,h}g(t)\right]g(t)^1.$$
(20)
In the case of $`U_{1,0}=\frac{d}{dx}`$ we obtain the usual (matrix) Riccati equation.

It is easy to see that the limit of Eq. (19) when $`q1`$ and $`h0`$ is the usual Riccati equation (3) with the same $`A(t),\mathrm{},D(t).`$ Hence we can consider Eq. (19) as a discretization of the usual Riccati equation (3).

As explained in Section 2, in order to prove the existence of a superposition formula we need to prove Statement 2’ from Section 2 in our case, i. e. we need to prove that any solution of Eq. (19) belonging to some class of solutions has the form (17), where the point $`W_0`$ is the initial condition $`W(t_0)=W_0,`$ and the curve $`g(t)`$ (18) is uniquely determined by the condition $`g(t_0)=I,`$ i. e. $`g(t)`$ is the same for all solutions.

As was already explained, we need a uniqueness theorem to prove this. This requires some assumptions on equation (19) and on the class of considered solutions.

## 4 Superposition formulas

The theory of superposition formulas for the discrete Riccati equation (19) depends on the values of $`q`$ and $`h.`$ We will consider two very different cases: the case $`|q|1`$ and the case $`q=1,h0.`$ We will not consider other cases since the case $`q=1,h=0`$ is simply that of differential equations, and the case $`|q|=1,q1`$ would merit a separate treatment.

### 4.1 Superposition formulas in the case $`|q|1`$

The key observation in this case is the following. Let $`T`$ be the transformation $`T:tqt+h`$ of the complex plane. Then the set of points $`T^n(t_0),n,`$ has an accumulation point $`\frac{h}{1q}.`$ The well-known uniqueness theorem for holomorphic functions  says that if a function $`f,`$ holomorphic in a disc, is equal to zero on a subset possessing an accumulation point inside the disc, then $`f`$ is identically zero. This provides us with a uniqueness theorem for solutions of a difference equation in the class of functions meromorphic in the entire complex plane and holomorphic in a neighborhood of the accumulation point. Let us consider two solutions of Eq. (19) from this class. Let us consider a disc such that both solutions are holomorphic in this disc. For simplicity let us consider the case $`|q|<1.`$ Let our solutions be equal in some point $`t_0`$ inside the disc. Then we shall prove that two solutions are equal also in the points $`T^n(t_0),n.`$ It follows from the above mentioned uniqueness theorem for holomorphic functions that two solutions are identical.

On the other hand the class of functions meromorphic in the entire complex plane and holomorphic in a neighborhood of the accumulation point is very natural since if $`g(t)`$ is holomorphic then the solution $`W(t)`$ (17) is, in general, a meromorphic function.

This leads us to the following theorem.

###### Theorem 3

Let $`|q|1.`$ The discrete Riccati equation (19) with holomorphic coefficients $`A(t),B(t),C(t),D(t)`$ has a superposition formula in the class of functions meromorphic in the entire complex plane and holomorphic in a neighborhood of the point $`\frac{h}{1q}.`$

Proof. Let us remark that the change of coordinates $`t=t^{}+\frac{h}{1q}`$ permits us to reduce Eq. (19) to an equation of the same form but with $`h=0.`$ Hence we shall consider only this case, the case of $`q`$-derivative. Let us also remark that the change of coordinates $`t=\frac{1}{q}t^{}`$ permits us to reduce further to an equation of the same form but with $`|q|<1.`$

Let us consider the formula (20) as a difference equation for the $`N\times N`$-matrix $`g(t).`$ We can rewrite it as

$$U_{q,0}g(t)=\left(\begin{array}{cc}B(t)& A(t)\\ D(t)& C(t)\end{array}\right)g(t).$$
(21)
It is a homogeneous linear $`q`$-difference equation for $`g(t).`$ Let us also fix some $`t_0`$ and consider the initial problem $`g(t_0)=I.`$

It can be proven by methods similar to those used in Ref.  that there exists a unique solution $`g(t)`$ holomorphic in the entire complex plane. Let us prove that $`g(t)`$ is non-degenerate, i. e. $`g(t)GL(N).`$ Suppose that in some point $`t_1`$ we have $`detg(t_1)=0.`$ We can rewrite equation (21) as

$$g(qt)=((q1)t\left(\begin{array}{cc}B(t)& A(t)\\ D(t)& C(t)\end{array}\right)+I)g(t).$$

We see that if $`g(t)`$ is degenerate then $`g(qt)`$ is degenerate. It follows that for the holomorphic function $`detg(t)`$ we have $`detg(q^nt_1)=0,n.`$ This implies $`g(t)0,`$ but it contradicts the initial condition $`g(t_0)=I.`$ Hence $`g(t)GL(N).`$

Let $`W_1(t)`$ be a solution of Eq. (19) meromorphic in the entire complex plane and holomorphic in a neighborhood of the point $`0.`$

Let us choose some $`R>0`$ such that the following conditions hold in the disc $`|t|<R:`$

1. $`W_1(t)`$ is holomorphic;
2. The matrix

$$I(q1)t(C(t)+D(t)W_1(t))$$

is non-degenerate.

Such $`R`$ exists since $`C(t),D(t)`$ and $`W_1(t)`$ are holomorphic in zero.

Let us choose some $`t_0`$ such that $`|t_0|<R.`$ Let $`W_0=W_1(t_0)`$ and let

$$g(t)=\left(\begin{array}{cc}M(t)& N(t)\\ P(t)& Q(t)\end{array}\right)GL(N)$$

be the solution of (21) with the initial condition $`g(t_0)=I.`$ Let us consider the function

$$W_2(t)=g(t)W_0=(M(t)W_0+N(t))(P(t)W_0+Q(t))^1.$$
(22)
It follows from Theorem 2 that the function (22) satisfies Eq. (19) and the initial condition $`W_2(t_0)=W_0.`$ Let us note that $`W_2(t)`$ is meromorphic in the entire complex plane.

We can rewrite Eq. (19) in the form

$$W(qt)[I(q1)t(C(t)+D(t)W(t))]=(q1)t(A(t)+B(t)W(t))+W(t).$$

Since the matrix $`I(q1)t(C(t)+D(t)W(t))`$ is invertible in the disc $`|t|<R`$ and $`|q|<1,`$ we have

$$W_1(qt_0)=[(q1)t(A(t_0)+B(t_0)W_0)+W_0][I(q1)t(C(t_0)+D(t_0)W_0)]^1.$$

The same is true for $`W_2(qt),`$ so $`W_2(qt_0)=W_1(qt_0).`$ We can prove in analogous way that $`W_2(q^nt)=W_1(q^nt),n.`$

Since $`W_1(t)`$ is holomorphic in a neighborhood of $`0`$, we have

$$\underset{n\mathrm{}}{lim}W_2(q^nt_0)=\underset{n\mathrm{}}{lim}W_1(q^nt_0)=W_1(0)$$

and this implies that $`W_2(t)`$ cannot have a pole in $`0.`$ Since both $`W_1(t)`$ and $`W_2(t)`$ are holomorphic in some neighborhood of $`0`$ and their values are the same in the sequence of points having $`0`$ as an accumulation point, we obtain $`W_1(t)W_2(t).`$ Thus any solution of Eq. (19) meromorphic in the entire complex plane and holomorphic in a neighborhood of zero has the form (22).

This is “nearly” Statement 2’ from Section 2. “Nearly” because the initial point $`t_0`$ cannot be arbitrary, but should satisfy the condition $`|t_0|<R`$ and this $`R`$ depends on the solution. However it permits us to find a superposition formula. Indeed, let $`W_1(t),\mathrm{},W_m(t)`$ be particular solutions meromorphic in the entire complex plane and holomorphic in a neighborhood of zero. For each $`W_i`$ we can find its own constant $`R_i`$ and take $`R=\mathrm{min}(R_1,\mathrm{},R_m).`$

Let us choose an initial parameter value $`t_0`$ such that $`|t_0|<R.`$ It follows that $`W_1(t),\mathrm{},W_m(t)`$ have the form $`W_i(t)=g(t)W_i(t_0),`$ where $`g(t)`$ is the same for all $`i.`$ As in Section 2 we have the system of equations for $`g(t).`$ If $`m`$ is sufficiently large and $`W_i(t)`$ satisfy some additional independence conditions necessary for $`g(t)`$ to be expressed in terms of $`W_i(t)`$ and $`W_i(t_0),`$ we have

$$g(t)=F(u_1(t),\mathrm{},u_m(t);u_1(t_0),\mathrm{},u_m(t_0)).$$

It follows that the general solution holomorphic in the disc $`|t|<|t_0|+\epsilon `$ has the form

$$u(t)=g(t)u_0=F(u_1(t),\mathrm{},u_m(t);u_1(t_0),\mathrm{},u_m(t_0))u(t_0),$$
(23)
the initial condition $`u(t_0)`$ plays a role of the arbitrary constant in the superposition formula (23). Since the function $`F`$ does not depend on the choice of $`t_0,`$ we can consider a general solution holomorphic in a disc of radius less then $`R.`$ If we decrease $`R`$ and repeat our construction of the superposition formula, the resulting superposition formula will be the same. This finishes the proof. $`\mathrm{}`$

As we already mentioned in Section 2, the resulting superposition formula is the same as in the case of usual Riccati equation (3). These superposition formulas can be found in .

### 4.2 Superposition formulas in the case $`q=1,`$ $`h0`$

By rescaling $`t=ht^{}`$ we can transform Eq. (19) into the same equation with $`q=1,h=1.`$

This case is very different from the case $`|q|1.`$ It is difficult to find a natural class of solutions defined in the entire complex plane. Since the set $`\{t_0+nh,n\}`$ has no accumulation point, the uniqueness theorem for holomorphic functions cannot be applied. In this situation we have, for example, a uniqueness theorem for entire functions of growth order $`1`$ and of normal type, but it is very unnatural. We already mentioned in the beginning of Subsection 4.1 that solutions of the Riccati equations are in general meromorphic. But the situation becomes very natural if we consider functions defined only on $`,`$ i. e. in integer points.

###### Theorem 4

The discrete Riccati equation

$$W(n+1)W(n)=A(n)+B(n)W(n)+W(n+1)C(n)+W(n+1)D(n)W(n)$$
(24)
such that the matrices

$$\left(\begin{array}{cc}I+B(n)& A(n)\\ D(n)& IC(n)\end{array}\right)$$
(25)
are non-degenerate for all $`n`$ has a superposition formula in the class of solutions defined on all $`.`$

Proof. Let us consider the formula (20) as a difference equation for the $`N\times N`$-matrix $`g(t).`$ We can rewrite it as

$$g(n+1)=\left(\begin{array}{cc}I+B(n)& A(n)\\ D(n)& IC(n)\end{array}\right)g(n).$$

Let us also fix some $`n_0`$ and consider the initial problem $`g(n_0)=I.`$ Since the matrix (25) is non-degenerate,

$$g(n)=\left(\begin{array}{cc}M(n)& N(n)\\ P(n)& Q(n)\end{array}\right)$$

is uniquely defined for all $`n`$ and $`g(n)GL(N).`$

Let us consider a solution $`W_1(n)`$ defined on all $`.`$ Let us use the same $`n_0`$ and let $`W_0`$ be the initial condition $`W_0=W_1(n_0).`$ We can consider a function

$$W_2(n)=g(n)W_0=(M(n)W_0+N(n))(P(n)W_0+Q(n))^1.$$
(26)
It is also a solution of Eq. (24) with the initial condition $`W_2(n_0)=W_0.`$

Eq. (24) can be rewritten in the following forms:

$`W(n+1)\left[IC(n)D(n)W(n)\right]`$ $`=`$ $`A(n)+B(n)W(n)+W(n),`$
$`\left[I+B(n)+W(n+1)D(n)\right]W(n)`$ $`=`$ $`W(n+1)A(n)W(n+1)C(n).`$
It follows that if for all $`n`$ we have

$$det[IC(n)D(n)W_1(n)]0,det[I+B(n1)+W_1(n)D(n1)]0,$$
(27)
then for all $`n`$ we can find $`W_1(n)`$ and $`W_2(n)`$ starting from the initial condition $`W_0=W_1(n_0)=W_2(n_0)`$ and we can see that $`W_1(n)W_2(n).`$ We see that all the solutions satisfying the conditions (27) are of the form (26). This implies in the same way as in the proof of Theorem 3 that equation (24) has a superposition formula in the class of solutions defined on all $``$ and satisfying (27). Moreover, a general solution satisfies conditions (27) since they represent only a countable number of inequalities (labeled by $`n`$). $`\mathrm{}.`$

In the case of the differential Riccati equations we have a superposition formula in the class of $`C^1()`$-solutions. However, for some particular choice of the constants this formula gives us solutions outside this class, for example with poles. The same situation arises in the case of equation (24). We have a superposition formula in the class of solutions defined on all $`,`$ but for some particular choices of the constants this formula gives us solutions which are singular in some points.

For example, let us consider

$$g(t)=\left(\begin{array}{cc}1+t^2& t^2\\ t^2& 1t^2\end{array}\right)GL(2).$$

This curve gives us solutions

$$w(t)=\frac{(u_01)t^2+u_0}{(u_01)t^2+1}$$

of the equation

$$\dot{w}(t)=2t+4tw(t)2tw(t)^2$$

with initial condition $`w(0)=u_0.`$ For $`u_01`$ these solutions belong to $`C^1()`$ and we can write a superposition formula. If we take three particular solutions $`w_0,w_1,w_2`$ corresponding to $`u_0=1,u_0=2`$ and $`u_0=3,`$ respectively, we obtain from (2) the formula for the general solution belonging to $`C^1():`$

$$w(t)=\frac{2t^2+4c}{2t^2+2c}.$$

If, for example, $`c=4,`$ we obtain a solution

$$w(t)=\frac{t^2}{t^21},$$

which is singular at $`t=\pm 1.`$

In the difference case the same curve gives us solutions

$$w(n)=\frac{(u_01)n^2+u_0}{(u_01)n^2+1}$$

of the equation

$$w(n+1)w(n)=2n1+(2n+1)w(n)+(2n+1)w(n+1)(2n+1)w(n+1)w(n).$$

We obtain the same superposition formula

$$w(n)=\frac{2n^2+4c}{2n^2+2c}$$

valid for the solutions defined on all $`.`$ It also gives us solutions with singularities. For example for $`c=4`$ we have

$$w(n)=\frac{n^2}{n^21},$$

which is singular at $`n=\pm 1.`$ We see that the superposition formula can be used to avoid singularities in numerical computations, as already pointed out in Ref. .

## 5 Conclusions

We have constructed very natural discrete matrix Riccati equations with superposition formulas. This was done for a family of difference operators general enough to include both $`q`$-derivatives and standard discrete derivatives.

The key idea was considering superposition formulas in an appropriate class of solutions. Which class depends on the type of difference operator used. It turns out that in the case of the $`q`$-derivative it is natural to consider solutions meromorphic in the entire complex plane and holomorphic in a neighborhood of zero, and in the case of the standard discrete derivative, it is natural to consider solutions defined only in integer points.

The discretization of Riccati equations of this article should be compared with a previously proposed discretization of linearizable equations, including the Riccati one . There the curve in a Lie algebra, i. e. $`\xi (t)`$ of Eq. (8), was replaced by a curve in the corresponding Lie group $`G`$ and this gave a linearizable mapping. E. g. for the Riccati equation (1), the corresponding discretization would be

$$w(t+h)=[g_{11}(t)w(t)+g_{12}(t)][g_{21}(t)w(t)+g_{22}(t)]^1.$$
(28)
Here $`h`$ is, on one hand, a lattice spacing, on the other hand a group parameter. Thus, the whole class of Riccati equations is replaced by the above class of homographic mappings. If equation (28) has constant coefficients, then the relation between the Riccati mapping and the Riccati equation is simple:

$$g=\left(\begin{array}{cc}g_{11}& g_{12}\\ g_{21}& g_{22}\end{array}\right)=e^{Lh},Lsl(2),gSL(2),$$

i. e.: $`g_{11}=1+hb+\mathrm{},`$ $`g_{12}=ha+\mathrm{},`$ $`g_{21}=hd+\mathrm{},`$ $`g_{22}=1hc+\mathrm{},`$ and for $`h0`$ (28) reduces to

$$\dot{w}(t)=a+(b+c)w(t)+dw^2(t).$$

For variable coefficients the relation between a specific Riccati mapping and a specific Riccati equation is difficult to establish (i. e. it is necessary to integrate the Riccati equation explicitly). For all details we refer to the original articles . We also mention that discrete integrable mappings of the type (28) have been used to discretize the Pinney equation .

The relation between the discrete and continuous equations of this article, on the other hand, is quite straightforward. Given the ODE (3), we give the discrete equation (19) explicitly, involving the same functions $`A(t),\mathrm{},D(t).`$ Thus equation (19) can serve as the basis for numerical approximations for differential equations. This would complement a previous use of superposition formulas in numerical analysis .

We would like to comment on the relation between the discrete Riccati equations and the Runge-Kutta method. Statement 2 of Section 2 states that in the class of $`C^1()`$-solutions the differential Riccati equation (3) is equivalent to the linear matrix equation

$$\dot{g}(t)=\left(\begin{array}{cc}B(t)& A(t)\\ D(t)& C(t)\end{array}\right)g(t)$$
(29)
with the initial condition $`g(t_0)=I.`$ Applying the Runge-Kutta method to Eq. (29), we obtain the difference equation

$$U_{1,h}g(t)=\left(\begin{array}{cc}B(t)& A(t)\\ D(t)& C(t)\end{array}\right)g(t)$$

which is in a similar way equivalent to the difference Riccati equation defined by the operator $`U_{1,h}.`$ As it was explained before this will also give us singular solutions.

Equivalently we can say that all $`C^1()`$-solutions of the Riccati equation (3) have the form $`W(t)=X(t)Y(t)^1,`$ where $`X(t)`$ and $`Y(t)`$ are solutions of the linear system

$$\frac{d}{dt}\left(\begin{array}{c}X(t)\\ Y(t)\end{array}\right)=\left(\begin{array}{cc}B(t)& A(t)\\ D(t)& C(t)\end{array}\right)\left(\begin{array}{c}X(t)\\ Y(t)\end{array}\right).$$
(30)
Applying the Runge-Kutta method to Eq. (30) gives us the equation

$$U_{1,h}\left(\begin{array}{c}X(t)\\ Y(t)\end{array}\right)=\left(\begin{array}{cc}B(t)& A(t)\\ D(t)& C(t)\end{array}\right)\left(\begin{array}{c}X(t)\\ Y(t)\end{array}\right).$$
(31)
Introducing the inhomogeneous coordinates (the matrix elements of $`W=XY^1`$) we transform (31) into the discrete Riccati equation (19) with $`q=1,h0.`$

In other words, the $`h`$-discretization of the matrix Riccati equation, presented in this article, not only preserves the superposition formula. It is also equivalent to applying the Runge-Kutta method to the associated linear system (30).

## Acknowledgments

A. P. is grateful to the Centre de recherches mathématiques, Université de Montréal, for its hospitality. The research of P. W. was partly supported by research grants from NSERC of Canada, and FQRNT du Québec.


format: text

text:
# A System with Multiple Liquid-Liquid Critical Points

## Acknowledgments

We thank V. V. Brazhkin, A. Geiger, G. Franzese, T. Loerting, and A. Skibinsky for helpful discussions and the NSF Chemistry Division (CHE-0096892) for support.


text:
# Complex Band Structures and Decay Length in Polyethylene Chains

## I Introduction

The conductance of interfaces, tunnel junctions or molecular chains connecting two tips is of paramount importance for areas such as nanoelectronics or molecular electronics. It is well known that in order to understand the electronic and transport properties of surfaces and interfaces, the complex electronic band structure (CBS) of the constituent materials  maybe as important as the real band structure. However, ab-initio calculation of CBS are still rare.

Very recently the Bloch functions with complex wavevector $`𝐤`$ have been revived for the calculation of the ballistic conductance of carbon nanotubes,  because they provide an effective way to deal with open quantum systems and transmission coefficients within the Landauer theory. Furthermore, the spin dependent conductance of ferromagnetic/insulator/ferromagnetic junctions has been explained in relation with the complex bands of the insulator.  Inside an isolated surface or tip, the electronic wave functions can be written as a superposition of Bloch functions, including both real and complex $`𝐤`$, namely both propagating and evanescent waves. The vacuum tunneling between two tips arises as the overlap in the vacuum gap of two opposite evanescent waves. The electronic states of a nanowire, a nanotube, a long polymer bridged between two tips, or of an insulating thin film sandwiched between two metallic contacts, can all be expanded in Bloch states with both real and complex $`𝐤`$. A detailed analysis of the CBS itself helps to explain the properties of the tunneling current. For instance, in the theory of Scanning Tunnelling Microscope (STM), the exponential decay of the current with the tip surface distance, depends on the metal work function which controls the decay length of the evanescent metal states in the vacuum.

This concept was applied in Ref.  to calculate the non-resonant tunneling parameter $`\beta `$ of n-alkane chains measured accurately in recent experiments.  The resistance of alkanethiol chains of variable monomer numbers was measured by adsorbing the n-alkane chains on a gold (111) surface and contacting the other extreme with an Atomic Force Microscopy (AFM) tip coated with gold. The pressure on the chains is kept constant by controlling the force exerted by the tip on the n-alkane film. The resistance depends strongly on the contact force, but at fixed force, it depends exponentially on the chain length ($`R=R_0e^{\beta N}`$, where $`N`$ is the number of CH2 monomers). The reported values of $`\beta `$ range from 0.8 to 1.2 per monomer, depending on the contact force. Note that in vacuum, $`\beta =\sqrt{\frac{8m}{\mathrm{}^2}\mathrm{\Phi }_{Au}}a=3.0`$ where ($`a=1.25`$ Å is the length of a monomer, $`\mathrm{\Phi }_{Au}=5.3`$ eV is the gold work function, $`m`$ is the electron mass, and $`\mathrm{}`$ the Planck’s constant). The metal wavefunctions decay exponentially in the polymer and, at low tip-surface voltage, the current is due to non-resonant tunneling because the Fermi energy of the gold contacts lies within the gap of the n-alkane chain.  The decay length $`\beta =2(\mathrm{Im}k)a`$ can be derived from the CBS of a polyethylene chain (the polymeric form of the n-alkane). For sufficiently long chains, only the slowest decaying state dominates the tunneling current.  In the CBS, this state belongs to a loop that joins, at imaginary $`k`$, the highest occupied molecular orbital (HOMO) and an antibonding unoccupied molecular orbital (that we call here LUMO) of the polyethylene chain. Taking the Fermi level at the point where $`\frac{dϵ}{dk}=0`$ a value of $`\beta `$ close to 1.0 was found in literature,  based on CBS calculations done with two different local basis sets (a minimal one, $`s`$ for H and $`sp^3`$ for C, and an extended one, $`sp^3`$ for H and $`sp^3d^5`$ for C). There, the decay length $`\beta `$ was well converged with respect to the basis size, but the imaginary bands show a huge basis set dependence. This is the aspect that we address here. A local basis should work very well for the localized bonding and antibonding states of molecules, but problems could arise when considering saturated, wide gap molecular wires such as n-alkane chains. Here the antibonding empty states often fall above the zero of the electrostatic potential in the vacuum. It is known that, experimentally  as well as in DFT,  polyethylene is a negative electron affinity material. While holes in the highest valence band are intra-chain bonding states, electrons in the lowest conduction band are delocalized inter-chain states while the empty antibonding states lie above the vacuum level. In a monolayer of n-alkane chains attached to the gold surface, both intrachain and interchain states appear. The latter cannot easily be described and are sometimes missing altogether in local basis calculations. Ignoring these interchain states or vacuum states in conductante calculations may still work or it may not: our point is that they cannot be a priori neglected. Vacuum states between the HOMO and the LUMO, in general, will influence the complex band structure of the polymer in a way that needs to be addressed.

To clarify this point, we decided to study the complex electronic structure of an isolated periodic one-dimensional polyethylene chain and to compare the CBS obtained with a plane wave basis with those given by a local basis. Our results confirm a very strong dependence of CBS upon the choice of the basis. In the plane wave basis, the decay length for non-resonant tunneling is shown to be close to that obtained by the local basis, but, connecting to the free-electron-like states present in the middle of the HOMO-LUMO gap, at imaginary $`k`$, a continuum of free-electron-like parabolic states appears. It is very difficult to describe these delocalized states with a local basis and in general they are substituted by a few basis dependent states. Depending on the chemical potential, they might be relevant or not to conductance. For instance, they represent vacuum tunneling processes for very short nanocontacts. Whatever the case, in order to discuss the tunneling current, it is useful to keep in mind that they exists and that they might be important in some cases.

This paper is organized as follow. In Section II, we describe the method used to get the CBS of the polyethylene with a plane wave basis.  In Section III, we describe the geometry of the isolated chain and of the bulk polyethylene used in the calculation and the other technical details. In Section IV, we present our calculated CBS, compare with previous calculations, and discuss the consequences of our results on the tunneling current of n-alkane chains. Finally, Section V contains our conclusions.

## II Method

In this Section, we describe briefly the method used to compute the CBS of a periodic system. The Bloch theorem allows us to label the one electron states with a wavevector $`𝐤`$ and a band index $`n`$. We are interested in Bloch functions with both real and complex $`𝐤`$. Within DFT in the local density approximation (LDA) these Bloch functions are solutions of the one-electron Kohn-Sham equations. To deal with ultrasoft pseudopotentials,  we need to consider the generalized form of these equations  (atomic units $`e^2=2m=\mathrm{}=1`$ are used):

$$E\widehat{S}|\psi _𝐤=\left[^2+V_{\mathrm{eff}}+\widehat{V}_{NL}\right]|\psi _𝐤,$$
(1)
where $`V_{\mathrm{eff}}`$ is the effective local potential (see Ref. ), $`V_{NL}`$ is the nonlocal part of the ultrasoft pseudopotential:

$$\widehat{V}_{NL}=\underset{Imn}{}D_{mn}^I|\beta _m^I\beta _n^I|,$$
(2)
constructed using the set of projector functions $`\beta _m^I`$ associated with atom $`I`$ and $`\widehat{S}`$ is the overlap operator:

$$\widehat{S}=1+\underset{Imn}{}q_{mn}^I|\beta _m^I\beta _n^I|.$$
(3)
The coefficients $`q_{mn}^I`$ are the integrals of the augmentation functions defined in Ref. . It is convenient to rewrite Eq. (1) in the form:

$$E|\psi _𝐤=\left[^2+V_{\mathrm{eff}}\right]|\psi _𝐤+\underset{Imn}{}(D_{mn}^IEq_{mn}^I)|\beta _m^I\beta _n^I|\psi _𝐤.$$
(4)
In this equation we take the energy $`E`$ as a fixed external parameter and the values of $`𝐤`$ (in general complex numbers) are determined so that $`|\psi _𝐤`$ obeys the Bloch periodicity condition. This problem, formulated within the ultrasoft pseudopotential scheme is not substantially different from the norm conserving one solved in Ref. .

We solve this equation in a supercell geometry along the $`x`$ and $`y`$ directions. A wave propagating at energy $`E`$ will have $`𝐤=(𝐤_{},k)`$. Due to the supercell geometry $`\psi _{(𝐤_{},k)}`$ will have the Bloch form along $`x`$ and $`y`$:

$$\psi _{(𝐤_{},k)}(𝐫_{}+𝐑_{},z)=e^{i𝐤_{}𝐑_{}}\psi _{(𝐤_{},k)}(𝐫_{},z),$$
(5)
where $`𝐤_{}`$ is a real vector whereas $`k`$ can be complex. Different $`𝐤_{}`$ do not mix and can be considered separately. From now on, we omit the $`𝐤_{}`$ on $`\psi _{(𝐤_{},k)}`$.

As in Ref.  we write the solution of the integro-differential Eq. (4) as a linear combination:

$$\psi _k(𝐫)=\underset{n}{}a_{n,k}\psi _n(𝐫)+\underset{Im}{}c_{Im,k}\psi _{Im}(𝐫),$$
(6)
where $`\psi _n`$ are linearly independent solutions of the homogeneous equation:

$$E\psi _n(𝐫)=\left[^2+V_{\mathrm{eff}}(𝐫)\right]\psi _n(𝐫),$$
(7)
and $`\psi _{Im}`$ is a particular solution of the inhomogeneous equation:

$$E\psi _{Im}(𝐫)=\left[^2+V_{\mathrm{eff}}(𝐫)\right]\psi _{Im}(𝐫)+\underset{𝐑_{}}{}e^{i𝐤_{}𝐑_{}}\beta _m^I(𝐫\tau \tau \tau ^I𝐑_{}).$$
(8)
Both $`\psi _n`$ and $`\psi _{Im}`$ are $`(x,y)`$ periodic as in Eq. (5). Summation over $`Im`$ in Eq. 6 involves all the projectors in the unit cell and the coefficients $`c_{Im,k}`$ are determined by:

$$c_{Im,k}=\underset{n}{}[D_{mn}^IEq_{mn}^I][\beta _n^I(𝐫\tau \tau \tau ^I)]^{}\psi _k(𝐫)d^3r.$$
(9)
Note that the set of wave functions $`\psi _n`$ is infinite. In practice the expansion $`\psi _n(𝐫)=_𝐆_{}\psi _n(𝐆_{},z)e^{i(𝐤_{}+𝐆_{})𝐫_{}}`$ over two-dimensional plane waves is used. If one considers only $`N_{2D}`$ plane waves with $`𝐆_{}^2E_{\mathrm{cut}}`$ the number of $`\psi _n`$ becomes finite and equals to $`2N_{2D}`$. As in Ref., to find the functions $`\psi _n`$ we discretize the unit cell along the $`z`$ axis by dividing it into slices and replacing the $`V_{\mathrm{eff}}(𝐫)`$ in each slice by a $`z`$-independent potential. The $`\psi _n`$ in the slices is a linear combination of two exponentials with coefficients obtained by a matching procedure.

The allowed values of $`k`$ for a given energy $`E`$ can be determined by imposing periodicity along $`z`$ to the Bloch function and to its $`z`$ derivative. If $`d`$ is the cell length along $`z`$ we have:

$$\psi _k(𝐆_{},d)=e^{ikd}\psi _k(𝐆_{},0).$$
(10)
$$\psi _k^{}(𝐆_{},d)=e^{ikd}\psi _k^{}(𝐆_{},0).$$
(11)
Inserting Eq. 6 into Eqs. 9-11 one can show that the last three equations are equivalent to the generalized eigenvalue problem:

$$AX=e^{ikd}BX,$$
(12)
where $`A`$ and $`B`$ are known matrices. We solve this problem to obtain in general a complex $`k`$ and coefficients $`X=\{a_{n,k},c_{Im,k}\}`$ for the generalized Bloch state $`\psi _k`$ at a given energy $`E`$ and $`𝐤_{}`$.

Our CBS calculations proceed in two steps. First we compute with a standard electronic structure code (PWSCF)  the ground state electronic structure of the system and we obtain the effective local potential $`V_{\mathrm{eff}}(𝐫)`$ and the screened coefficients $`D_{mn}^I`$. In a second step, we use the potential $`V_{\mathrm{eff}}(𝐫)`$ and the screened coefficients $`D_{mn}^I`$ to compute the CBS applying the scheme described above.

## III Geometrical model and technical details

The n-alkane chains \[X(CH2)nY\] are made of CH2 monomers with singly bonded carbon atoms and two terminating groups (indicated with X and Y). If grafted via a thiolic terminating group (X=SH) on the Au(111) surface, they can spontaneously assemble into a dense and ordered monolayer film to form an insulating layer.  We model the n-alkane chain with a periodic polyethylene chain \[poly-CH2\] as depicted in Fig.1a. Since we shall compare our results with Ref. , we choose the same geometry, C-C-C and H-C-H angles are those of the perfect tetrahedral structure (109.5) and the C-C bond lengths and the C-H bond lengths are respectively equal to 1.54 Å and 1.10 Å.

The chain is modelled with a tetragonal supercell. The chain is parallel to the $`z`$ axis and periodically repeated along the $`x`$ and $`y`$ directions with a chain-chain distance $`a`$ ($`10.58`$ Å) large enough to avoid replica interactions (tests with a larger supercell gave no relevant change of the results). The length of the unit cell $`c`$ in the $`z`$ direction is $`c=2.51`$ Å(experimental $`c=2.50`$ Å). In order to estimate the influence of the chain-chain interaction, we also calculated bulk polyethylene. The rombohedral unit cell edges are $`a=4.93`$ Å, $`b=7.40`$ Å, $`c=2.534`$ Å, as in Ref. , and the internal coordinates were allowed to relax. We found $`d_{CC}=1.51`$ Å, $`d_{CH}=1.11`$ Å, $`\alpha _{CCC}=114^{}`$ and $`\alpha _{HCH}=105^{}`$ (See Fig. 1b) to be compared with experimental values $`d_{CC}=1.53`$ Å, $`d_{CH}=1.09`$ Å, $`\alpha _{CCC}=113.3^{}`$ and $`\alpha _{HCH}=109.5^{}`$.

We calculated the electronic structure of these systems with the PWSCF code  using the local density approximation (LDA) with the parameterization of Perdew and Zunger.  The electron-ion interaction was described by ultrasoft pseudopotentials.  The wavefunctions were expanded in a plane waves basis with an energy cut-off of $`25`$ Ry. The cut-off for the charge density is $`150`$ Ry. The integrations on the Brillouin zone are done with a uniform mesh of $`40`$ $`k`$ points for the isolated chain and with a mesh $`3\times 3\times 10`$ $`𝐤`$ points for the bulk.

The CBS calculation was subsequently carried out with the technique described above. The unit cell length $`c`$, in the $`z`$ direction, is divided into 15 slices and the integrals are performed with 11 points in each slice (the slice thickness is 0.167 Å which is sufficient to describe accurately the slowly varying effective potential). The two dimensional cut-off energy is taken as in the electronic structure calculation. We consider only the $`𝐤_{}=\mathrm{\Gamma }`$ point. This is a very good approximation for the isolated chain, while it gives the complex band structure of bulk polyethylene only along one line in reciprocal space.

## IV Results and discussion

The CBS calculated for the periodic isolated polyethylene chain and for the bulk polyethylene are shown in Fig. 2a and 2b respectively. The lowest unoccupied state, identifiable with the vacuum zero in the limit of infinite interchain spacing, is taken as the zero of the energy. In each Figure, the right panel shows the bands at the real $`k`$, the Bloch modes propagating along the $`z`$ direction, whereas the left panel displays the bands obtained for complex $`k`$, Bloch modes evanescent (or exponentially diverging) in the $`z`$ direction (note that for each solution $`k`$, at a given energy $`E`$, also the complex conjugate $`k^{}`$ is a solution). Bands with purely imaginary $`k`$ are represented with small dots, while bands with complex $`k`$ are projected on the imaginary $`k`$ plane and indicated by small circles.

The energy bands at real $`k`$ of polyethylene (both for the bulk structure and for an isolated chain) have been reported in several papers.  The valence bands, C-C and C-H bonding states, have been also compared with experimental angular-resolved photoemission data.  In Table I, we compare our results with previous DFT calculations. All theoretical band widths are in substantial agreement with each other, weakly dependent on the basis set used or the actual geometry (isolated chain or bulk polyethylene). Our conduction band widths are 14.0 eV (chain) and 13.6 eV (bulk) in good agreement with previous calculations. All these results are also in reasonable agreement with the experimental value of 16.2 eV. In all theoretical calculations, however, the absolute position of the occupated bands of the chain with respect to the vacuum level, is too high if compared with experiment. In the bulk, according to the photoemission data of Ref. , the HOMO level lies 10 eV below vacuum (8.8 eV below vacuum, according to Ref. ). DFT calculations however place the valence band top about 6.0 eV below vacuum. Thus, for example, a shift of 4.5 eV was needed in Ref.  to bring theory and experiment in closer agreement.

The theoretical band gaps appear to be in disagreement with each other. For the chain there are two sets of values: depending on the exact geometry or exchange and correlation functional, gaps calculated with plane waves cluster around 5-6 eV,  while gaps calculated with local basis are around 8-10 eV.  We find an energy gap of 5.14 eV (for the isolated chain) which becomes 6.4 eV (in bulk polyethylene) in agreement with previous plane waves results.  Süle et al.  have found that the band gap depends on the basis-set quality. Variations of several electron volts have been attributed to the appearance in the gap of new states belonging to the continuum when the basis-set size is increased. The nature of the conduction band states obtained with plane waves has been discussed in detail by Serra et al.. In bulk polyethylene, the charge density associated to the $`\mathrm{\Gamma }`$ point wavefunction at the conduction band minimum shows clear maxima between the chains. Increasing the interchain distance, this state evolves to a free-electron-like state, or perhaps a surface state.  In addition a slab calculation has shown that, within GGA, polyethylene is a negative electron affinity material.  According to experiment (Ref. ) the electron conduction band lies at an energy $`A`$ above the vacuum level, and $`A`$ is about 0.5 eV above the vacuum. It is a characteristic feature of DFT in the LDA approximation to yield gaps that are generally too small, with the exception of bonding anti-bonding gaps. Comparing our complex bands with those obtained via a localized basis,  it is possible to identify the reason of the different theoretical values for the band gap. In the chain, the gap calculated by plane waves measures the position of the HOMO with respect to the vacuum level, while the gap calculated by local basis is very close to the energy of the HOMO with respect to an antibonding empty state (that we call LUMO). Actually free-electron-like states are poorly described by a local basis set and they barely appear within the HOMO-LUMO gap. Instead in our calculation, the lowest conduction bands are parabolic, indicating free-electron-like states which, in bulk polyethylene, indicate interchain states.  In a truly isolated polyethylene chain the vacuum zero energy identifies the onset of free-electron-like states. The LUMO, which is higher in energy than the vacuum zero should appear as a resonant level in the density of states at positive energy. Since we are using a supercell geometry the density of free-electron-like states will however depend on the box size.

We now analyze the band structure at imaginary $`k`$ and at complex $`k`$ of the isolated chain and compare our results with those of Ref. . At imaginary $`k`$ the CBS is formed mainly by free-electron-like parabolas. The presence of the polyethylene chain splits the degeneracy of free-electron bands. In our calculation these bands are discrete, but at infinite box size, above the free-electron parabola departing from the $`\mathrm{\Gamma }`$ point at vacuum zero energy, they are a continuum of imaginary $`k`$ states. The continuum of parabolic bands is not well reproduced with a local basis and actually the states that one obtains at imaginary $`k`$ depends on the basis as shown in Ref.  where a minimal and an extended basis give totally different results. At the $`\mathrm{\Gamma }`$ point parabolic bands depart also from all the valence states except the HOMO. A loop starts from the HOMO and joins the LUMO of the molecule which is at $`9.3`$ eV above the HOMO, 4.13 eV higher, in DFT, than the vacuum. We note that our HOMO-LUMO gap is correctly close to the local basis gap.

For complex $`k`$ there are loops joining maxima and minima of the bands. The loop between $`12.3`$ eV and $`9.9`$ eV closes the gap between the C-C $`\sigma `$-bonding band and the C-H bonding states. Since these states are well reproduced by local basis, the shape and the size of this loop is in good agreement with Ref. . Two additional loops, also in good agreement with the local basis calculation, join maxima in the valence bands with minima in the conductions states. In the energy gap, these states have a decay length shorter than the loop joining the HOMO and the LUMO of the molecule.

The imaginary and complex $`k`$ bands of bulk polyethylene are shown in Fig. 1b. From the comparison of this result with the isolated chain, we can estimate the effect of the reduced interchain distance on the complex bands. As expected, the density of free-electron-like states is much reduced. Furthermore, since there are four monomers per unit cell, all the bands are doubled with small splittings in the valence bands and much larger splittings in the conduction bands. The loop which, in the isolated chain, joins the HOMO and the LUMO at imaginary $`k`$ is now split in two loops through an avoided crossing with a parabolic band. The two loops that, at complex $`k`$, join valence and conduction bands are both split and, since the gap is now indirect, in a small energy window below the LUMO, the states in one of these loops have a longer decay lenght than those in the loop which joins the valence and conduction band at the $`\mathrm{\Gamma }`$ point.

In a typical experiment where the tunneling decay length $`\beta `$ of an n-alkane chain is measured, the molecule is adsorbed on an Au(111) substrate by terminating it with a thiol group. As shown in Ref.  the Fermi level of the gold substrate is between the HOMO-LUMO gap. Its real position is not known experimentally; theoretically it depends on the basis set used in the calculation. With plane waves, it is at about 3.5 eV  above the HOMO of the polymer. It is reasonable to suppose that the presence of the thiol group, as well as the termination of the chain are perturbations which decay much faster than the metal tunneling states. With this hypothesis, the typical decay length of the tunneling current can be approximated with the longest decay length of the evanescent Bloch states of the polyethylene chain at the Fermi energy. The important complex band determining the tunneling decay length is therefore the loop joining the HOMO and LUMO states since it has the smallest $`\mathrm{Im}k`$ at the Fermi level. Our value for $`\mathrm{Im}k`$ is 0.36 Å-1, corresponding to $`\beta =0.9`$ per monomer, in good agreement with experiment and with the calculation of Refs. , the small difference with Ref.  being due only to the choice of the Fermi level. Our CBS shows that local basis and plane waves give the same decay lenght for the non-resonant tunneling regime if the chemical potential is much lower than the vacuum. According to a local basis calculation, increasing the substrate tip voltage, one would however expect to move the chemical potential of the metal and to enter in a resonant tunneling regime where the current pass through the LUMO of the molecule.  According to our plane wave LDA calculation there is instead a range of voltages before the resonant tunneling regime, where vacuum tunneling should start to contribute to the current. The LDA underestimates the voltage necessary for this phenomenum to occur, since the HOMO is only 5.14 eV below the vacuum, while the correct gap should be of 8.8 eV but, in negative electron affinity molecules, vacuum states do start to play a role before the LUMO. The LDA error in positioning the HOMO of the molecule must be kept in mind before comparing with experiment the conductance calculated with a plane wave basis. Nonetheless, negative electron affinity molecules such as the n-alkanes must be treated carefully expecially with local bases because the incorrect description of vacuum states may lead, not only to a complete lack of free-electron-like states in the HOMO-LUMO gap, but also to the presence of many fictitious states at imaginary $`k`$, that are totally dependent on the choice of the basis set.

## V Conclusion

In this paper, the CBS of a periodic polyethylene chain has been calculated using a plane wave basis set and compared with the CBS obtained with a local basis.  We have shown that the tunneling decay length $`\beta `$ of n-alkane chains, here replaced by the polyethylene, is independent of the choice of the basis, at least as long as at low voltages the tunneling is dominated by a slowly decaying state which belongs to the loop joining, at imaginary $`k`$, the HOMO and LUMO states. Instead, our results show that the choice of the basis could be critical in large gap polymers with negative electron affinity particularly at large voltages. Future calculations of the conductance must be done carefully because of the presence of free-electron-like vacuum states within the HOMO-LUMO gap, signifying vacuum tunneling besides intra-chain tunneling. With a plane wave basis, the free-electron-like modes are taken into account but the density functional approximation puts the HOMO too close to the vacuum. With a local basis that treats correctly the bonding and antibonding states of the polymer the density functional problem remains, and on top of that the free-electron-like states are replaced by basis dependent levels.

###### Acknowledgements.

This project is sponsored by MIUR COFIN2001, INFM/F, INFM/G, Iniziativa Trasversale Calcolo Parallelo, and by EU through TMR FULPROP, Contract ERBFMRXCT970155. Part of the calculations have been done on the IBM SP3 at Cineca in Bologna.


text:
# Resonant nature of phonon-induced damping of Rabi oscillations in quantum dots

## Abstract

Optically controlled coherent dynamics of charge (excitonic) degrees of freedom in a semiconductor quantum dot under influence of lattice dynamics (phonons) is discussed theoretically. We show that the dynamics of the lattice response in the strongly non-linear regime is governed by a semiclassical resonance between the phonon modes and the optically driven dynamics. We stress the importance of the stability of intermediate states for the truly coherent control.

One of the most challenging proposals exploiting the atomic-like properties of quantum dots (QDs) jacak98a ; bimberg99  is to implement the quantum information processing schemes bouwmeester00 ; alber01  using confined carrier states in these artificial systems. One of the proposed solutions is to use the charge degrees of freedom which may be efficiently controlled optically even on sub-picosecond scale biolatti00 ; derinaldis02 ; chen01 ; piermarocchi02  which may look promising in view of the measured 1 ns lifetime of excitons in QDs borri01 ; bayer02 . The first steps towards this goal include experimental proofs of controlled coherent dynamics in these structures bonadeo98 , observation of Rabi oscillations stievater01 ; kamada01 ; htoon02 ; zrenner02 ; borri02a , demonstration of entanglement between states of interacting dots bayer01  and performing a quantum logic gate based on a bi-exciton system in a QD li03 .

However, the coherent dynamics of quantum confined carrier states is very sensitive to interaction with the macroscopic number of degrees of freedom of the outside world. If the coherent dynamics involves an excited state of the confined exciton bonadeo98 ; kamada01 ; htoon02  the coherence time is limited by relaxation to the ground state (usually $`30`$ ps). The measured signal may further be affected by the particularities of the detection technique htoon02 . Another approach consists in measuring the dot occupation upon driving by pulses with constant length and varying amplitudes stievater01 ; zrenner02 ; borri02a . So far, it has always turned out that such pulse-area-dependent Rabi oscillations deviate from the ideal ones, the discrepancy being larger for stronger pulses. In principle, this might also be explained by experimental conditions or environmental perturbation: scattering by weakly localized excitons around an interface fluctuation QD (further confirmed by increasing decay for stronger pulses) stievater01 , tunneling to leads in the photodiode structure (on $`10`$ ps timescale) zrenner02 , or dipole moment distribution in the QD ensemble borri01 .

One might believe that all the perturbation comes from sources that may be removed or minimized by technology improvement and by optimizing the experimental conditions and hence produce no fundamental obstacle to arbitrarily perfect quantum control over the excitonic states. However, in every case the QDs are inherently coupled to the surrounding crystal lattice. The perturbing effect of lattice modes (phonons) has been observed experimentally as a fast ($`1`$ ps) partial decay of optical polarization induced by an ultra-fast laser pulse borri01 . The theoretical analysis krummheuer02 ; vagov02a ; jacak03b  shows that the experimentally observed effects may be quantitatively accounted for by invoking the carrier-phonon coupling. Indeed, under suitable excitation conditions in a high-quality sample the carrier density may be kept low enough to eliminate Coulomb scattering effects while the coupling to the electromagnetic modes is manifested by the exciton radiative decay on timescales many order of magnitude longer than these relevant here. The theory shows that the coherence decay should be viewed as a trace of coherent lattice dynamics (due to lattice inertia), persisting even at zero temperature when the lattice is initially at its ground state and cannot perturb the stationary ground state of the carrier subsystem. Due to strong reservoir memory on timescales relevant for these processes, they cannot be fully understood within the Markovian approximations. In particular, the idea of a “decoherence time”, with which the control dynamics competes, is misleading alicki02a .

The recent theoretical study forstner03  on optical Rabi flopping of excitons in QDs driven by finite-length optical pulses shows that exponential damping models fail to correctly describe the system kinetics. The appropriate description yields much less damping, especially for long pulses. It turns out that the lowest “quality” of the Rabi oscillation is obtained for pulse durations of a few ps, while for longer durations the damping is again decreased.

In this paper we propose both qualitative and quantitative explanation of the mechanism leading to the phonon-induced damping reported in the theoretical and experimental studies. We show that the carrier-phonon interaction responsible for the damping of the oscillations has a resonant character: While in the linear limit the system response depends only on the spectral decomposition of the pulse, the situation is different when a strong pulse induces an oscillating charge distribution in the system. In a semi-classical picture, this would act as a driving force for the lattice dynamics. If the induced carrier dynamics is much faster than phonon oscillations the lattice has no time to react until the optical excitation is done. The subsequent dynamics will lead to exciton dressing, accompanied by emission of phonon packets, and will partly destroy coherence of superposition states borri01 ; krummheuer02 ; vagov02a ; jacak03b  but cannot change the exciton occupation number. In the opposite limit, the carrier dynamics is slow enough for the lattice to follow adiabatically. The optical excitation may then be stopped at any stage without any lattice relaxation incurred, hence with no coherence loss. The intermediate case corresponds to modifying the charge distribution in the QD with frequencies resonant with the lattice modes which leads to increased interaction with phonons and to decrease of the carrier coherence (see Ref. axt99,  for a simple, single-mode model).

The above ideas are supported by the formal analysis in the following. First, we describe the theoretical method used; then the description is applied to a specific self-assembled InAs/GaAs quantum dot to get quantitative estimations; finally, we conclude with some summary and possible consequences of the results.

Following the earlier studies krummheuer02 ; vagov02a ; vagov03  successfully accounting for the observed carrier dynamics on picosecond timescales, we describe the system by the RWA Hamiltonian (in the rotating frame)

$$H=H_\mathrm{X}+H_{\mathrm{ph}}+H_{\mathrm{int}}.$$

Here $`H_\mathrm{X}=f(t)(a+a^{})`$ describes the interaction of the ground-state exciton with the resonant laser pulse modulated by the real envelope function $`f(t)`$ ($`a,a^{}`$ are excitonic operators), $`H_{\mathrm{ph}}=_𝒌\omega _𝒌b_𝒌^{}b_𝒌`$ is the lattice Hamiltonian ($`b_𝒌^{},b_𝒌`$ refer to the phonon mode $`𝒌`$ whose energy is $`\omega _𝒌`$), and

$$H_{\mathrm{int}}=a^{}a\underset{𝒌,s}{}\left(f_𝒌b_𝒌^{}+f_𝒌^{}b_𝒌\right)$$

is the interaction Hamiltonian ($`f_𝒌`$ are the coupling constants for the individual phonon modes). Due to momentum conservation, only long wavelength phonon modes (compared to dot size) are effectively coupled jacak02a ; jacak03a . For the timescales relevant here, the damping effect is caused mostly by acoustic phonons vagov03 ; forstner03 . In the calculations, we include deformation potential coupling to longitudinal acoustical (LA) phonons, which describes the modification of the band structure due to lattice compression and is different for electrons and holes. Since no cancellation is involved, its effect depends only weakly on the exact wavefunction geometry. We assume that transitions to higher states, including possible inter-valley transitions, are well separated in energy, compared to the spectral characteristics of the dynamics considered here and may be neglected. Also bi-exciton transitions are assumed to be eliminated by a suitable polarization choice. Piezoelectric coupling to acoustic phonons has been shown to be negligible if the electron and hole function overlapkrummheuer02  and is neglected here (but might contribute if the wavefunction geometry is different). We have verified that the Fröhlich coupling to the longitudinal optical (LO) phonons is negligible, due not only to charge cancelation but mostly to large LO phonon frequencies, compared to the timescales involved.

Initially, the state of the system is described by the density matrix $`\varrho =\varrho _0\varrho _\mathrm{l}(T)`$, where $`\varrho _0`$ is the state with no exciton and $`\varrho _\mathrm{l}(T)`$ is the thermal equilibrium of lattice modes at the temperature $`T`$. The dynamics of the system is calculated by the perturbation expansion of the evolution operator in Born approximation. The lattice degrees of freedom are then traced out leading to the reduced density matrix for the carrier subsystem, which contains all the information accessible by optical methods. The details of the procedure may be found in alicki02a . As a result, the reduced density matrix after time $`t`$ may be written as $`\varrho (t)=U(t)(\varrho _0+\varrho _1)U^{}(t)`$, where $`U(t)`$ is the driven evolution for the unperturbed carrier system and $`\varrho _1`$ describes the perturbation to the carrier dynamics coming from phonons. The matrix elements of the latter, in the basis of empty dot ($`|0`$) and one exciton ($`|1`$) states, may be written as

$`0|\varrho _1|0=1|\varrho _1|1`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega {\displaystyle \frac{R(\omega )}{\omega ^2}}S_{00}(\omega ),`$ (1a)
$`0|\varrho _1|1=1|\varrho _1|0^{}`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega {\displaystyle \frac{R(\omega )}{\omega ^2}}S_{01}(\omega ).`$ (1b)
An additional component describing a unitary correction, e.g. light coupling renormalization and energy shifts, may be canceled by an appropriate modification of $`H_\mathrm{X}`$ (these effects may lead e.g. to intensity dependence of the observed Rabi frequencyforstner03 ). The function $`R(\omega )`$ is the spectral density of the reservoir, fully characterizing the lattice properties at this order of perturbation treatment, defined as

$$R(\omega )=\underset{𝒌}{}|f_𝒌|^2\left[(n_𝒌+1)\delta (\omega \omega _𝒌)+n_𝒌\delta (\omega +\omega _𝒌)\right].$$

It depends on the material parameters and QD size via the coupling constants $`f_𝒌`$ (see e.g. jacak03b  for explicit formulas). The functions $`S_{00}(\omega ),S_{01}(\omega )`$ contain the complete information on the optically controlled carrier dynamics. If we assume that final state is measured after a time long compared to the phonon oscillation periods, they may be written as

$`S_{00}(\omega )`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left[\mathrm{sin}^2\alpha +|K_\mathrm{s}(\omega )|^2\right],`$
$`S_{01}(\omega )`$ $`=`$ $`{\displaystyle \frac{1}{8}}\left[\mathrm{sin}2\alpha +2\mathrm{R}\mathrm{e}K_\mathrm{s}(\omega )K_\mathrm{c}^{}(\omega )\right],`$
where

$`K_\mathrm{s}(\omega )`$ $`=`$ $`{\displaystyle _s^t}𝑑\tau e^{i\omega \tau }{\displaystyle \frac{d}{d\tau }}\mathrm{sin}2F(\tau ),`$ (2a)
$`K_\mathrm{c}(\omega )`$ $`=`$ $`{\displaystyle _s^t}𝑑\tau e^{i\omega \tau }{\displaystyle \frac{d}{d\tau }}\mathrm{cos}2F(\tau ),`$ (2b)
where $`s,t`$ are the initial and final time of the evolution, $`F(t)=_s^t𝑑\tau f(\tau )`$ is the rotation (on the Bloch sphere) performed on the exciton state up to time $`t`$ and $`\alpha `$ is the total angle of rotation. For a given $`\alpha `$, the functions (2a,b) actually depend only on $`\omega \tau _\mathrm{p}`$, where $`\tau _\mathrm{p}`$ is the pulse duration.

In a typical experiment stievater01 ; zrenner02 ; borri02a  one measures the average dot occupation, $`n=1|\varrho (\mathrm{})|1`$, after a pulse of fixed length but variable amplitude. The simulation of such an experiment within the presented model is shown in Fig. 1 for Gaussian pulses \[$`\tau _\mathrm{p}`$ is the full width at half maximum of the pulse envelope $`f(t)`$\]. On the grounds of earlier numerical calculations jacak03b , we have assumed the exciton ground state to be of nearly product character, with Gaussian electron and hole wavefunctions with in-plane localization widths of 4.9 nm and 4.0 nm, respectively and with growth-direction width of 1 nm (this corresponds to a dot of roughly 20 nm diameter and 4 nm height). The oscillations are almost perfect for very short pulses ($`1`$ ps), then loose their quality for longer pulse durations ($`10`$ ps). Although this might be expected from any simple decoherence model, the striking feature is that the effect dramatically grows for higher oscillations, despite the fact that the whole process has exactly constant duration. Even more surprising is the improvement of the quality of oscillations for long pulses ($`50`$ nm) where, in addition, the first oscillation is nearly perfect.

The formulas (1a,b) quantify the idea of resonance between the induced dynamics and lattice modes: in Fig. 2 the phonon spectral density is compared to the non-linear frequency characteristics of the optically controlled exciton dynamics for $`\tau _\mathrm{p}=1`$ ps (the characteristics for other durations is easily obtained by scaling). According to (1a,b), the overlap of these spectral characteristics with the phonon spectral density gives the perturbation of the coherent carrier dynamics.

For growing number of rotations $`n`$, the nonlinear pulse spectrum $`S_{00}(\omega )`$ develops a series of maxima of increasing strength (Fig. 2a,b,c). The position of the last and highest maximum corresponds approximately to $`2\pi n/\tau _\mathrm{p}`$, in accordance with the semiclassical resonance concept. However, spectral components are also present at all the frequencies $`2\pi n^{}/\tau _\mathrm{p}`$, $`n^{}<n`$, which is due to the turning on/off of the pulse. It is interesting to note that for high $`n`$, the low-frequency part of $`S_{00}(\omega )`$ does not evolve with $`n`$ (Fig. 2d). It is now clear that there are two ways of minimizing the overlap: either the pulse must be so short that all the maxima of $`S_{00}(\omega )`$ are pushed to the right into the exponentially vanishing tail of the reservoir spectral density $`R(\omega )`$, or the pulse must be very long, to “squeeze” the spectral function near $`\omega =0`$ and thus reduce its area. In the latter case, the maxima developing with growing number of oscillations will eventually overlap with $`R(\omega )`$ destroying the coherence.

Although it might seem that speeding up the process is the preferred solution, it is clear that this works only because no high frequency features are included into the present model. In reality, speeding up the dynamics is limited e.g. by the presence of excited states and non-adiabatically enhanced LO phonon coupling fomin98 . Moreover, it turns out that the resulting dynamics, even within this model, is actually not fully coherent. It has been shown vagov02a  that superposition of states created by an ultra-short $`\pi /2`$ pulse becomes corrupted, preventing a second pulse (after some delay time $`\mathrm{\Delta }t`$) from generating the final state of $`n=1`$ with unit efficiency. In order to prove the fully coherent character of carrier dynamics it is necessary to demonstrate the stability of the intermediate state in a two-pulse experiment. The simulations of such an experiment are shown in Fig. 1d. A short pulse ($`\tau _\mathrm{p}=0.1`$ ps) creates a superposition of bare states (surrounded by non-distorted lattice) which then decohere due to dressing processes krummheuer02 ; jacak03a . As a result, the exciton cannot be created by the second pulse with unit probability vagov02a . For a longer pulse ($`\tau _\mathrm{p}=1.0`$ ps), the lattice partly manages to follow the evolution of charge distribution during the optical operation and the destructive effect is smaller. Finally, if the carrier dynamics is slow compared to the lattice response times ($`\tau _\mathrm{p}10.0`$ ps), the lattice distortion follows adiabatically the changes in the charge distribution and the superposition created by the first pulse is an eigenstate of the interacting carrier-lattice system, hence does not undergo any decoherence and the final effect is the same for any delay time (its quality limited by decoherence effects during pulsing). In fact, splitting the $`\pi `$–pulse into two corresponds to slowing down the carrier dynamics which, in the absence of decoherence during delay time, improves the quality of the final state, as seen in Fig 1d.

The above analysis shows that damping of pulse-area-dependent Rabi oscillations due to interaction with lattice modes is a fundamental effect of non-Markovian character: it is due to a semiclassical resonance between optically induced confined charge dynamics and lattice modes. The destructive effect may be minimized both by speeding up and slowing down the dynamics. However, in the former case, the system passes through unstable (decohering) states. Moreover, fast operation on a real system induces many undesirable effects: transitions to higher states, bi-exciton generation or resonant LO phonon dynamics. On the other hand, for slow operation, the number of “good” oscillations is limited. Thus, it is impossible to perform an arbitrary number of fully coherent Rabi oscillations on an exciton confined in a quantum dot.

By increasing the pulse duration as $`\tau _\mathrm{p}\alpha ^2`$, the phonon effect on the exciton dynamics may be kept constant. However, in this case the achievable number of oscillations is strongly restricted by the exciton lifetime and other (thermally activated) processes. Eliminating the radiative losses e.g. by using stimulated Raman adiabatic passage instead of a simple optical excitation troiani03  seems to be a promising direction from this point of view.

The model presented above accounts for the decrease of the quality of Rabi oscillations in the short duration range observed in the experiment borri02a . It predicts, however, that this trend is reversed for longer pulse durations. The quantitative value of 96% for the first maximum of the oscillations with $`\tau _\mathrm{p}=1`$ ps agrees very well with the experimental result zrenner02  although the following extrema are much worse in reality than predicted here. This suggests an increased lattice response at higher frequencies which may be due to more complicated wavefunction geometry heitz99  or electric-field induced charge separation leading to strong piezoelectric effects krummheuer02 . Nevertheless, as far as it may be inferred from the experiment zrenner02 , the decrease of the oscillation quality seems to saturate after one full Rabi rotation, as predicted by the model calculations.

In conclusion, we have studied the carrier-lattice dynamics for optically induced Rabi oscillations of exciton occupation in a quantum dot. We have shown that the lattice response is resonantly driven by a combination of the (linear) pulse spectrum and the Rabi frequency. This semiclassical interpretation allows one to qualitatively predict the damping effect based on the general knowledge of the spectral properties of the lattice modes and leads to the conclusion that the resonant lattice response precludes performing a large number of fully coherent oscillations in this system.

We are grateful to T. Kuhn for reading the manuscript and many useful comments. P.M. is grateful to M. Horodecki for discussions and to the Alexander von Humboldt Foundation for generous support. Supported by the Polish Ministry of Scientific Research and Information Technology under the (solicited) grant No PBZ-MIN-008/P03/2003 and by the Polish KBN under grant No. PB 2 PO3B 085 25.