text:
# Commutation Technique for interacting close-to-boson excitons

## 1 Exact semiconductor hamiltonian and one-exciton states

The semiconductor conduction electrons correspond to S states with a $`s=\pm 1/2`$ spin along an arbitrary $`z`$ direction. Let us call $`a_{𝐤,s}^{}`$ the creation operator for such a free conduction electron with momentum $`𝐤`$ and spin $`s`$. For the free valence electrons, the situation is much more complex : According to the Kohn-Luttinger representation (12), the valence energy matrix (13) is given by

$$H_𝐤=Ak^2+B(𝐤.𝐈_{3/2})^2,$$
(1)
if one neglects warping, $`𝐈_{3/2}`$ being the $`3/2`$ angular momentum matrix. The $`H_𝐤`$ diagonalization, which is done by taking the $`z`$ axis along the $`𝐤`$ direction, generates the so-called light and heavy valence electrons. If we now add Coulomb interaction, we have shown (14) that this Coulomb interaction is not diagonal between heavy and light valence electrons, a fact which seriously complicates the problem since it mixes the $`m=\pm 3/2`$ with the $`m=\pm 1/2`$ valence states. In quantum wells, one can forget about this complexity because the $`\pm 3/2`$ and $`\pm 1/2`$ valence states having different hole masses, they are shifted differently by the confinement, so that the heavy-light hole Coulomb coupling usually gives a negligible contribution due to the splitting energy induced by the confinement, which appears in denominators. In bulk material however, the only consistent way to forget about these heavy-light hole Coulomb couplings is to assume $`B=0`$, which corresponds to take all the valence electrons with the same mass. In doing so, $`H_𝐤`$ is diagonal for any $`𝐤`$ direction, so that the valence electron states are simply characterized by a quantum number $`m=\pm 3/2,\pm 1/2`$, along a direction which is now arbitrary. Let us call $`b_{𝐤,m}^{}`$ the creation operator for such a free hole with momentum $`𝐤`$ and angular momentum $`m`$.

The exact electron-hole semiconductor hamiltonian then reads

$$H=H_0+V_{\mathrm{eh}}+V_{\mathrm{ee}}+V_{\mathrm{hh}},$$
(2)
where

$$H_0=\underset{𝐤,s}{}(\mathrm{\Delta }+ϵ_𝐤^e)a_{𝐤,s}^{}a_{𝐤,s}+\underset{𝐤,m}{}ϵ_𝐤^hb_{𝐤,m}^{}b_{𝐤,m},ϵ_𝐤^e=\frac{\mathrm{}^2k^2}{2m_e}ϵ_𝐤^h=\frac{\mathrm{}^2k^2}{2m_h},$$
(3)
$`\mathrm{\Delta }`$ being the band gap. For quantum wells, $`m=\pm 3/2`$, while for bulk materials $`m=\pm 3/2,\pm 1/2`$, as a unique hole mass $`m_h`$ is then assumed for the two hole bands. Within this approximation, the Coulomb scattering of a hole is diagonal with respect to its “spin” $`m`$, so that a hole $`(𝐤,m)`$ is scattered into $`(𝐤+𝐪,m)`$ with the *same* $`m`$. The hole-hole, electron-hole and electron-electron interactions thus read

$`V_{\mathrm{hh}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐪\mathrm{𝟎},𝐤,𝐤^{},m,m^{}}{}}V_𝐪b_{𝐤+𝐪,m}^{}b_{𝐤^{}𝐪,m^{}}^{}b_{𝐤^{},m^{}}b_{𝐤,m},`$ (4)
$`V_{\mathrm{eh}}`$ $`=`$ $`{\displaystyle \underset{𝐪\mathrm{𝟎},𝐤,𝐤^{},s,m}{}}V_𝐪a_{𝐤+𝐪,s}^{}b_{𝐤^{}𝐪,m}^{}b_{𝐤^{},m}a_{𝐤,s},`$ (5)
$`V_{\mathrm{ee}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐪\mathrm{𝟎},𝐤,𝐤^{},s,s^{}}{}}V_𝐪a_{𝐤+𝐪,s}^{}a_{𝐤^{}𝐪,s^{}}^{}a_{𝐤^{},s^{}}a_{𝐤,s},`$ (6)
where $`V_𝐪=\frac{4\pi e^2}{𝒱q^2}`$ in 3D and $`V_𝐪=\frac{2\pi e^2}{𝒮q}`$ in 2D, with $`𝒱`$ being the sample volume and $`𝒮`$ the well area.

From the free electron and free hole operators $`a_{𝐤,s}^{}`$ and $`b_{𝐤,m}^{}`$, one can construct the exciton operators $`B_i^{}`$ as

$$B_i^{}B_{\nu _i,𝐐_i,s_i,m_i}^{}=\underset{𝐤_i}{}𝐤_i|x_{\nu _i}a_{𝐊_i^e,s_i}^{}b_{𝐊_i^h,m_i}^{},$$
(7)
where we have set

$$𝐊_i^e=𝐤_i+\alpha _e𝐐_i,𝐊_i^h=𝐤_i+\alpha _h𝐐_i,\alpha _{e,h}=\frac{m_{e,h}}{m_e+m_h}.$$
(8)
$`𝐐_i`$ and $`\nu _i`$ are the $`i`$ exciton center of mass momentum and relative motion quantum parameter : In 3D, $`\nu _i`$ corresponds to the quantum numbers $`(n_i,l_i,m_i)`$, while in 2D it corresponds to two quantum numbers only $`(n_i,m_i)`$. $`𝐤|x_{\nu _i}`$ is the $`i`$ exciton relative motion wave function in $`𝐤`$ space. In a similar way, it is possible to write the free electron-hole pair operator in terms of exciton operators :

$$a_{𝐤_e,s_i}^{}b_{𝐤_h,m_i}^{}=\underset{\nu _i}{}x_{\nu _i}|\alpha _h𝐤_e\alpha _e𝐤_hB_{\nu _i,𝐤_e+𝐤_h,s_i,m_i}^{},$$
(9)
as checked by inserting Eq (9) into Eq (7). Let us stress that this sum contains the bound states as well as the diffusive states.

The one-exciton state $`B_i^{}|v`$, where $`|v`$ is the electron-hole vacuum state, is eigenstate of $`H`$. Using Eqs (2-9), and the following Schrödinger equation for the electron-hole relative motion in $`𝐤`$ space,

$$(ϵ_𝐤^e+ϵ_𝐤^h)𝐤|x_{\nu _i}\underset{𝐪\mathrm{𝟎}}{}V_𝐪𝐤+𝐪|x_{\nu _i}=ϵ_{\nu _i}𝐤|x_{\nu _i},$$
(10)
where $`ϵ_{\nu _i}`$ is the energy of the relative motion state $`|x_{\nu _i}`$, it is indeed easy to check that

$$HB_i^{}|v=(H_0+V_{\mathrm{eh}})B_i^{}|v=E_iB_i^{}|v,E_i=\mathrm{\Delta }+ϵ_{\nu _i}+\frac{\mathrm{}^2𝐐_i^2}{2(m_e+m_h)}.$$
(11)
Let us again stress that in 3D this result strongly relies on the fact that the Coulomb interaction is diagonal between holes, i. e. it does not contain terms in $`b_{𝐤^{}𝐪,m^{}}^{}b_{𝐤^{},m}`$ with $`m^{}m`$.

In many problems dealing with excitons, the photons play an important rôle. As a $`\sigma _+`$ photon with spin $`(J=1,M=\pm 1)`$ creates an exciton with the same *total* angular momentum, it may appear appropriate to use, instead of the $`(s_i,m_i)`$ angular momentum variables of the separate electron and hole, the $`(J_i,M_i)`$ total angular momentum variables of the exciton, with $`(J_i=2,M_i=\pm 2,\pm 1,0)`$ or $`(J_i=1,M_i=\pm 1,0)`$. One can easily go from one set of operators to the other by

$`B_{\nu _i,𝐐_i,J_i,M_i}^{}`$ $`=`$ $`{\displaystyle \underset{s_i,m_i}{}}s_i,m_i|J_i,M_iB_{\nu _i,𝐐_i,s_i,m_i}^{},`$ (12)
$`B_{\nu _i,𝐐_i,s_i,m_i}^{}`$ $`=`$ $`{\displaystyle \underset{J_i,M_i}{}}J_i,M_i|s_i,m_iB_{\nu _i,𝐐_i,J_i,M_i}^{}.`$ (13)
Only two of these $`B_{\nu _i,𝐐_i,J_i,M_i}^{}`$ states are coupled to light, namely $`(J_i=1,M_i=\pm 1)`$, while the six other exciton states correspond to the so-called “dark” excitons. However, even if these dark excitons are not coupled to light, they are generated by Coulomb scatterings so that we must keep these 8 operators $`B_{\nu _i,𝐐_i,J_i,M_i}^{}`$ anyway in order to possibly describe the exciton-exciton scatterings properly. Moreover, as the Coulomb interaction is diagonal within the $`(s,m)`$ quantum numbers, but not within the $`(J,M)`$’s, it is in fact far simpler to work with the 8 operators $`B_{\nu ,𝐐,s,m}^{}`$ for all processes dealing with Coulomb interactions, and just at the beginning and the end of the calculations, use Eqs (12-13) to transform these $`B_{\nu ,𝐐,s,m}^{}`$ operators into the $`B_{\nu ,𝐐,J,M}^{}`$ operators coupled to the light. This is why all over this work, we will use the $`B_{\nu ,𝐐,s,m}^{}`$ operators only.

## 2 Coulomb creation operator $`V_i^{}`$ and direct Coulomb scattering $`\xi _{lnij}^{\mathrm{dir}}`$

In standard problems with interactions, one usually divides the system hamiltonian $`H`$ into $`H_0+V`$, where $`H_0`$ is the so-called “non-interacting” part, i. e. the part of $`H`$ which can be diagonalized, while $`V`$ is the interaction part, i. e. the part of $`H`$ which cannot be handled exactly, but is hoped to be small enough to be treated as a perturbation. In addition, when $`H`$ can be written as $`H_0+V`$, the $`H_0`$ eigenstates form an orthogonal basis which can be used to expand any state of the system.

For interacting excitons, even if we could guess that the energies $`E_i`$ of individual excitons should appear in the “non-interacting” part of an hypothetical exciton hamiltonian, one cannot divide the electron-hole Coulomb interaction given in Eq (5), into a part which would bind a specific electron to a specific hole to form the $`i`$ exciton, and a “rest” which would make this exciton to interact with other excitons : The electrons, as well as the holes, being indistinguishable particles, such a formal separation is indeed impossible.

A separation, similar in its spirit to the splitting $`H=H_0+V`$, is nevertheless possible through our commutation technique. If $`H`$ is the *exact* semiconductor hamiltonian, given in Eqs (2-6), and if $`B_i^{}`$ is the *exact* exciton creation operator, given in Eq (7), we find that their commutator reads

$$[H,B_i^{}]=E_iB_i^{}+V_i^{},$$
(14)
where $`V_i^{}`$ is given by

$`V_i^{}={\displaystyle \underset{𝐪\mathrm{𝟎},\nu _l}{}}V_𝐪\gamma _{li}(𝐪)B_{\nu _l,𝐐_i+𝐪,s_i,m_i}^{}W_𝐪,`$
$`W_𝐪={\displaystyle \underset{𝐩,s}{}}a_{𝐩𝐪,s}^{}a_{𝐩,s}{\displaystyle \underset{𝐩,m}{}}b_{𝐩𝐪,m}^{}b_{𝐩,m}.`$ (15)
The coefficient $`\gamma _{li}(𝐪)`$, characterizes the scattering of a $`\nu _i`$ exciton into a $`\nu _l`$ state under a $`𝐪`$ Coulomb excitation. Using the following relation,

$$\underset{𝐤}{}x_\nu |𝐤+\alpha 𝐪𝐤|x_\nu ^{}=x_\nu |e^{i\alpha 𝐪.𝐫}|x_\nu ^{},$$
(16)
it is easy to show that it is given by

$$\gamma _{li}(𝐪)=x_{\nu _l}|e^{i\alpha _h𝐪.𝐫}e^{i\alpha _e𝐪.𝐫}|x_{\nu _i}.$$
(17)
The analogy with the usual separation $`H=H_0+V`$ is quite transparent if we note that Eq (14) also reads

$$HB_i^{}=B_i^{}(H+E_i)+V_i^{}.$$
(18)
By considering the state $`B_i^{}|\varphi `$, where $`|\varphi `$ is any electron-hole state, Eq (18) leads to

$$HB_i^{}|\varphi =E_iB_i^{}|\varphi +B_i^{}H|\varphi +V_i^{}|\varphi .$$

In the first term, the contribution of the $`i`$ exciton to the energy of the system is just $`E_i`$ as if this exciton were not interacting with the other electrons or holes of $`|\varphi `$. The second term corresponds to $`H`$ acting on $`|\varphi `$ *independently* of the presence of the $`i`$ exciton : The operator $`(H+E_i)`$, on the right hand side of $`B_i^{}`$, thus plays the rôle of the $`H_0`$ part of the hamiltonian for usual problems in which $`H`$ can be written as $`H_0+V`$. The third term is there because the $`i`$ exciton does in fact interact with $`|\varphi `$. This operator $`V_i^{}`$ thus describes all possible Coulomb interactions between the $`i`$ exciton and the rest of the system. It has to be seen as the formal equivalent of $`V`$ in $`H=H_0+V`$. It is however important to stress that, while the usual $`V`$’s conserve particles, i. e. contain the same number of creation operators $`a^{}`$ and destruction operators $`a`$, the operator $`V_i^{}`$ is not a real potential in the sense that it contains one additional electron-hole pair creation operator $`B^{}`$ (see Eq (15)). This is why we call it Coulomb *creation* operator. In addition, usual two-body potentials have prefactors which depend on four indices which are the ones of two initial and two final states. We here see that the prefactors appearing in $`V_i^{}`$ depend on two indices only, $`l`$ and $`i`$.

In order to cope with these difficulties and generate a “scattering” which depends on four indices, we can push the commutation technique one step further and calculate $`[V_i^{},B_j^{}]`$. Using Eqs (7-9) and (15), we find

$$[V_i^{},B_j^{}]=\underset{𝐪\mathrm{𝟎},\nu _l,\nu _n}{}V_𝐪\gamma _{li}(𝐪)\gamma _{nj}(𝐪)B_{\nu _l,𝐐_i+𝐪,s_i,m_i}^{}B_{\nu _n,𝐐_j𝐪,s_j,m_j}^{},$$
(19)
which formally reads

$$[V_i^{},B_j^{}]=\underset{l,n}{}\xi _{lnij}^{\mathrm{dir}}B_l^{}B_n^{}.$$
(20)
If we compare Eq (19) to Eq (20), and symmetrize the result with respect to $`(l,n)`$ (which will appear convenient afterwards), we find that the direct Coulomb scattering $`\xi _{lnij}^{\mathrm{dir}}`$ can be written as

$$\xi _{lnij}^{\mathrm{dir}}=\xi _{nlij}^{\mathrm{dir}}=\xi _{lnji}^{\mathrm{dir}}=\frac{1}{2}[\mathrm{\Delta }^{\mathrm{dir}}({}_{nj}{}^{li})\widehat{\xi }^{\mathrm{dir}}({}_{nj}{}^{li})+(ln)].$$
(21)
It contains an angular momentum contribution and an orbital contribution which appear as two independent factors. For direct processes, this angular momentum contribution,

$$\mathrm{\Delta }^{\mathrm{dir}}({}_{nj}{}^{li})=\delta _{s_l,s_i}\delta _{m_l,m_i}\delta _{s_n,s_j}\delta _{m_n,m_j},$$
(22)
just says that the $`l`$ exciton has the same electron and hole momenta as the $`i`$ exciton and similarly for the $`n`$ and $`j`$ excitons. Let us stress that this angular momentum part would be much more complicated if the $`(J,M)`$ variables for the exciton total angular momentum were used ; this is why it is indeed appropriate to keep these $`(s,m)`$ variables as long as we deal with Coulomb scatterings even if they are not the good variables for semiconductor-photon interaction.

The orbital part of this direct Coulomb scattering is given by

$$\widehat{\xi }^{\mathrm{dir}}({}_{nj}{}^{li})=\delta _{𝐐_l+𝐐_n,𝐐_i+𝐐_j}V_{𝐐_l𝐐_i}\gamma _{li}(𝐐_l𝐐_i)\gamma _{nj}(𝐐_n𝐐_j).$$
(23)
It of course contains the fact that the center of mass momenta have to be conserved, $`𝐐_i+𝐐_j=𝐐_l+𝐐_n`$, in the scattering. It also contains the two factors which characterize the scatterings of one exciton from a $`\nu _i`$ state to a $`\nu _l`$ state and the other exciton from a $`\nu _j`$ state to a $`\nu _n`$ state, under the Coulomb excitation $`𝐐_l𝐐_i=(𝐐_n𝐐_j)`$.

The appendix A contains the explicit calculation of these $`\gamma _{li}(𝐪)`$ factors for the lowest S states.

Even if the above expression of $`\widehat{\xi }^{\mathrm{dir}}({}_{nj}{}^{li})`$ is quite transparent and extremely convenient for explicit calculations, it will appear useful, when we will generate exchange Coulomb processes, to note that this direct Coulomb scattering is also equal to

$`\widehat{\xi }^{\mathrm{dir}}({}_{nj}{}^{li})={\displaystyle \underset{𝐤_l,𝐤_n,𝐤_i,𝐤_j}{}}x_{\nu _l}|𝐤_lx_{\nu _n}|𝐤_n𝐤_i|x_{\nu _i}𝐤_j|x_{\nu _j}`$
$`\times {\displaystyle \underset{𝐪\mathrm{𝟎}}{}}V_𝐪[\delta _{𝐊_l^e,𝐊_i^e+𝐪}\delta _{𝐊_l^h,𝐊_i^h}(eh)][\delta _{𝐊_n^e,𝐊_j^e𝐪}\delta _{𝐊_n^h,𝐊_j^h}(eh)].`$ (24)
In appendix B, we show that this $`\widehat{\xi }^{\mathrm{dir}}({}_{nj}{}^{li})`$ is nothing but the *direct* part of the Coulomb interaction widely quoted (7,15) in the effective boson exciton hamiltonian, namely

$`\widehat{\xi }^{\mathrm{dir}}({}_{nj}{}^{li})={\displaystyle }d𝐫_ed𝐫_hd𝐫_e^{}d𝐫_h^{}\varphi _l^{}(𝐫_e,𝐫_h)\varphi _n^{}(𝐫_e^{},𝐫_h^{})`$
$`\times [V_{ee^{}}+V_{hh^{}}V_{eh^{}}V_{e^{}h}]\varphi _i(𝐫_e,𝐫_h)\varphi _j(𝐫_e^{},𝐫_h^{})=\left[\widehat{\xi }^{\mathrm{dir}}({}_{jn}{}^{il})\right]^{},`$ (25)
where $`\varphi _i(𝐫_e,𝐫_h)`$ is the whole wave function of the $`i=(\nu _i,𝐐_i)`$ exciton,

$$𝐫_e,𝐫_h|B_i^{}|v=\varphi _i(𝐫_e,𝐫_h)=\frac{1}{\sqrt{𝒱}}e^{i𝐐_i.(\alpha _e𝐫_e+\alpha _h𝐫_h)}𝐫_e𝐫_h|x_{\nu _i}.$$
(26)
By these two successive commutators $`[H,B_i^{}]`$ and $`[V_i^{},B_j^{}]`$, we have found a formal way to generate a scattering of two excitons $`(i,j)`$ into two other excitons $`(l,n)`$. If we look at its expression given in Eq (25), we see that the coefficient $`\widehat{\xi }^{\mathrm{dir}}({}_{nj}{}^{li})`$ obtained by this procedure turns out to have a very simple physical meaning : It corresponds to all electron-electron, hole-hole and electron-hole Coulomb interactions *between* $`(i,j)`$ and $`(l,n)`$, when the initial and final excitons are made with the same electron-hole pairs $`(e,h)`$ and $`(e^{},h^{})`$.

The scattering processes corresponding to $`\widehat{\xi }^{\mathrm{dir}}({}_{nj}{}^{li})`$ are shown in Fig. (1) : One of the electrons, or holes, suffers a $`(𝐐_l𝐐_i)`$ Coulomb excitation, while another particle has a $`(𝐐_l𝐐_i)`$ excitation, the excitons before and after the scattering being made with the *same* electron-hole pairs. This is why we call it a *direct* Coulomb scattering, by contrast with the exchange Coulomb scatterings which will appear below, in which the excitons before and after scattering are made with different pairs.

## 3 Boson departure operator $`D_{ij}`$ and bare exchange coefficient $`\lambda _{lnij}`$

Even if the above procedure allows to *formally* extract from the whole electron-hole interaction the part which corresponds to a Coulomb interaction *between* excitons, such a formal extraction seeming not obvious at first, we are far from having picked all the physics which controls the interactions between excitons. Another rather subtle origin of these interactions comes from the fermionic character of the excitons. From handwaving arguments, we can say that two excitons feel each other *even in the absence of any Coulomb interaction*, because they are made of two fermions and the fermions of two excitons must be in different states. This condition in itself produces a “link” between excitons. There is no need of Coulomb forces. It is thus very likely that the exciton-exciton interaction has to contain a purely fermionic contribution, quite different from the Coulomb interaction dressed by exchange processes as thought by everyone up to now.

In this quite tricky determination of the correct exciton-exciton interaction, it is however highly necessary to formalize the above handwaving argument. We want to find a coefficient, which depends on four exciton indices $`(lnij)`$, and which originates from the fermionic character of the excitons only. This coefficient has of course to be linked to the fact that the excitons are not real bosons. The most direct way to have this property appearing, is to start with the commutator of two (exact) exciton creation operators. It leads to the boson departure operator $`D_{ij}`$ defined as

$$D_{ij}=\delta _{ij}[B_i,B_j^{}].$$
(27)
Using Eq (7) and the standard commutation rules for the fermion operators $`a^{}`$ and $`b^{}`$, we find that this operator reads

$$D_{ij}=\underset{𝐤_i,𝐤_j}{}x_{\nu _i}|𝐤_i𝐤_j|x_{\nu _j}\left[\delta _{m_i,m_j}\delta _{𝐊_i^h,𝐊_j^h}a_{𝐊_j^e,s_j}^{}a_{𝐊_i^e,s_i}+\delta _{s_i,s_j}\delta _{𝐊_i^e,𝐊_j^e}b_{𝐊_j^h,m_j}^{}b_{𝐊_i^h,m_i}\right].$$
(28)
$`D_{ij}0`$ if we have both $`s_is_j`$ and $`m_im_j`$.

From this boson departure operator, we can construct a four index coefficient by taking the commutator of this $`D_{ij}`$ operator with an exciton operator. Using Eqs (7,9,28), we find

$$[D_{li},B_j^{}]=2\underset{n}{}\lambda _{lnij}B_n^{},$$
(29)
where the bare exchange coefficient $`\lambda _{lnij}`$ can here again be split into an angular momentum part and an orbital part :

$$\lambda _{lnij}=\lambda _{nlij}=\lambda _{lnji}=\frac{1}{2}[\mathrm{\Delta }^{\mathrm{exch}}({}_{nj}{}^{li})\widehat{\lambda }({}_{nj}{}^{li})+(ln)].$$
(30)
The angular momentum part,

$$\mathrm{\Delta }^{\mathrm{exch}}({}_{nj}{}^{li})=\delta _{s_l,s_i}\delta _{m_l,m_j}\delta _{s_n,s_j}\delta _{m_n,m_i},$$
(31)
just says that the spins of the $`l`$ and $`i`$ electrons are the same while the angular momentum of the $`l`$ hole is not the one of the $`i`$ hole, as in $`\mathrm{\Delta }^{\mathrm{dir}}({}_{nj}{}^{li})`$, but the one of the $`j`$ hole. The orbital part appears as

$$\widehat{\lambda }({}_{nj}{}^{li})=\underset{𝐤_l,𝐤_n,𝐤_i,𝐤_j}{}x_{\nu _l}|𝐤_lx_{\nu _n}|𝐤_n𝐤_i|x_{\nu _i}𝐤_j|x_{\nu _j}\delta _{𝐊_l^e,𝐊_i^e}\delta _{𝐊_l^h,𝐊_j^h}\delta _{𝐊_n^e,𝐊_j^e}\delta _{𝐊_n^h,𝐊_i^h}.$$
(32)
In appendix B, we show that this orbital part $`\widehat{\lambda }({}_{nj}{}^{li})`$ can be rewritten in $`𝐫`$ space as

$$\widehat{\lambda }({}_{nj}{}^{li})=d𝐫_ed𝐫_hd𝐫_e^{}d𝐫_h^{}\varphi _l^{}(𝐫_e,𝐫_h)\varphi _n^{}(𝐫_e^{},𝐫_h^{})\varphi _i(𝐫_e,𝐫_h^{})\varphi _j(𝐫_e^{},𝐫_h)=\left[\widehat{\lambda }({}_{jn}{}^{il})\right]^{}.$$
(33)
From the above expression, we see that, in this orbital part, the electrons and the holes forming the excitons are exchanged, the $`l`$ and $`i`$ excitons being made with the same electron but different holes. Note that the $`(ln)`$ term of Eq (30) restores the $`(e,h)`$ symmetry in $`\lambda _{lnij}`$, which is broken in $`\widehat{\lambda }({}_{nj}{}^{li})`$.

Eq (33) shows in a transparent way that this $`\widehat{\lambda }({}_{nj}{}^{li})`$ exchange coefficient can exist because of the composite nature of the excitons which can be formed in different ways, the electrons and the holes being indistinguishable fermions. An equally transparent link between this coefficient and the composite nature of the excitons can be obtained by considering the two-exciton operator $`B_i^{}B_j^{}`$, each $`B^{}`$ being given by Eq (7), and by binding the electrons and holes of these two $`B^{}`$’s in a different way through Eq (9). This gives

$$B_i^{}B_j^{}=\underset{l,n}{}\lambda _{lnij}B_l^{}B_n^{},$$
(34)
with the same prefactor $`\lambda _{lnij}`$ as the one of Eq (29).

The electron-hole exchange appearing in the first term of Eq (30) is shown in Fig. (2). It corresponds to cross the holes when forming the $`(i,j)`$ or $`(l,n)`$ excitons : It is indeed a bare exchange process in the sense that it does not contain any Coulomb interaction, by contrast with the Coulomb exchange processes which will appear later on.

The four $`\delta `$’s appearing in $`\widehat{\lambda }({}_{nj}{}^{li})`$ given in Eq (32) impose the expected exciton total momentum conservation $`𝐐_i+𝐐_j=𝐐_l+𝐐_n`$ as well as $`𝐤_i+𝐤_j=𝐤_l+𝐤_n`$. By expliciting these four $`\delta `$’s, we can rewrite $`\widehat{\lambda }({}_{nj}{}^{li})`$ in a compact form,

$$\widehat{\lambda }({}_{nj}{}^{li})=\delta _{𝐐_l+𝐐_n,𝐐_i+𝐐_j}F_{lnij}(\alpha _e(𝐐_l𝐐_i),\alpha _h(𝐐_n𝐐_i)),$$
(35)
where $`F_{lnij}(𝐩,𝐩^{})`$ is a sum over one momentum only,

$$F_{lnij}(𝐩,𝐩^{})=\underset{𝐤}{}x_{\nu _l}|𝐤\frac{𝐩+𝐩^{}}{2}x_{\nu _n}|𝐤+\frac{𝐩+𝐩^{}}{2}𝐤+\frac{𝐩𝐩^{}}{2}|x_{\nu _i}𝐤\frac{𝐩𝐩^{}}{2}|x_{\nu _j}.$$
(36)
Although less transparent for physical understanding, the above expression of $`\widehat{\lambda }({}_{nj}{}^{li})`$ is more convenient for explicit calculations. When the four excitons are in S states, the function $`F_{lnij}(𝐩,𝐩^{})`$ a priori depends on three independent parameters, $`p`$, $`p^{}`$ and the angle $`\theta `$ between $`𝐩`$ and $`𝐩^{}`$. Some interesting values of these parameters are :

(i) $`p=0`$ or $`p^{}=0`$, which corresponds to $`𝐐_l=𝐐_i`$ or $`𝐐_n=𝐐_i`$ : The final excitons have the same momenta as the initial excitons ;

(ii) $`\theta =0`$ or $`\theta =\pi `$, which corresponds to $`𝐩`$ parallel to $`𝐩^{}`$. $`𝐐_l𝐐_i`$ is parallel to $`𝐐_n𝐐_i𝐐_j𝐐_l`$ in the particular case of $`𝐐_i=𝐐_j`$ or $`𝐐_l=𝐐_n`$, i. e. when the two initial or the two final excitons have the same momentum. In this case, the two independent parameters can be taken either as $`\alpha _e|𝐐_l𝐐_i|`$ and $`\alpha _h|𝐐_n𝐐_i|`$, or better as the momentum transfer $`|𝐐_l𝐐_i|`$ and the mass ratio $`m_e/m_h`$.

In appendix C, we have calculated $`F_{lnij}(𝐩,𝐩^{})`$ when all the $`\nu `$’s are equal to the 1S ground state and $`𝐐_i=𝐐_j=𝐐_l=𝐐_n`$. We get

$`F_{1s1s1s1s}(\mathrm{𝟎},\mathrm{𝟎})`$ $`=`$ $`(33\pi /2)a_x^3/𝒱\mathrm{in}3\mathrm{D}`$ (37)
$`=`$ $`(4\pi /5)a_x^2/𝒮\mathrm{in}2\mathrm{D},`$
$`a_x`$ being the 3D exciton Bohr radius.

We have also numerically computed $`\widehat{\lambda }({}_{nj}{}^{li})`$ in the 2D case, when $`\nu _l=\nu _n=\nu _i=\nu _j=1\mathrm{S}`$, $`𝐐_i=𝐐_j=\mathrm{𝟎}`$, $`𝐐_l=𝐐_n=𝐩`$. Fig. (3) shows

$$\mathrm{\Lambda }(p)=\frac{𝒮}{\lambda _{2D}^2}\widehat{\lambda }({}_{𝐩\mathbf{0}}{}^{𝐩\mathbf{\hspace{0.17em}0}}),$$
(38)
for three values of $`m_e/m_h`$, as a function of $`P=p\lambda _{2D}`$, with $`\lambda _{2D}=a_x/2`$ being the 2D Bohr radius. We see that $`\mathrm{\Lambda }(p)`$ weakly depends on the mass ratio, and that it is significant for values of $`p`$ such that $`p\lambda _{2D}<4`$ only.

## 4 Coulomb exchange scatterings

From Eq (33), we see that the $`\lambda _{lnij}`$ coefficient is dimensionless. So that it has to be “cooked” with quantities homogeneous to an energy in order to possibly appear in an exciton-exciton scattering. From the two energy-like quantities we have yet found in this problem, namely the bare energies of the excitons $`E_i`$ and the direct Coulomb scattering $`\xi _{lnij}^{\mathrm{dir}}`$, there are of course various ways to construct a $`\lambda _{lnij}`$ dependent scattering. In the next section, we will show that two specific combinations of the bare exchange coefficient $`\lambda _{lnij}`$ and the direct Coulomb scattering $`\xi _{lnij}^{\mathrm{dir}}`$ appear in a natural way, namely

$$\xi _{lnij}^{\mathrm{right}}=\underset{p,r}{}\xi _{lnpr}^{\mathrm{dir}}\lambda _{prij},$$
(39)
$$\xi _{lnij}^{\mathrm{left}}=\underset{p,r}{}\lambda _{lnpr}\xi _{prij}^{\mathrm{dir}}.$$
(40)
They are shown in Fig. (4).

The simplest way to calculate these sums is to use Eqs (30-32) for $`\lambda _{lnij}`$ and Eqs (21,22,24) for $`\xi _{lnij}^{\mathrm{dir}}`$. We find

$$\xi _{lnij}^{\mathrm{right}}=\frac{1}{2}(\mathrm{\Delta }^{\mathrm{exch}}({}_{nj}{}^{li})\widehat{\xi }^{\mathrm{right}}({}_{nj}{}^{li})+(ln)),$$
(41)
where the angular momentum conservation part $`\mathrm{\Delta }^{\mathrm{exch}}({}_{nj}{}^{li})`$ is just the one appearing in the bare exchange coefficient $`\lambda _{lnij}`$ (Eq (31)), while the orbital part of this right exchange Coulomb term is given by

$`\widehat{\xi }^{\mathrm{right}}({}_{nj}{}^{li})={\displaystyle \underset{p,r}{}}\widehat{\xi }^{\mathrm{dir}}({}_{nr}{}^{lp})\widehat{\lambda }({}_{rj}{}^{pi})`$
$`={\displaystyle \underset{𝐤_l,𝐤_n,𝐤_i,𝐤_j}{}}x_{\nu _l}|𝐤_lx_{\nu _n}|𝐤_n𝐤_i|x_{\nu _i}𝐤_j|x_{\nu _j}`$
$`\times {\displaystyle \underset{𝐪\mathrm{𝟎}}{}}V_𝐪[\delta _{𝐊_l^e,𝐊_i^e+𝐪}\delta _{𝐊_l^h,𝐊_j^h}\delta _{𝐊_l^h,𝐊_j^h+𝐪}\delta _{𝐊_l^e,𝐊_i^e}][\delta _{𝐊_n^e,𝐊_j^e𝐪}\delta _{𝐊_n^h,𝐊_i^h}\delta _{𝐊_n^h,𝐊_i^h𝐪}\delta _{𝐊_n^e,𝐊_j^e}].`$ (42)
Note that the second term of each bracket does *not* correspond to $`(eh)`$ as for $`\widehat{\xi }^{\mathrm{dir}}({}_{nj}{}^{li})`$ given in Eq (24).

The set of four $`\delta `$’s appearing in $`\widehat{\xi }^{\mathrm{right}}({}_{nj}{}^{li})`$ again impose the expected momentum conservation $`𝐐_i+𝐐_j=𝐐_l+𝐐_n`$. By expliciting these four $`\delta `$ functions, it is possible to rewrite $`\widehat{\xi }^{\mathrm{right}}({}_{nj}{}^{li})`$ in a more compact form in terms of sums over $`𝐪`$ and one $`𝐤`$ only. In addition to the function $`F_{lnij}(𝐩,𝐩^{})`$ appearing already in the bare exchange coefficient (Eq (36)), two other functions enter $`\widehat{\xi }^{\mathrm{right}}({}_{nj}{}^{li})`$. They are defined by

$$G_{lnij}^{(i)}(𝐩,𝐩^{};𝐪)=\underset{𝐤}{}x_{\nu _l}|𝐤\frac{𝐩+𝐩^{}}{2}x_{\nu _n}|𝐤+\frac{𝐩+𝐩^{}}{2}𝐤+\frac{𝐩𝐩^{}}{2}+𝐪|x_{\nu _i}𝐤\frac{𝐩𝐩^{}}{2}|x_{\nu _j},$$
(43)
$$G_{lnij}^{(j)}(𝐩,𝐩^{};𝐪)=\underset{𝐤}{}x_{\nu _l}|𝐤\frac{𝐩+𝐩^{}}{2}x_{\nu _n}|𝐤+\frac{𝐩+𝐩^{}}{2}𝐤+\frac{𝐩𝐩^{}}{2}|x_{\nu _i}𝐤\frac{𝐩𝐩^{}}{2}+𝐪|x_{\nu _j}.$$
(44)
(We can see that $`F_{lnij}(𝐩,𝐩^{})`$ is also either $`G_{lnij}^{(i)}(𝐩,𝐩^{};\mathrm{𝟎})`$ or $`G_{lnij}^{(j)}(𝐩,𝐩^{};\mathrm{𝟎})`$). In terms of these three functions, the right exchange Coulomb term reads

$`\widehat{\xi }^{\mathrm{right}}({}_{nj}{}^{li})=\delta _{𝐐_l+𝐐_n,𝐐_i+𝐐_j}{\displaystyle \underset{𝐪\mathrm{𝟎}}{}}V_𝐪`$
$`\times [F_{lnij}(\alpha _e(𝐐_l𝐐_i)𝐪,\alpha _h(𝐐_n𝐐_i))+F_{lnij}(\alpha _e(𝐐_l𝐐_i),\alpha _h(𝐐_n𝐐_i)+𝐪)`$
$`G_{lnij}^{(i)}(\alpha _e(𝐐_l𝐐_i),\alpha _h(𝐐_n𝐐_i);𝐪)G_{lnij}^{(j)}(\alpha _e(𝐐_l𝐐_i),\alpha _h(𝐐_n𝐐_i);𝐪)].`$ (45)
This $`\widehat{\xi }^{\mathrm{right}}({}_{nj}{}^{li})`$ coefficient can be rewritten in a form much more transparent for the physical understanding, although much less convenient for calculations. In appendix B, we do show that

$`\widehat{\xi }^{\mathrm{right}}({}_{nj}{}^{li})={\displaystyle }d𝐫_ed𝐫_hd𝐫_e^{}d𝐫_h^{}\varphi _l^{}(𝐫_e,𝐫_h)\varphi _n^{}(𝐫_e^{},𝐫_h^{})`$
$`\times \left[V_{ee^{}}+V_{hh^{}}V_{eh^{}}V_{e^{}h}\right]\varphi _i(𝐫_e,𝐫_h^{})\varphi _j(𝐫_e^{},𝐫_h).`$ (46)
The physical meaning of this $`\widehat{\xi }^{\mathrm{right}}({}_{nj}{}^{li})`$ becomes now clear : It contains the electron-electron and hole-hole interactions *between* two excitons when these excitons are made with their electrons and holes coupled in a different way. With respect to the electron-hole interactions, the situation is however more subtle as these interactions are *between* the excitons on one side, but *inside* the excitons of the other side, due to the exchange of the electrons or holes making the excitons of the two sides. This $`\widehat{\xi }^{\mathrm{right}}({}_{nj}{}^{li})`$ turns out to be nothing but the exchange Coulomb interaction which appears in the exciton-exciton scattering coefficient of the effective boson exciton hamiltonian widely quoted (7).

The coefficient $`\xi _{lnij}^{\mathrm{left}}`$ reads as $`\xi _{lnij}^{\mathrm{right}}`$ except that $`V_{eh^{}}+V_{e^{}h}`$ is replaced by $`V_{eh}+V_{e^{}h^{}}`$ in Eq (46), so that the electron-hole Coulomb interactions are now *between* the $`(i,j)`$ excitons and *inside* the $`(l,n)`$ excitons. As a consequence, the two coefficients $`\xi ^{\mathrm{left}}`$ and $`\xi ^{\mathrm{right}}`$ verify

$$\xi _{lnij}^{\mathrm{left}}=\left[\xi _{ijln}^{\mathrm{right}}\right]^{},$$
(47)
with a similar relation between the orbital parts $`\widehat{\xi }^{\mathrm{left}}`$ and $`\widehat{\xi }^{\mathrm{right}}`$. This can also be seen from Eqs (25,33,39,40).

From $`\xi _{lnij}^{\mathrm{right}}`$ and $`\xi _{lnij}^{\mathrm{left}}`$, it is possible to construct a quantity symmetrical with respect to the electron-hole interactions, by taking

$$\xi _{lnij}^{\mathrm{exch}}=\frac{1}{2}\left(\xi _{lnij}^{\mathrm{right}}+\xi _{lnij}^{\mathrm{left}}\right).$$
(48)
Using Eq (46), the corresponding orbital part of this Coulomb exchange coefficient reads

$`\widehat{\xi }^{\mathrm{exch}}({}_{nj}{}^{li})={\displaystyle }d𝐫_ed𝐫_hd𝐫_e^{}d𝐫_h^{}\varphi _l^{}(𝐫_e,𝐫_h)\varphi _n^{}(𝐫_e^{},𝐫_h^{})`$
$`\times [V_{ee^{}}+V_{hh^{}}{\displaystyle \frac{1}{2}}(V_{eh}+V_{e^{}h^{}}+V_{eh^{}}+V_{e^{}h}]\varphi _i(𝐫_e,𝐫_h^{})\varphi _j(𝐫_e^{},𝐫_h).`$ (49)
It contains the four possible electron-hole Coulomb interactions between two electrons and two holes. This $`\xi _{lnij}^{\mathrm{exch}}`$ could be a reasonable scattering for the interaction part of an hypothetical effective exciton hamiltonian : As $`\xi _{lnij}^{\mathrm{exch}}=\left(\xi _{ijln}^{\mathrm{exch}}\right)^{}`$, its contribution would be hermitian, which is not the case neither for $`\xi _{lnij}^{\mathrm{right}}`$ nor for $`\xi _{lnij}^{\mathrm{left}}`$ alone.

## 5 Matrix element of $`H`$ between two-exciton states

In some physical problems such as the one considered in section VI of this paper, we are led to consider the matrix elements of the exact hamiltonian $`H`$ between two-exciton states, namely $`v|B_lB_nHB_i^{}B_j^{}|v`$. Using Eqs (14,20) and $`V_i^{}|v=0`$, as obvious from Eq (15), we find

$$HB_i^{}B_j^{}|v=(E_i+E_j)B_i^{}B_j^{}|v+\underset{p,r}{}\xi _{prij}^{\mathrm{dir}}B_p^{}B_r^{}|v.$$
(50)
So that this matrix element reads

$$v|B_lB_nHB_i^{}B_j^{}|v=(E_i+E_j)v|B_lB_nB_i^{}B_j^{}|v+\underset{p,r}{}\xi _{prij}^{\mathrm{dir}}v|B_lB_nB_p^{}B_r^{}|v.$$
(51)
To go further, we must calculate $`v|B_lB_nB_i^{}B_j^{}|v`$. This is easily done using Eqs (27) and (29). Since $`D_{ij}|v=0`$, as obvious from Eq (28), we get

$$v|B_lB_nB_i^{}B_j^{}|v=\delta _{li}\delta _{nj}+\delta _{ni}\delta _{lj}2\lambda _{lnij}.$$
(52)
The two first terms are naïve : They differ from zero for $`(l,n)=(i,j)`$, i. e. when the excitons on both sides are identical. The last term is more subtle. It has a fermionic origin : In addition to the fact that the coefficient $`\lambda _{lnij}`$ is directly related to the boson deviation operator $`D_{li}`$ through Eq (29), its fermionic origin can also be traced back to Eq (34) which says that a product of two exciton operators writes as a sum of products of any two other $`B^{}`$’s so that $`B_i^{}B_j^{}`$ in fact contains a “piece” of any $`B_l^{}B_n^{}`$, even if $`(l,n)(i,j)`$. In other words, the $`B_i^{}B_j^{}|v`$ states do not form an orthogonal basis for two-pair states.

By inserting Eq (52) into Eq (51), we get

$$v|B_lB_nHB_i^{}B_j^{}|v=(E_i+E_j)(\delta _{li}\delta _{nj}+\delta _{lj}\delta _{ni}2\lambda _{lnij})+2(\xi _{lnij}^{\mathrm{dir}}\xi _{lnij}^{\mathrm{left}}).$$
(53)
This calculation thus produces the Coulomb exchange term $`\xi _{lnij}^{\mathrm{left}}`$. A similar calculation done with $`H`$ acting on the left gives

$$v|B_lB_nHB_i^{}B_j^{}|v=(E_l+E_n)(\delta _{li}\delta _{nj}+\delta _{lj}\delta _{ni}2\lambda _{lnij})+2(\xi _{lnij}^{\mathrm{dir}}\xi _{lnij}^{\mathrm{right}}).$$
(54)
Here appears $`\xi _{lnij}^{\mathrm{right}}`$. From these two expressions of the $`H`$ matrix element, we deduce that $`\xi _{lnij}^{\mathrm{right}}`$ and $`\xi _{lnij}^{\mathrm{left}}`$ are linked by

$$\xi _{lnij}^{\mathrm{left}}\xi _{lnij}^{\mathrm{right}}=(E_l+E_nE_iE_j)\lambda _{lnij},$$
(55)
so that they are equal for $`E_l+E_n=E_i+E_j`$ only.

It is possible to check that Eq (55) remains valid when $`\xi `$ and $`\lambda `$ are replaced by their orbital parts $`\widehat{\xi }`$ and $`\widehat{\lambda }`$. Using this equation, we can easily calculate the difference $`\widehat{\xi }^{\mathrm{left}}\widehat{\xi }^{\mathrm{right}}`$ and compare it to $`\widehat{\xi }^{\mathrm{left}}`$ (or $`\widehat{\xi }^{\mathrm{right}}`$). Ciuti et al (11) and Rochat et al (9) have calculated $`\widehat{\xi }^{\mathrm{left}}({}_{nj}{}^{li})`$ for 2D scatterings between $`\nu _i=\nu _j=1\mathrm{S}`$, $`𝐐_i=𝐐_j=\mathrm{𝟎}`$, and $`\nu _l=\nu _n=1\mathrm{S}`$, $`𝐐_l=𝐐_n=𝐩`$. Fig. (4) of reference (9) precisely shows

$$g^{\mathrm{exch}}(p)=\frac{\pi ^2𝒮}{4e^2\lambda _{2D}}\widehat{\xi }^{\mathrm{left}}({}_{𝐩\mathbf{0}}{}^{𝐩\mathbf{\hspace{0.17em}0}}),$$
(56)
as a function of $`P=p\lambda _{2D}`$, with $`\lambda _{2D}`$ being the 2D Bohr radius. On our Fig. (5) is shown

$$\mathrm{\Delta }g(p)=\frac{\pi ^2𝒮}{4e^2\lambda _{2D}}[\widehat{\xi }^{\mathrm{left}}({}_{𝐩\mathbf{0}}{}^{𝐩\mathbf{\hspace{0.17em}0}})\widehat{\xi }^{\mathrm{right}}({}_{𝐩\mathbf{0}}{}^{𝐩\mathbf{\hspace{0.17em}0}})]=\frac{\pi ^2}{2}\alpha _e\alpha _hp^2\lambda _{2D}^2\mathrm{\Lambda }(p),$$
(57)
with $`\mathrm{\Lambda }(p)`$ given by Eq (38). We see that, except for very small $`p\lambda _{2D}`$, this difference is of the order of $`g^{\mathrm{exch}}(p)`$, unless $`m_e<<m_h`$.

We can rewrite the $`H`$ matrix element between two-exciton states in a more symmetrical form with respect to $`(l,n)`$ and $`(i,j)`$ by taking half the sum of Eqs (53) and (54). If in addition, we consider this $`H`$ matrix element between two-exciton states *normalized*, but still not orthogonal, namely

$$|\psi _{ji}=|\psi _{ij}=\frac{B_i^{}B_j^{}|v}{v|B_iB_jB_j^{}B_i^{}|v^{1/2}}=\frac{1}{𝒩_{ij}}B_i^{}B_j^{}|v,$$
(58)
where $`𝒩_{ij}=(1+\delta _{ij}2\lambda _{ijij})^{1/2}`$ due to Eq (52), we find, using Eqs (53-54),

$$\psi _{ln}|H|\psi _{ij}=\frac{1}{2}(E_l+E_n+E_i+E_j)\delta _{lnij}+\frac{2}{𝒩_{ln}𝒩_{ij}}\left[\xi _{lnij}^{\mathrm{dir}}\xi _{lnij}^{\mathrm{exch}}(1\delta _{lnij})\eta _{lnij}\right].$$
(59)
$`\delta _{lnij}`$ is equal to 1 when the excitons on both sides are identical i. e. $`(l,n)=(i,j)`$, and 0 otherwise, while $`\eta _{lnij}`$ is given by

$$\eta _{lnij}=\frac{1}{2}(E_l+E_n+E_i+E_j)\lambda _{lnij}.$$
(60)
The first term of this $`H`$ matrix element just corresponds to the energies of the *non-interacting* $`(i,j)`$ or $`(l,n)`$ excitons if, *in addition*, the two-exciton states $`|\psi _{ij}`$ are assumed to be orthogonal. The last terms of this $`H`$ matrix element come from interactions. They have three origins. The first term, $`\xi _{lnij}^{\mathrm{dir}}`$, is a direct Coulomb scattering between excitons made on both sides with the same electrons and holes. The second term, $`\xi _{lnij}^{\mathrm{exch}}`$, is an exchange Coulomb scattering in which the holes (or electrons) making the excitons are exchanged. The last term, $`\eta _{lnij}`$, has a purely fermionic origin : It is directly linked to the fact that the $`|\psi _{ij}`$ and $`|\psi _{ln}`$ states are not orthogonal for $`(l,n)(i,j)`$. It does not exist in diagonal scatterings. Let us however stress that, as $`\eta _{lnij}`$ depends on the *sum* of the four exciton energies, it is band gap dependent ; so that it is very unlikely that such a $`\eta _{lnij}`$ term appears in a physical scattering. We are going to come back to this problem in the next paragraph.

It might be useful to mention that the diagonal matrix element of $`H`$ reduces to

$$\psi _{ij}|H|\psi _{ij}=E_i+E_j+(\xi _{ijij}^{\mathrm{dir}}\xi _{ijij}^{\mathrm{exch}})/(1+\delta _{ij}2\lambda _{ijij}),$$
(61)
with $`\xi _{ijij}^{\mathrm{exch}}`$ possibly replaced by $`\xi _{ijij}^{\mathrm{right}}`$ or $`\xi _{ijij}^{\mathrm{left}}`$ as these three scatterings are equal for diagonal processes. $`(E_i+E_j)`$ is just what we would expect for the expectation value of the hamiltonian between two *non-interacting* boson-excitons $`i`$ and $`j`$. The other terms come from scattering processes between the $`i`$ and $`j`$ excitons resulting from Coulomb interaction and Pauli exclusion.

## 6 Exciton-exciton scattering rate

Let us end this work by considering a physical quantity in which these $`H`$ matrix elements between two-exciton states may appear in a direct way. Ciuti et al (CSPQS) have proposed (11) to calculate the elastic Coulomb scattering of 1s excitons by 1s excitons, namely

$$(1\mathrm{s},𝐐,S)+(1\mathrm{s},𝐐^{},S^{})(1\mathrm{s},𝐐+𝐪,S_f)+(1\mathrm{s},𝐐^{}𝐪,S_f^{}),$$
(62)
(see their Eq (3)), from the matrix element of the exact hamiltonian $`H`$ between these two-exciton states. Let us rewrite this Eq (62) as $`(i)+(j)(l)+(n)`$. In the following we will drop the angular momentum parts for simplicity, as they are unimportant for the problem we raise.

The one-exciton wave function $`\varphi _𝐐(𝐫_e,𝐫_h)`$ given in their Eq (1) is just our $`\varphi _i(𝐫_e,𝐫_h)`$ (see Eq (26)), while the two-exciton wave function $`\varphi _{\mathrm{𝐐𝐐}^{}}^{SS^{}}`$ given in their Eq (4) reads, within our notations,

$$\varphi _{ij}(𝐫_e,𝐫_h,𝐫_e^{},𝐫_h^{})=\frac{1}{2}[\varphi _i(𝐫_e,𝐫_h)\varphi _j(𝐫_e^{},𝐫_h^{})+(𝐫_e,𝐫_h𝐫_e^{},𝐫_h^{})(𝐫_e𝐫_e^{})(𝐫_h𝐫_h^{})],$$
(63)
so that it is just the wave function of the state $`B_i^{}B_j^{}|v`$.

It is possible to show directly from Eq (63), i. e. their Eq (4), that

$$\varphi _{ln}|\varphi _{ij}=\delta _{li}\delta _{nj}+\delta _{lj}\delta _{ni}2\lambda _{lnij},$$
(64)
with

$$\lambda _{lnij}=\frac{1}{2}d𝐫_ed𝐫_hd𝐫_e^{}d𝐫_h^{}[\varphi _l^{}(𝐫_e,𝐫_h)\varphi _n^{}(𝐫_e^{},𝐫_h^{})+(ln)]\varphi _i(𝐫_e,𝐫_h^{})\varphi _j(𝐫_e^{},𝐫_h),$$
(65)
which is nothing but our Eqs (30,33) and (52). As a consequence, the $`|\varphi _{ln}`$ states are *not* normalized. Eq (64) also shows that the initial and final states of the exciton-exciton scattering considered in Eq (62) are *not* orthogonal for $`𝐪\mathrm{𝟎}`$ or $`𝐪𝐐^{}𝐐`$. So that, even if for such $`𝐪`$, the $`\delta `$ terms of Eq (64) give zero, this matrix element is not zero due to the $`\lambda _{lnij}`$ term.

According to CSPQS, the scattering amplitude corresponding to the process of Eq (62) should be equal to $`H_{SS^{}}^{S_fS_f^{}}(𝐐,𝐐^{},𝐪)=\varphi _{ij}|H|\varphi _{ln}`$ (see their Eq (6)), where $`H`$ is the exact two-electron and two-hole semiconductor hamiltonian. Their Eq (6) is nothing but $`v|B_iB_jHB_l^{}B_n^{}|v`$ written in $`𝐫`$ space. Our equations (53,54) immediately give

$`v|B_iB_jHB_l^{}B_n^{}|v`$ $`=`$ $`2\left[\xi _{ijln}^{\mathrm{dir}}\xi _{ijln}^{\mathrm{left}}(E_l+E_n)\lambda _{ijln}\right]`$ (66)
$`=`$ $`2\left[\xi _{ijln}^{\mathrm{dir}}\xi _{ijln}^{\mathrm{right}}(E_i+E_j)\lambda _{ijln}\right],`$ (67)
as the $`\delta `$ terms give zero for $`𝐪\mathrm{𝟎}`$ and $`𝐪𝐐^{}𝐐`$. If we consider the results given by CSPQS in their Eq (7), we see that the first two terms given by their Eqs (12-14), exactly correspond to the two terms of $`(2\xi _{ijln}^{\mathrm{dir}})`$ (see our Eqs (21,25)), while the last two terms (their Eq (16) and its “hole” equivalent) correspond to the two terms of $`(2\xi _{ijln}^{\mathrm{right}})`$ (see our Eqs (41,46)). Consequently they should have found the result given in Eq (67). Even if calculations in $`𝐫`$ space are quite cumbersome, it is in fact possible to check the existence of the missing term $`(E_i+E_j)\lambda _{ijln}`$ directly from a (tedious) calculation of the right hand side of their Eq (6). The origin of their mistake probably comes from the fact that they have considered the two states $`|\varphi _{ln}`$ and $`|\varphi _{ij}`$ as orthogonal, which is not true (see Eq (64)).

The correct value of the $`\varphi _{ij}|H|\varphi _{ln}`$ matrix element given in Eqs (66-67) leads us to question the validity of the whole procedure to determine the exciton-exciton scattering rate.

(i) First the $`|\varphi _{ij}`$ wave functions are not normalized, so that it would be reasonable to use normalized states and relate the exciton-exciton scattering not to $`\varphi _{ij}|H|\varphi _{ln}`$ but to $`\varphi _{ij}|H|\varphi _{ln}/\varphi _{ij}|\varphi _{ij}^{1/2}\varphi _{ln}|\varphi _{ln}^{1/2}`$. This would add a prefactor to all scatterings (see Eqs (58-59)).

(ii) There is however a much more dramatic problem with CSPQS procedure. The existence of the last terms of Eqs (66-67) is highly non physical for an exciton-exciton scattering : As it contains the exciton energy $`E_l+E_n`$ (or $`E_i+E_j`$), which is essentially equal to twice the band gap, such a scattering would be band gap dependent. We could of course get rid of this problem by deciding to use for $`H`$ an hamiltonian without the band gap i. e. by dropping $`\mathrm{\Delta }`$ in Eq (3).

Actually the exciton-exciton scattering rate cannot be related to this matrix element for more fundamental reasons. In usual problems dealing with interactions, the hamiltonian can be written as $`H=H_0+V`$. The Fermi golden rule then says that the transition rate between two different $`H_0`$ eigenstates $`|i`$ and $`|f`$ results from the possible interaction between these two states through $`|f|V|i|^2`$. When these states are eigenstates of $`H_0`$, they form an orthogonal basis, so that this transition rate is *also* equal to $`|f|H|i|^2`$, since $`f|H_0|i=0`$ for $`|f|i`$. All the difficulties here come from the fact that $`H`$ cannot be written as $`H_0+V`$ : There is no hamiltonian for which the two-exciton states $`|\varphi _{ln}`$ and $`|\varphi _{ij}`$ are the (exact) eigenstates, so that (i) these two-exciton states are not orthogonal for $`(l,n)(i,j)`$ ; (ii) the matrix elements of $`H`$ between these states have no reason to be equal to the matrix element of “the” interacting potential $`V`$ – which, anyway, cannot be formally extracted from $`H`$. Consequently there is no reason to believe that the exciton-exciton scattering rate is given by $`|\varphi _{ln}|H|\varphi _{ij}|^2`$.

## 7 Conclusion

We have reconsidered our commutation technique designed to deal with interacting close-to-boson particles and extended it to excitons with spin degrees of freedom. Although more cumbersome, the expressions of the two important parameters of this commutation technique, namely $`\xi _{lnij}^{\mathrm{dir}}`$ and $`\lambda _{lnij}`$, are still rather transparent. When the “spins” of the electron and the hole making the exciton are used – instead of the total kinetic momentum of the exciton –, these two parameters appear as a product of a spin part and an orbital part which have both a *very simple physical meaning* (see Eqs (22) and (25) for $`\xi _{lnij}^{\mathrm{dir}}`$ and Eqs (31) and (33) for $`\lambda _{lnij}`$). The first one $`\xi _{lnij}^{\mathrm{dir}}`$ corresponds to all Coulomb interactions between the $`(i,j)`$ or $`(l,n)`$ excitons when, on both sides, these excitons are made with the same electron-hole pairs $`(e,h)`$ and $`(e^{},h^{})`$. The second one $`\lambda _{lnij}`$ has a purely fermionic origin and is simply related to the fact that the $`(i,j)`$ and $`(l,n)`$ excitons can be made with different pairs, $`(e,h)`$ $`(e^{},h^{})`$ and $`(e,h^{})`$ $`(e^{},h)`$. This commutation technique allows to easily calculate any matrix elements between exciton states in an *exact* way. In this paper we have explicitly calculated the matrix elements of $`H`$ between two-exciton states and we have shown that they cannot be related to the scattering rate of two excitons as previously proposed.

APPENDIX A : Calculation of the $`\gamma _{li}(𝐪)`$

Using Eq (17), $`\gamma _{li}(𝐪)`$ reads

$$\gamma _{li}(𝐪)=\beta _{\nu _l\nu _i}(\alpha _ha_x𝐪)\beta _{\nu _l\nu _i}(\alpha _ea_x𝐪),$$
(68)
with $`\beta _{\nu \nu ^{}}(𝐮)`$ defined by

$$\beta _{\nu \nu ^{}}(𝐮)=x_\nu |e^{i𝐮.𝐫/a_x}|x_\nu ^{}=\beta _{\nu ^{}\nu }^{}(𝐮),$$
(69)
and $`a_x`$ chosen to be the 3D Bohr radius.

From Eq (69), we see that $`\beta _{\nu \nu ^{}}(\mathrm{𝟎})=\delta _{\nu \nu ^{}}`$, the value $`𝐮=\mathrm{𝟎}`$ being obtained either for $`𝐪=\mathrm{𝟎}`$ or for $`\alpha _{e,h}=0`$, i. e. $`m_h`$ (or $`m_e`$) infinite.

We can note that, if the $`\nu `$ and $`\nu ^{}`$ states have the same parity, $`\beta _{\nu \nu ^{}}(𝐮)=\beta _{\nu \nu ^{}}(𝐮)`$, so that $`\gamma _{li}(𝐪)=0`$ for $`m_e=m_h`$ : There is no direct scattering towards a same parity state if the electron and hole masses are equal.

We can also note that the scattering $`\gamma _{li}(𝐪)`$ depends on $`q=|𝐪|`$ only if the $`(l,i)`$ states are S states.

Let us now calculate some values of these $`\gamma _{li}(𝐪)`$, for 3D and 2D systems, when the $`(l,i)`$ states are S states.

1) 3D case

In 3D, the exciton relative motion wave function reads

$$𝐫|x_\nu =a_x^{3/2}\phi _\nu ^{(3D)}(\rho ,\theta ,\phi ),$$
(70)
where $`\rho =r/a_x`$. The wave functions of the lowest energy S states are given by

$$\phi _{1s}^{(3D)}=e^\rho /\sqrt{\pi },\phi _{2s}^{(3D)}=(2\rho )e^{\rho /2}/\sqrt{32\pi }.$$
(71)
In 3D, Eq (69) reads for S states,

$$\beta _{\nu \nu ^{}}^{(3D)}(u)=_0^+\mathrm{}\rho ^2𝑑\rho _0^{2\pi }𝑑\phi _0^\pi \mathrm{sin}\theta d\theta e^{iu\rho \mathrm{cos}\theta }\phi _\nu ^{(3D)}(\rho ,\theta ,\phi )\phi _\nu ^{}^{(3D)}(\rho ,\theta ,\phi ).$$
(72)
By inserting Eq (71) into Eq (72), we find

$`\beta _{1s,1s}^{(3D)}(u)`$ $`=`$ $`{\displaystyle \frac{1}{(1+u^2/4)^2}},`$ (73)
$`\beta _{2s,2s}^{(3D)}(u)`$ $`=`$ $`{\displaystyle \frac{(1u^2)(12u^2)}{(1+u^2)^4}},`$ (74)
$`\beta _{1s,2s}^{(3D)}(u)`$ $`=`$ $`{\displaystyle \frac{2^{17/2}u^2}{3^6(1+4u^2/9)^3}}.`$ (75)
From these $`\beta _{\nu \nu ^{}}^{(3D)}(u)`$, we can easily deduce the $`\gamma _{li}^{(3D)}(𝐪)`$ according to Eq (68), and the direct scattering $`\xi _{lnij}^{\mathrm{dir}}`$ according to Eq (23).

2) 2D case

In 2D, the exciton relative motion wave function reads

$$𝐫|x_\nu =a_x^1\phi _\nu ^{(2D)}(\rho ,\phi ).$$
(76)
The lowest energy S state wave functions are given by

$$\phi _{1s}^{(2D)}=\left[2^{3/2}/\sqrt{\pi }\right]e^{2\rho },\phi _{2s}^{(2D)}=\left[(2/3)^{3/2}/\sqrt{\pi }\right](14\rho /3)e^{2\rho /3}.$$
(77)
For 2D S states, Eq (69) reads

$$\beta _{\nu \nu ^{}}^{(2D)}(u)=_0^+\mathrm{}\rho 𝑑\rho _0^{2\pi }𝑑\phi e^{iu\rho \mathrm{cos}\phi }\phi _\nu ^{(2D)}(\rho ,\phi )\phi _\nu ^{}^{(2D)}(\rho ,\phi ).$$
(78)
Inserting Eq (77) into Eq (78), we find

$`\beta _{1s,1s}^{(2D)}(u)`$ $`=`$ $`{\displaystyle \frac{1}{(1+u^2/16)^{3/2}}},`$ (79)
$`\beta _{2s,2s}^{(2D)}(u)`$ $`=`$ $`{\displaystyle \frac{1(27u^2/16)+(81u^4/256)}{(1+9u^2/16)^{7/2}}},`$ (80)
$`\beta _{1s,2s}^{(2D)}(u)`$ $`=`$ $`{\displaystyle \frac{3^{7/2}u^2}{2^9(1+9u^2/64)^{5/2}}}.`$ (81)
From these $`\beta _{\nu \nu ^{}}^{(2D)}(u)`$, we can deduce the $`\gamma _{li}^{(2D)}(𝐪)`$ according to Eq (68) and the direct scattering $`\xi _{lnij}^{\mathrm{dir}}`$ according to Eq (23). For direct scatterings within the 1S states only, we recover the result given in ref. (11).

APPENDIX B : Expressions of $`\widehat{\xi }_{lnij}^{\mathrm{dir}}`$, $`\widehat{\lambda }_{lnij}`$ and $`\widehat{\xi }_{lnij}^{\mathrm{right}}`$ in $`𝐫`$ space

In order to show the equivalence between Eqs (24) and (25), it is simpler to start from Eq (25). Let us consider the $`V_{ee^{}}`$ term. By taking the Fourier transform of $`V_{ee^{}}`$, namely

$$V_{ee^{}}=\underset{𝐪\mathrm{𝟎}}{}V_𝐪e^{i𝐪.(𝐫_e𝐫_e^{})},$$
(82)
and the Fourier transform of the relative motion part of the wave function given in Eq (26), namely

$$\varphi _i(𝐫_e,𝐫_h)=\frac{1}{𝒱}\underset{𝐤_i}{}𝐤_i|x_{\nu _i}e^{i𝐊_i^e.𝐫_e}e^{i𝐊_i^h.𝐫_h},$$
(83)
with $`𝐊_i^e=𝐤_i+\alpha _e𝐐_i`$ and $`𝐊_i^h=𝐤_i+\alpha _h𝐐_i`$, we can rewrite the $`V_{ee^{}}`$ term of Eq (25) as

$`\widehat{\xi }_{ee^{}}^{\mathrm{dir}}({}_{nj}{}^{li})={\displaystyle \underset{𝐤_l,𝐤_n,𝐤_i,𝐤_j}{}}x_{\nu _l}|𝐤_lx_{\nu _n}|𝐤_n𝐤_i|x_{\nu _i}𝐤_j|x_{\nu _j}{\displaystyle \underset{𝐪\mathrm{𝟎}}{}}V_𝐪{\displaystyle \frac{1}{𝒱^4}}{\displaystyle }d𝐫_ed𝐫_hd𝐫_e^{}d𝐫_h^{}`$
$`\times e^{i\left[(𝐊_l^e+𝐊_i^e+𝐪).𝐫_e+(𝐊_l^h+𝐊_i^h).𝐫_h+(𝐊_n^e+𝐊_j^e𝐪).𝐫_e^{}+(𝐊_n^h+𝐊_j^h).𝐫_h^{}\right]}.`$ (84)
The integral being equal to $`\delta _{𝐊_l^e,𝐊_i^e+𝐪}\delta _{𝐊_l^h,𝐊_i^h}\delta _{𝐊_n^e,𝐊_j^e𝐪}\delta _{𝐊_n^h,𝐊_j^h}`$, we immediately recover the first term of Eq (24).

By transforming in the same way the other terms in $`V_{hh^{}}`$, $`V_{eh^{}}`$, $`V_{e^{}h}`$, it is easy to derive Eq (24) from Eq (25).

In order to show the equivalence of Eqs (32) and (33), we can proceed similarly : By inserting Eq (83) into Eq (33) and by performing the integral, we immediately obtain Eq (32).

We do the same to show the equivalence of Eqs (42) and (46) : We insert Eqs (82,83) into Eq (46) and perform the integral. This immediately gives the first term of Eq (42). The same transformation of the terms in $`V_{hh^{}}`$, $`V_{eh^{}}`$, $`V_{e^{}h}`$ leads to the three other terms of Eq (42).

APPENDIX C : Explicit calculation of $`F_{1s1s1s1s}(\mathrm{𝟎},\mathrm{𝟎})`$

From Eq (36), we get

$$F_{1s1s1s1s}(\mathrm{𝟎},\mathrm{𝟎})=\underset{𝐤}{}|𝐤|x_{1s}|^4.$$
(85)
The Fourier transforms of the relative motion wave functions $`𝐤|x_{1s}`$ are respectively given, in 3D and 2D, by

$`𝐤|x_{1s}^{(3D)}`$ $`=`$ $`{\displaystyle \frac{8\sqrt{\pi }a_x^{3/2}}{\sqrt{𝒱}(1+a_x^2k^2)^2}},`$ (86)
$`𝐤|x_{1s}^{(2D)}`$ $`=`$ $`{\displaystyle \frac{\sqrt{2\pi }a_x}{\sqrt{𝒮}\left(1+\frac{a_x^2k^2}{4}\right)^{3/2}}},`$ (87)
$`a_x`$ being the 3D exciton Bohr radius.

In 3D we thus obtain

$$F_{1s1s1s1s}^{(3D)}(\mathrm{𝟎},\mathrm{𝟎})=\frac{2^{11}a_x^3}{𝒱}_0^+\mathrm{}\frac{x^2dx}{(1+x^2)^8}=\frac{33\pi a_x^3}{2𝒱}.$$
(88)
In 2D we get

$$F_{1s1s1s1s}^{(2D)}(\mathrm{𝟎},\mathrm{𝟎})=\frac{2\pi a_x^2}{𝒮}_0^+\mathrm{}\frac{xdx}{\left(1+\frac{x^2}{4}\right)^6}=\frac{4\pi a_x^2}{5𝒮}.$$
(89)
REFERENCES

(1) See for instance : H. HAUG and S. KOCH, Quantum Theory of the Optical and Electronic Properties of Semiconductors, World Scientific, London (1990).

(2) V.M. AXT, K. VICTOR and A. STAHL, Phys. Rev. B 53, 7244 (1996).

(3) W. SCHÄFER, D.S. KIM, J. SHAH, T.C. DAMEN, J.E. CUNNINGHAM, K.W. GOOSSEN, L.N. PFEIFFER and K. KÖHLER, Phys. Rev. B 53, 16429 (1996).

(4) V.M. AXT and S. MUKAMEL, Rev. Mod. Phys. 70, 145 (1998).

(5) T. MEIER and S.W. KOCH, Phys. Rev. B 59, 13202 (1999).

(6) See for instance : A. KLEIN and E.R. MARSHALEK, Rev. Mod. Phys. 63, 375 (1991).

(7) See for instance : H. HAUG and S. SCHMITT-RINK, Progr. Quant. Electron. 9, 3 (1984) ; more precisely their Eqs (5. 7,9).

(8) M. COMBESCOT and O. BETBEDER-MATIBET, Europhys. Lett. 58, 87 (2002).

(9) Recently, a “kinetic correction” has been added to the usual boson-exciton effective hamiltonian by G. ROCHAT, C. CIUTI, V. SAVONA, C. PIERMAROCCI and A. QUATTROPANI (Phys. Rev. B 61, 13856 (2000)). It may look like a fermionic term. It is however intrinsically incorrect as it depends on free electron and free hole energies ($`E_k^c`$ and $`E_k^v`$ with their notations), while it should depend on exciton energies $`E_i`$ (see our Eq (56)).

(10) A recent work by S. BEN-TABOU DE LEON and B. LAIKHTMAN (Phys. Rev. B 63, 125306-1 (2001)) uses a new bosonisation procedure. It is still incorrect as it relies on a separation between a subspace with all the excitons in the ground state and a subspace with at least one exciton in an excited state. Such a separation is meaningless due to our Eq (34).

(11) C. CIUTI, V. SAVONA, C. PIERMAROCCI, A. QUATTROPANI and P. SCHWENDIMANN, Phys. Rev. B 58, 7926 (1998).

(12) J. M. LUTTINGER and W. KOHN, Phys. Rev. 97, 969 (1955).

(13) J. M. LUTTINGER, Phys. Rev. 102, 1030 (1956).

(14) M. COMBESCOT and P. NOZIERES, Solid State Com. 10, 301 (1972).

(15) In 2D, and in the particular case $`\nu _l=\nu _n=\nu _i=\nu _j=1S`$, by setting $`𝐐_l𝐐_i=𝐪`$, we can check that $`\widehat{\xi }^{\mathrm{dir}}`$ is nothing but the quantity $`H_{dir}(q)`$ introduced by Ciuti et al (11) in their Eq (12). By inserting our Eqs (68) and (79) into our Eq (23), we recover the result they give in their Eqs (19-20).

FIGURE CAPTIONS

Fig. (1)

Direct Coulomb scattering of an $`i`$ exciton into an $`l`$ exciton, while a $`j`$ exciton is scattered into an $`n`$ exciton. As Coulomb interaction conserves the angular momenta of electrons and holes, we must have $`s_l=s_i`$, $`s_n=s_j`$, $`m_l=m_i`$, $`m_n=m_j`$. As it also conserves momenta, we must have $`𝐐_i+𝐐_j=𝐐_l+𝐐_n`$, between the centers of mass of the “in” and “out” excitons. Finally this direct Coulomb scattering contains the quantities $`\gamma _{li}(𝐪)`$ and $`\gamma _{nj}(𝐪)`$, with $`𝐪=𝐐_l𝐐_i`$, which characterize the scattering of a $`\nu _i`$ exciton into a $`\nu _l`$ state and a $`\nu _j`$ exciton into a $`\nu _n`$ state, under a $`𝐪`$ excitation (see Eq (17)).

Fig. (2)

Bare exchange coefficient : The $`i`$ and $`l`$ excitons are made with the same electron but different holes ; we do have $`s_l=s_i`$ but $`m_l=m_j`$.

Fig. (3)

The dimensionless exchange coefficient $`\mathrm{\Lambda }(p)`$, given by Eq (38), as a function of $`P=p\lambda _{2D}`$, for 2D excitons and three different values of $`m_e/m_h`$. Thick solid line : $`m_e/m_h=0`$. Thin solid line : $`m_e/m_h=0.5`$. Dashed line : $`m_e/m_h=1`$.

Fig. (4)

Right exchange Coulomb scattering (a), and left exchange Coulomb scattering (b). In the right exchange scattering, the exchange of holes making the $`(i,j)`$ excitons is made *before* the direct Coulomb interaction which then scatters the $`(p,r)`$ excitons into $`(l,n)`$ states. Note that the electron-hole interactions are *between* the $`(l,n)`$ excitons but *inside* the $`(i,j)`$ excitons. Exchange and Coulomb processes conserving spins, we do have $`s_l=s_i`$, $`s_n=s_j`$ along with $`m_l=m_j`$, $`m_n=m_i`$, the holes being crossed in the process. Note that the $`(ln)`$ change appearing in the definition of $`\xi _{lnij}^{\mathrm{dir}}`$ and $`\lambda _{lnij}`$ generates a similar term in which the electrons are crossed instead of the holes. This restores the electron-hole symmetry.

Fig. (5)

The dimensionless difference between left and right exchange scatterings $`\mathrm{\Delta }g(p)`$, given by Eq (57), as a function of $`P=p\lambda _{2D}`$, for 2D excitons and three different values of $`m_e/m_h`$. Thick solid line : $`m_e/m_h=0`$. Thin solid line : $`m_e/m_h=0.5`$. Dashed line : $`m_e/m_h=1`$.


format: text

text:
# Statistics of finite-time Lyapunov exponents in a random time-dependent potential

## I Introduction

In this work, we give a uniform description of the complete asymptotic statistics of the finite-time Lyapunov exponent for a particle moving in a random time-dependent potential. The Lyapunov exponent $`\lambda _{\mathrm{}}`$ characterizes the sensitivity of trajectories to small perturbations of the initial conditions and plays a fundamental role in the characterization of systems which display deterministic chaos ott . The Lyapunov exponent is defined in the joint limits of vanishing initial perturbation and infinitely large times. In a hyperbolic Hamiltonian system $`\lambda _{\mathrm{}}`$ may be obtained from any non-periodic trajectory, because for arbitrarily long times the trajectories uniformly explore the complete phase space.

A widely studied generalization of $`\lambda _{\mathrm{}}`$ is the finite-time Lyapunov exponent ott ; fujisaka ; grassberger0 ; Badii ; tel ; grassberger ; sepulveda ; Vaienti ; romeiras ; vulpiani ; amitrano ; eckhardt ; ernst ; adrover ; Prasad ; Diakonos ; yamada , which is defined for finite stretches (time interval $`t`$) of trajectories (generalizations to finite perturbations also exist aurell ). The sensitivity of the dynamics to initial perturbations is given by the stability matrix map $`M`$, which is the linearization of the map of initial coordinates to final coordinates. In terms of elements $`M_{ij}`$ of $`M`$, the (complex) finite-time Lyapunov exponent may then be defined as

$$\lambda =\frac{1}{t}\mathrm{ln}M_{ij}.$$
(1)
In contrast to $`\lambda _{\mathrm{}}`$, $`\lambda `$ is not a unique number independent of the initial conditions, but a fluctuating quantity with a distribution function $`P(\lambda ;t)`$ (defined by uniformly sampling all initial conditions in phase space). This distribution function determines, e. g., the generalized entropy and dimension spectra of dynamical systems ott , and more practically the weak-localization correction to the conductance al1  and the shot-noise suppression al2 ; oberholzer  in mesoscopic systems. Finite-time Lyapunov exponents also determine the wavefront stability of acoustic and electromagnetic wave propagation through a random medium, in the ray-acoustics/ray-optics regime of short wave lengths (for a recent application see Refs. wolfson ; brown ). Moreover, they have shifted into the focus of attention due to recent advances in the understanding of the role of the Lyapunov exponents for quantum-chaotic wave propagation jalabert ; cucchietti01a ; cerruti ; prosen ; jacquod ; cucchietti01 ; jacquod02 : It has been observed that under certain conditions the Lyapunov exponent can be extracted from the decay of the overlap of two wavefunctions which are propagated by two slightly different Hamiltonians (the so-called Loschmidt echo). Since the overlap is studied as a function of time, the distribution of the finite-time Lyapunov exponent is directly relevant for these investigations. This extends also to related semiclassical time scales, like to the Ehrenfest time $`(\mathrm{log}\mathrm{})/\lambda `$, which is a semiclassical estimate of the diffraction time of wave packets due to the chaotic classical dynamics.

In the limit of infinite time $`t`$ the distribution function $`P(\lambda ;t)`$ in a completely chaotic phase space tends to the limiting form $`P(\lambda ;\mathrm{})=\delta (\lambda \lambda _{\mathrm{}})`$. For large but finite $`t`$ the bulk of the distribution function can be approximated by a Gaussian centered around $`\lambda _{\mathrm{}}`$, with the width vanishing $`t^{1/2}`$ as $`t\mathrm{}`$. However, many of the properties determined by $`P(\lambda ;t)`$ (like the generalized entropy and dimension spectra) cannot be calculated from the Gaussian bulk of the distribution function ott .

In this paper we investigate for a particle moving in a one-dimensional random time-dependent potential how $`P(\lambda ;t)`$ approaches the limiting distribution function $`P(\lambda ;\mathrm{})=\delta (\lambda \lambda _{\mathrm{}})`$ for large times. Our approach uniformly applies both to the bulk as well as to the far tail $`\lambda \lambda _{\mathrm{}}`$ of the distribution function. We find that the cumulant-generating function of $`P(\lambda ;t)`$,

$$\eta (\xi )=\mathrm{ln}\mathrm{exp}(\xi t\lambda )=\underset{n=1}{\overset{\mathrm{}}{}}\lambda ^n\frac{(\xi t)^n}{n!},$$
(2)
(where the average $``$ is over initial conditions and $``$ denotes the cumulants), takes the asymptotic form

$$\eta (\xi )=\mu (\xi )t/t_c+𝒪(t^0),$$
(3)
with $`\mu (\xi )`$ a universal function (within the statistical model) and

$$t_c=\lambda _{\mathrm{}}^1\mu ^{(1)}$$
(4)
a system-specific time-scale which can be determined from the infinite-time-Lyapunov exponent and the constant $`\mu ^{(1)}=d\mu /d\xi |_{\xi =0}`$ (by definition, $`d\eta /d\xi |_{\xi =0}=\lambda _{\mathrm{}}t`$). The function $`\mu (\xi )`$ is given by the leading eigenvalue of a second-order differential equation in which $`\xi `$ appears as a parameter. This eigenvalue can be calculated perturbatively in $`\xi `$, which generates the cumulants of $`\lambda `$. The values of $`\mu `$ at integer $`\xi `$ determine the asymptotic growth rates $`(1/t)\mathrm{ln}M_{ij}^\xi =\mu (\xi )/t_c`$ of positive moments of elements of the stability matrix. We find that these values are given by the leading eigenvalue of finite-dimensional matrices.

A random time-dependent potential is often considered as a statistical model for the ergodic properties of hyperbolic chaotic motion, in the spirit of the early work of Chirikov chirikov . The time dependence of the potential may be considered to mimic the dependence of the potential in the eigentime along the trajectory. In the context of finite-time Lyapunov exponents there have been indications that a statistical description is usually valid for the chaotic background of its distribution Diakonos , while system-specific deviations may exists in some exceptional cases even in the bulk of the distribution function Prasad . While the statistical model considered in this work is tailored to a specific class of Hamiltonian systems, it can be modified straightforwardly to other classes of chaotic systems (this is briefly described at the end of this paper).

The problem of finite-time Lyapunov exponents in the random time-dependent potential is equivalent to the problem of wave-localization in a random one-dimensional potential Anderson ; Borland ; Landauer ; Mott ; reviews , because the equations of motion for the matrix elements $`M_{ij}`$ are formally equivalent to the corresponding Schrödinger equation vulpiani ; crisanti . Indeed, the Fokker-Planck equation employed in this work is based on the phase formalism described, e.g., in Ref. Frisch ; Halperin ; lifshitz . Hence, the asymptotic statistics of the finite-time Lyapunov exponent presented in this work directly is of interest and can be transferred to this field of research. A number of additional areas of application of our method come into scope if one considers the vast arena of problems which can be analyzed by products of random matrices, since the finite-time Lyapunov exponents are a valuable way to characterize the eigenvalues of these products vulpiani .

The plan of this paper is as follows: In Sec. II we formulate the problem of finite-time Lyapunov exponents in the one-dimensional random time-dependent potential. In Sec. III we show how the cumulant-generating function can be related to the parameterized eigenvalue of a second-order differential equation, and that the cumulants can be calculated systematically. Positive moments of $`M_{ij}`$ are calculated in Sec. IV. We close the paper with discussion and conclusions in Sec. V.

## II Formulation of the problem

### II.1 Statistical model

Let us consider a time-dependent Hamiltonian system with one degree of freedom (canonically conjugated coordinates $`x`$, $`p`$), and the Hamiltonian given by

$$H=\frac{p^2}{2m}+V(x,t)+\frac{V_2}{2}x^2.$$
(5)
Here $`V(x,t)`$ is a time-dependent potential and $`m`$ is a mass. We also allow for an additional static potential with curvature $`V_2`$ acting in the background of the random potential (this potential is repulsive for $`V_2<0`$ and attractive for $`V_2>0`$).

We introduce the map $`_t`$ which propagates initial conditions $`(x_i,p_i)`$ over a time interval $`t`$ to the final coordinates $`(x_f,p_f)=_t(x_i,p_i)`$. The stability matrix $`M`$ is the linearization of the map $`_t`$ and describes the sensitivity of the final coordinates to a small perturbation of the initial conditions,

$$M=\frac{(x_f,p_f)}{(x_i,p_i)}=\left(\begin{array}{cc}M_{11}& M_{12}\\ M_{21}& M_{22}\end{array}\right).$$
(6)
Area preservation of the dynamics in phase space entails the property $`detM=1`$ of the stability matrix.

We are interested in the evolution of the stability matrix with given initial conditions and increasing time interval $`t`$. According to Hamilton’s equations of motion the stability matrix fulfills the differential equation

$$\frac{dM}{dt}=KM,K=\left(\begin{array}{cc}0& m^1\\ v& 0\end{array}\right),$$
(7)
where the function $`v(t)`$ in the matrix $`K`$ is given by

$$v=V_2\frac{d^2V}{dx^2}|_{(x,p)=(x_f,p_f)}.$$
(8)
This differential equation is supplemented by the initial conditions

$$M(0)=\mathrm{diag}(1,1),$$
(9)
corresponding to the identification of the initial and final coordinate systems for $`t=0`$.

In order to study the statistical behavior of the stability matrix we now assume that $`v(t)`$ is a randomly fluctuation function equivalent to Gaussian random $`\delta `$-correlated noise,

$$v(t)=V_2,v(t_1)v(t_2)=2D\delta (t_1t_2).$$
(10)
The condition of a vanishing mean of the time-dependent part of $`v`$ corresponds to the observation that the incidence of positive and negative curvature of the potential landscape along a typical chaotic trajectory should be identical. The $`\delta `$-function correlations are valid if the correlation time of the fluctuations is smaller than the mean free transport time in the random potential. The constant $`D`$ (similar to a diffusion constant, but not identical with conventional diffusion constants of motion in phase space) can be related to the strength of the temporal fluctuations of the potential $`V(x,t)`$. However, both $`D`$ as well as the mass $`m`$ can be eliminated from the subsequent analysis by rescaling quantities in the following way:

$`t=t_ct^{},v=(D/m)t_cv^{},V_2=(D/m)t_cV_2^{}`$
$`M_{12}=(t_c/m)M_{12}^{},M_{21}=(m/t_c)M_{21}^{},`$
$`M_{11}=M_{11}^{},M_{22}=M_{22}^{}.`$ (11)
Here we defined the characteristic time scale

$$t_c=m^{2/3}D^{1/3}.$$
(12)
\[In the course of our analysis we will see that this time scale also can be found from Eq. (4).\] The rescaled (primed) quantities fulfill Eqs. (7), (9), (10) with $`D=m=1`$. Also note that the rescaling leaves the property $`detM=1`$ invariant.

### II.2 Relation to one-dimensional localization

The set of linear first-order differential equations (7) can be decoupled by converting them into second-order differential equations. It is useful to note (as mentioned in the introduction) that the equations for the elements $`M_{11}`$ and $`M_{12}`$ are equivalent to the Schrödinger equation, at energy $`E=V_2/m`$, of a particle of mass $`\mathrm{}^2/2`$ in a one-dimensional random potential $`(v+V_2)/m`$ (of vanishing mean), with $`t`$ playing the role of the spatial coordinate,

$$\frac{d^2M_{11}}{dt^2}=\frac{v}{m}M_{11},\frac{d^2M_{12}}{dt^2}=\frac{v}{m}M_{12},$$
(13)
while the other matrix elements are directly related to them by

$$M_{21}=m\frac{dM_{11}}{dt},M_{22}=m\frac{dM_{12}}{dt}.$$
(14)
The problem of finite-time Lyapunov exponents hence is closely related to the problem of one-dimensional localization in a random potential, in which the Lyapunov exponent corresponds to the inverse decay length of the wave function.

## III Cumulants of the finite-time Lyapunov exponent

We now solve the problem of finding the probability distribution function of matrix elements $`M_{ij}`$ within the statistical model of chaotic dynamics, defined by the evolution equation (7) for $`M`$, with initial condition (9), and the statistical properties (10) of the random function $`v`$. For the sake of definiteness we will consider in this section the statistics of the upper diagonal element $`M_{11}`$. The results directly carry over to the other elements of $`M`$, as is discussed in Sec. IV.3.

### III.1 Cumulant-generating function as an eigenvalue

We introduce the quantities

$$u=\mathrm{ln}M_{11}^{},z=\frac{M_{21}^{}}{M_{11}^{}},$$
(15)
where the relation $`u=\lambda t`$ to the finite-time Lyapunov exponent $`\lambda `$ is established by Eq. (1) \[note that $`M_{11}=M_{11}^{}`$ in the rescaling Eq. (11)\]. According to Eqs. (7) and (11), $`u`$ and $`z`$ fulfill the differential equations

$$\frac{du}{dt^{}}=z,\frac{dz}{dt^{}}=v^{}z^2.$$
(16)
Note that the evolution equation of $`z`$ decouples from $`u`$ and can be interpreted as a Langevin equation. Hence the distribution $`P(z;t^{})`$ can be calculated from a Fokker-Planck equation, which was considered before in the context of wavefunction localization Halperin ; lifshitz ,

$`_t^{}P(z;t^{})`$ $`=`$ $`_zP(z;t^{}),`$ (17a)
$`_z`$ $`=`$ $`_z(z^2+V_2^{}+_z).`$ (17b)
For large $`t^{}`$ the distribution function $`P(z;t^{})`$ approaches the stationary solution Frisch ; Halperin ; lifshitz

$`P_{\text{stat}}(z)`$ $`=`$ $`\stackrel{~}{N}{\displaystyle _{\mathrm{}}^z}𝑑yK(y,z),`$ (18a)
$`K(y,z)`$ $`=`$ $`e^{(y^3z^3)/3+V_2^{}(yz)},`$ (18b)
$`\stackrel{~}{N}`$ $`=`$ $`\pi ^2[\mathrm{Ai}^2(V_2^{})+\mathrm{Bi}^2(V_2^{})]^1.`$ (18c)
Here $`\mathrm{Ai}`$ and $`\mathrm{Bi}`$ are Airy functions. The normalization constant is directly related to the integrated density of states in the localization problem Frisch ; Halperin ; lifshitz . For $`V_2^{}=0`$, $`\stackrel{~}{N}=3^{5/6}2^{1/3}\pi ^{1/2}/[\mathrm{\Gamma }(1/6)]`$. Because $`du/dt=z/t_c`$ it is clear lifshitz  that the infinite-time Lyapunov exponent can be obtained from $`\lambda _{\mathrm{}}=z/t_c`$; this relation will be demonstrated explicitely in Sec. III.3.

The Fokker-Planck equation for the joint distribution function $`P(u,z;t^{})`$ is given by

$$_t^{}P=z_uP+_zP.$$
(19)
This Fokker-Planck equation with $`V_2^{}=0`$ has been derived in Ref. al1  for the autonomous chaotic scattering of a particle from a dilute collection of scatterers (with more than one degree of freedom).

The joint distribution function $`P(u,z;t^{})`$ does not approach a stationary limit because $`u`$ runs away to infinitely large values. In order to analyze the behavior of the distribution function $`P(u,z;t^{})`$ for large times we convert the Fokker-Planck equation (19) into an eigenvalue problem which discriminates between the different time scales involved in this evolution. For this purpose, we introduce into Eq. (19) the ansatz

$$P(u,z;t^{})=_i\mathrm{}^{+i\mathrm{}}\frac{d\xi }{2\pi i}\underset{n=0}{\overset{\mathrm{}}{}}\mathrm{exp}(\mu _nt^{}\xi u)f_n(\xi ,z).$$
(20)
(The integration contour along the imaginary axis corresponds to an inverse Laplace transformation.) It follows that the functions $`f_n`$ fulfill the differential equation

$$\mu _nf_n(\xi ,z)=(\xi z+_z)f_n(\xi ,z),$$
(21a)
in which $`\mu _n`$ and $`\xi `$ appear as parameters. However, in order to obtain a meaningful probability distribution function (20) we have to impose boundary conditions on $`f_n(\xi ,z)`$ at $`z\pm \mathrm{}`$. It is convenient to express these boundary conditions by the requirement
$$𝒫_{\mathrm{}}^{\mathrm{}}𝑑zf_n(\xi ,z)z<\mathrm{}.$$
(21b)
Here $`𝒫`$ denotes the principal value with respect to the integration boundaries at $`\pm \mathrm{}`$. Condition (21b) follows from the behavior $`z(t^{}t_{\mathrm{}}^{})^1`$ of the solution of the differential equation (16) close to times $`t^{}t_{\mathrm{}}^{}`$ where $`|z|\mathrm{}`$ (and hence $`v^{}`$ can be ignored). In practical terms, the condition (21b) guarantees that the drift of $`u`$ remains finite for all times.

Eqs. (21) form an eigenvalue problem, since condition (21b) only can be fulfilled for a discrete set of numbers $`\mu _n`$—note that these eigenvalues depend on the parameter $`\xi `$. In the limit of large $`t^{}`$ only the largest eigenvalue $`\mu _0(\xi )\mu (\xi )`$ is relevant, because the other eigenvalues give rise to exponentially smaller contributions. This eigenvalue vanishes as $`\xi 0`$, i. e., $`\mu (0)=0`$, because the stationary distribution of $`z`$, Eq. (18), must be recovered for large times from Eq. (20) by integrating out $`u`$.

The moments of $`u`$ are given by

$`u^n`$ $`=`$ $`{\displaystyle _i\mathrm{}^{+i\mathrm{}}}{\displaystyle \frac{d\xi }{2\pi i}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑u\mathrm{exp}(\mu t^{}\xi u)f(\xi )u^n`$ (22)
$`=`$ $`\underset{\xi 0}{lim}_\xi ^n\mathrm{exp}(\mu (\xi )t^{})f(\xi ),`$
where the coefficients $`f(\xi )=_{\mathrm{}}^{\mathrm{}}𝑑zf_0(\xi ,z)`$ are determined, in principle, by the initial condition for $`P(u,z;t^{})`$ at $`t^{}=0`$. From Eq. (22) we obtain the moment-generating function

$$\chi (\xi )=\mathrm{exp}(\xi u)=\mathrm{exp}(\mu (\xi )t/t_c)f(\xi ),$$
(23)
where we re-introduced the original time variable $`t=t_ct^{}`$ by Eq. (11). The cumulant-generating function (2) hence takes the form of Eq. (3), including the corrections of order $`t^0`$,

$$\eta (\xi )=\mathrm{ln}\chi (\xi )=\mu (\xi )t/t_c+\mathrm{ln}f(\xi ).$$
(24)
The cumulants $`\lambda ^n`$ of the finite-time Lyapunov exponent are obtained by expanding the generating function $`\eta `$ in powers of $`\xi `$, see Eq. (2). In terms of the coefficients of the Taylor expansion

$$\mu =\underset{n=1}{\overset{\mathrm{}}{}}\xi ^n\mu ^{(n)}$$
(25)
\[which starts with the linear term in $`\xi `$ because $`\mu (0)=0`$\], according to Eqs. (2) and (24) the $`n`$th cumulant of $`\lambda `$ is then given by

$`\lambda ^n=n!\mu ^{(n)}t_c^1t^{1n}+𝒪(t^n).`$ (26)
This equation means that within the statistical model the cumulants are universal quantities in the leading order in $`t`$, in the sense that the initial conditions $`P(z,u;0)`$ only enter the next-order corrections. The only system-specific parameters which enter the cumulants are the time scale $`t_c`$ and the (rescaled) strength $`V_2^{}`$ of the static potential. Note that ratios of cumulants are even independent of the time scale $`t_c`$ (and hence of the parameters $`D`$ and $`m`$ of the statistical model).

The form (4) of $`t_c`$ follows from Eq. (26) when $`t_c`$ is expressed in terms of the infinite-time Lyapunov exponent $`\lambda _{\mathrm{}}`$ with help of the definition

$$\lambda _{\mathrm{}}=\underset{t\mathrm{}}{lim}\lambda =\mu ^{(1)}/t_c.$$
(27)
In terms of the bare quantities of the statistical model,

$$\lambda _{\mathrm{}}=\mu ^{(1)}D^{1/3}m^{2/3}.$$
(28)
In the next two sections we obtain general expressions for the expansion coefficients $`\mu ^{(n)}`$ and calculate explicitely the proportionality factor $`\mu ^{(1)}=d\mu /d\xi |_{\xi =0}`$ in (28), as well as the first few coefficients $`\mu ^{(2)}`$, $`\mu ^{(3)}`$, $`\mathrm{}`$, which determine, respectively, the variance and the leading non-Gaussian corrections (higher cumulants) of the fluctuations of the finite-time Lyapunov exponent around its limiting value $`\lambda _{\mathrm{}}`$.

### III.2 Recursion relations for the cumulants

We now show how the cumulants can be calculated from Eq. (26) by recursively solving a hierarchy of equations for coefficients $`\mu ^{(n)}`$ in the Taylor expansion of $`\mu (\xi )`$, Eq. (25).

In analogy to Eq. (25) let us also expand the function $`f_0(\xi ,z)`$ in powers of $`\xi `$,

$$f_0(\xi ,z)=\underset{n=0}{\overset{\mathrm{}}{}}\xi ^nf_0^{(n)}(z).$$
(29)
With Eqs. (25) and (29) the eigenvalue problem (21) can now be written order by order in powers of $`\xi ^n`$. For $`n=0`$ we recover the stationary variant (17) of the Fokker-Planck equation (19),

$$_zf_0^{(0)}(z)=0,$$
(30)
which is solved by the stationary solution $`f_0^{(0)}(z)=P_{\text{stat}}(z)`$, Eq. (18). For $`n>1`$ the differential equations are of the form

$$_zf_0^{(n)}(z)=zf_0^{(n1)}(z)+\underset{l=1}{\overset{n}{}}\mu ^{(l)}f_0^{(nl)}(z).$$
(31)
Let us assume that we have solved the hierarchy of equations up to order $`n1`$. In the next order $`n`$ both the unknown quantities $`f_0^{(n)}`$ as well as $`\mu ^{(n)}`$ appear. The unknowns can be separated by integrating the differential equation (31) over $`z`$ from $`\mathrm{}`$ to $`\mathrm{}`$: The integrated left-hand side vanishes because of condition (21b) of the eigenvalue problem. The integrated right-hand side can be rearranged to give $`\mu ^{(n)}`$,

$$\mu ^{(n)}=_{\mathrm{}}^{\mathrm{}}𝑑z[zf_0^{(n1)}(z)\underset{l=1}{\overset{n1}{}}\mu ^{(l)}f_0^{(nl)}(z)],$$
(32a)
which only involves quantities up to order $`n1`$. Subsequently, $`\mu ^{(n)}`$ can be inserted into Eq. (31). The function
$`f_0^{(n)}(z)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^z}𝑑y{\displaystyle _{\mathrm{}}^y}𝑑xK(y,z)`$ (32b)
$`\times [xf_0^{(n1)}(x)+{\displaystyle \underset{l=1}{\overset{n}{}}}\mu ^{(l)}f_0^{(nl)}(x)]`$
\[with the kernel $`K(y,z)`$ defined in Eq. (18)\] is then obtained by solving the resulting inhomogeneous differential equation with help of the partial solution $`f_0^{(0)}(z)`$ of its homogeneous counterpart, Eq. (30). This inhomogeneous part of the functions $`f_0^{(n)}(z)`$ is fixed by the requirement that $`f_0^{(0)}(z)`$ is normalized to 1. Adding the homogeneous solution to $`f_0^{(n)}(z)`$ in any order gives rise to additional terms in all higher orders, but these combine in such a way that they drop out of the calculation of the coefficients $`\mu ^{(n)}`$, which hence are uniquely determined by Eq. (32a).

The recursion relations (32) can be iterated to calculate successively all cumulants of $`\lambda `$.

### III.3 Explicit expressions and numerical values

According to Eq. (26), the two numbers $`\mu ^{(1)}`$ and $`\mu ^{(2)}`$ determine mean and variance of the distribution function of $`\lambda `$, which then is approximated by a Gaussian. The coefficient $`\mu ^{(1)}`$ has been obtained in Ref. lifshitz  from the Fokker-Planck equation (17) for arbitrary $`V_2`$. For the special case $`V_2=0`$, the two coefficients $`\mu ^{(1)}`$ and $`\mu ^{(2)}`$ have been obtained in Ref. al1  from the Fokker-Planck equation (19). However, the deviations from the Gaussian distribution function are not at all negligible for many chaotic systems, which is most clearly displayed in their generalized dimension and entropy spectra ott . As we have seen in the previous subsection III.2, our approach of reduction to the eigenvalue problem (21) allows to analyze the non-Gaussian deviations by the higher cumulants of $`\lambda `$. \[In next section IV, we show that one can even obtain from our analysis the positive moments of $`M_{11}`$, which are determined by the far tail $`\lambda \lambda _{\mathrm{}}`$ of $`P(\lambda ;t)`$, while the bulk of the distribution is essentially irrelevant for these moments.\]

Explicit expressions for the first few coefficients $`\mu ^{(1)}`$, $`\mu ^{(2)}`$, $`\mu ^{(3)}`$, and $`\mu ^{(4)}`$ result from Eq. (32a),

$`\mu ^{(1)}`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑zzf_0^{(0)}(z),`$ (33a)
$`\mu ^{(2)}`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑z(z\mu ^{(1)})f_0^{(1)}(z),`$ (33b)
$`\mu ^{(3)}`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑z[(z\mu ^{(1)})f_0^{(2)}(z)\mu ^{(2)}f_0^{(1)}(z)],`$ (33c)
$`\mu ^{(4)}`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}dz[(z\mu ^{(1)})f_0^{(3)}(z)\mu ^{(2)}f_0^{(2)}(z)`$ (33d)
$`\mu ^{(3)}f_0^{(1)}(z)],`$
where $`f_0^{(0)}(z)=P_{\text{stat}}(z)`$ is given by the stationary distribution function of $`z`$, Eq. (18), while the other functions follow from Eq. (32b),

$`f_0^{(1)}(z)`$ $`=`$ $`{\displaystyle _{z>y>x}}𝑑y𝑑xK(y,z)(\mu ^{(1)}x)f_0^{(0)}(x),`$ (34a)
$`f_0^{(2)}(z)`$ $`=`$ $`{\displaystyle _{z>y>x}}dydxK(y,z)[(\mu ^{(1)}x)f_0^{(1)}(x)`$ (34b)
$`\mu ^{(2)}f_0^{(0)}(x)],`$
$`f_0^{(3)}(z)`$ $`=`$ $`{\displaystyle _{z>y>x}}dydxK(y,z)[(\mu ^{(1)}x)f_0^{(2)}(x)`$ (34c)
$`\mu ^{(2)}f_0^{(1)}(x)\mu ^{(3)}f_0^{(0)}(x)].`$
The coefficient $`\mu ^{(1)}`$ is then given by lifshitz

$`\mu ^{(1)}={\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{dV_2^{}}}\mathrm{log}\stackrel{~}{N},`$ (35)
where $`\stackrel{~}{N}`$ is given in Eq. (18), while the cumulants for $`n2`$ can be obtained quickly by numerical integration of $`2n`$-fold integrals. The effort of integration can be greatly reduced down to the expense equivalent to a single integral, because the integrand factorizes. An efficient recursive scheme is described in the Appendix. In Fig. 1 we plot the coefficient $`\mu ^{(1)}`$ and the ratios $`n!\mu ^{(n)}/\mu ^{(1)}`$ for $`n=2,3,4`$ as a function of $`V_2^{}`$. The non-Gaussian corrections are largest around $`V_2^{}=0`$, while they become irrelevant for large negative or positive values of $`V_2^{}`$.

For strong confinement, $`V_2^{}1`$, the coefficients $`n!\mu ^{(n)}/\mu ^{(1)}\delta _{1n}+\delta _{2n}`$, with $`\delta _{mn}`$ the Kronecker symbol, and the Gaussian approximation

$$\mu _{\text{Gaussian}}(\xi )=\mu ^{(1)}\left(\xi +\frac{1}{2}\xi ^2\right)$$
(36)
becomes valid. \[In the context of wave localization, this corresponds to the well-known limit of a large Fermi energy $`EV_2^{}`$ (cf. Sec. II.2).\]

Analytical results can be found in the case $`V_2^{}=0`$ for the first two coefficients,

$`\mu ^{(1)}`$ $`=`$ $`{\displaystyle \frac{(3/2)^{1/3}\sqrt{\pi }}{\mathrm{\Gamma }(1/6)}},`$ (37a)
$`\mu ^{(2)}`$ $`=`$ $`{\displaystyle \frac{5\pi ^2}{18}}\stackrel{~}{N}{\displaystyle \frac{\pi }{2\sqrt{3}}}\stackrel{~}{N}_3F_2(1,1,{\displaystyle \frac{7}{6}};{\displaystyle \frac{3}{2}},{\displaystyle \frac{3}{2}};{\displaystyle \frac{3}{4}}),`$ (37b)
where $`\stackrel{~}{N}(V_2^{}=0)=3^{5/6}2^{1/3}\pi ^{1/2}/[\mathrm{\Gamma }(1/6)]`$, while $`{}_{3}{}^{}F_{2}^{}`$ is a generalized hypergeometric function. Incidentally, the numerical value given for $`\mu ^{(2)}`$ in Ref. al1  is wrong, but the analytic expression given in that paper is equivalent to Eqs. (33b) and (37b). In Tab. 1 we tabulate the numerical values of the first eight coefficients $`n!\mu ^{(n)}`$ for $`V_2^{}=0`$.

## IV Positive moments

### IV.1 Formally exact expressions

In view of Eqs. (15) and (23) we find that the exponential growth rates of the positive moments of $`M_{11}`$ are given by the eigenvalue $`\mu (\xi )`$ of Eq. (21) at integer values of $`\xi `$:

$$\frac{d\mathrm{ln}M_{11}^\xi }{dt}=\frac{\mu (\xi )}{t_c}.$$
(38)
As we will now show, for integer values of $`\xi `$ the eigenvalue problem (21) can be reduced to a matrix eigenvalue problem of finite dimension. For the first few moments the leading eigenvalue can be calculated explicitely, while for larger values it is formally given by the largest root of the corresponding characteristic polynomial.

In order to obtain a solution of the differential equation (21a), we write

$$f_0(\xi ,z)=_{\mathrm{}}^z𝑑yK(y,z)\frac{g(z)}{g(y)^2}$$
(39)
\[with the kernel $`K(y,z)`$ defined in Eq. (18)\], and obtain for $`g`$ the differential equation

$$(\mu \xi z)g=(z^2+V_2^{})_zg+_z^2g$$
(40)
(a triconfluent Heun’s equation with singularity at $`1/z=0`$). We introduce into this equation the polynomial ansatz

$$g(z)=\underset{n=0}{\overset{\xi }{}}c_nz^n.$$
(41)
Power matching results in the following recursion relation

$$(\xi n)c_n=\mu c_{n+1}+(n+2)[V_2^{}c_{n+2}(n+3)c_{n+3}]$$
(42a)
for the coefficients $`c_n`$, with initial conditions
$$c_\xi =1,c_{\xi 1}=\mu ,c_{\xi 2}=\mu ^2/2.$$
(42b)
For integer $`\xi `$ this recursion relation terminates. We obtain functions $`c_0(\mu )`$, $`c_1(\mu )`$, and $`c_2(\mu )`$ and an additional condition from the term in Eq. (40) which is constant in $`z`$,

$$p_\xi (\mu )=\mu c_0+V_2^{}c_12c_2=0,$$
(43)
where $`p_\xi (\mu )`$ is a polynomial of degree $`\xi +1`$.

The polynomial $`p_\xi (\mu )`$ can also be interpreted as the characteristic polynomial of the $`(\xi +1)\times (\xi +1)`$-dimensional matrix

$$\left(\begin{array}{cccccccc}0& V_2^{}& 12& 0& 0& \mathrm{}& \mathrm{}& \mathrm{}\\ \xi & 0& 2V_2^{}& 23& 0& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& \xi 1& 0& 3V_2^{}& 34& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& \xi 2& 0& 4V_2^{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& \xi 3& 0& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& 0& (1\xi )V_2^{}& (\xi 1)\xi \\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& 2& 0& \xi V_2^{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& 0& 1& 0\end{array}\right),$$
(44)
which is the matrix representation of the eigenvalue problem (21) in the space of the monomial expansion of $`g(z)`$.

The exponential growth rate $`\mu (\xi )`$ of the $`\xi `$th moment is given by the largest root of $`p_\xi (\mu )`$, or equivalently by the largest eigenvalue of the matrix (44). In subsection IV.2 we will see for the examples $`\xi =1,2`$ that the other roots show up in the transient behavior of the moments.

First we present results in absence of the static background potential, $`V_2^{}=0`$. The values for the first few moments are given in Tab. 2. Figure 2 shows the growth rates and the real part of the subleading eigenvalue for values of $`\xi `$ up to $`80`$. A log-normal statistics of $`M_{11}`$ (corresponding to a Gaussian statistics of the finite-time Lyapunov exponents) would result in the quadratic dependence Eq. (36) of $`\mu (\xi )`$ on $`\xi `$, while the plot shows a weaker (approximately linear) dependence for large $`\xi `$. This results from the influence of the terms $`\mu ^{(n)}\xi ^n`$ for $`n3`$ in the complete Taylor expansion of $`\mu `$, Eq. (25). Further note that the subleading eigenvalue stays at a finite distance to the leading eigenvalue (indeed, their distance increases with increasing $`\xi `$), as we have assumed before in restricting the attention to the leading eigenvalue $`\mu _0`$ of the eigenvalue problem (21).

For finite $`V_2^{}`$, the growth rate of the first moment

$$\mathrm{Re}\mu (1)=|\mathrm{Im}\sqrt{V_2^{}}|$$
(45)
vanishes in the case of confinement, $`V_2^{}>0`$. This will be confirmed by the direct computation in Sec. IV.2. The growth rate of the second moment is given by

$`\mu (2)`$ $`=`$ $`2^{1/3}\left(1+\sqrt{1+16V_{2}^{}{}_{}{}^{3}/27}\right)^{1/3}`$ (46)
$`+2^{1/3}\left(1\sqrt{1+16V_{2}^{}{}_{}{}^{3}/27}\right)^{1/3}`$
\[with the roots taken such that $`\mu (2)`$ is real\]. We plotted the real parts of the leading and subleading growth rates \[eigenvalues of matrix (44)\] for the first four moments in Fig. 3.

### IV.2 Direct computation of the first and second moment

In order to illustrate the results for the growth rates of the moments $`M_{11}^\xi `$ we compare the results for $`\xi =1`$ and $`\xi =2`$ to the exact results for all times (including the transient behavior). A formal solution of the differential equation (13) in terms of a series in the disorder potential is obtained by integrating Eq. (13) twice, under observation of the initial conditions $`M_{11}=1`$, $`dM_{11}/dt=0`$ for $`t=0`$, and iterating the resulting integral relation

$`M_{11}(t)`$ $`=`$ $`1+{\displaystyle _0^t}𝑑t_1{\displaystyle _0^{t_1}}𝑑s_1{\displaystyle \frac{v(s_1)}{m}}M_{11}(s_1)`$ (47)
$`=`$ $`1+{\displaystyle _0^t}𝑑t_1(tt_1){\displaystyle \frac{v(t_1)}{m}}M_{11}(t_1).`$
The formal solution is of the form

$$M_{11}(t_0)=1+\underset{n=1}{\overset{\mathrm{}}{}}\underset{k=1}{\overset{n}{}}_0^{t_{k1}}𝑑t_k(t_{k1}t_k)\frac{v(t_k)}{m},$$
(48)
where we introduced $`t_0=t`$ for notational convenience.

For the first moment we can average Eq. (48) directly. Because of the factors $`(t_{k1}t_k)`$ and the time ordering, the random function $`v`$ never appears instantaneously in second or higher order in any of the integrals. Hence we can replace $`v`$ by its average, given in Eq. (10). Consequently, the first moment is given by

$$M_{11}=\mathrm{cos}[(t/t_c)\sqrt{V_2^{}}]=\frac{1}{2}e^{(t/t_c)\sqrt{V_2^{}}}+\frac{1}{2}e^{(t/t_c)\sqrt{V_2^{}}}.$$
(49)
For $`V_2^{}=0`$ the first moment is constant and given by its initial value, $`M_{11}=1`$. This means that negative deviations $`M_{11}0`$, corresponding to inverse hyperbolic motion, cancel precisely the positive deviations $`M_{11}0`$ of hyperbolic motion. For negative $`V_2^{}`$ the first moment grows, while for positive $`V_2^{}`$ it oscillates and stays of order unity. In the decomposition of the cosine into the two exponentials, we identify in the exponents the two roots $`\pm \sqrt{V_2^{}}`$ of the characteristic polynomial $`p_{\xi =1}(\mu )=\mu ^2+V_2^{}`$ of the matrix (44) with $`\xi =1`$. For negative $`V_2^{}`$, the subleading exponent hence governs the transient behavior of the first moment.

For the second moment let us restrict for simplicity to the case $`V_2^{}=0`$. We group the functions $`v`$ in the two factors of $`M_{11}`$ in pairs and then invoke the delta-correlations of Eq. (10). Performing the time-ordered integrals we obtain

$`M_{11}^2=1+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(2/t_c^3)^n{\displaystyle \underset{k=1}{\overset{n}{}}}{\displaystyle _0^{t_{k1}}}𝑑t_k(t_{k1}t_k)^2`$
$`={\displaystyle \frac{1}{3}}[e^{\mu (2)t/t_c}+e^{(1)^{1/3}\mu (2)t/t_c}+e^{(1)^{1/3}\mu (2)t/t_c}].`$
(50)
The asymptotic growth rate of the second moment is given by the leading root $`\mu (2)=2^{2/3}`$ of the characteristic polynomial $`p_2(\mu )=\frac{1}{2}\mu ^32`$, which is in accordance to Tab. 2. The second and third exponent are the other two roots of this polynomial.

### IV.3 Equivalence of matrix elements

So far we mainly studied the statistics of the upper diagonal element $`M_{11}`$ of the stability matrix $`M`$. At this point now we can discuss how the results for the cumulant-generating function and the positive moments can be transferred to the other elements of $`M`$.

The differential equation (14) for $`M_{22}`$ can be integrated similarly as the one for $`M_{11}`$, from which we obtain analogously to Eq. (48) the formal solution

$$M_{22}(t_0)=1+\underset{n=1}{\overset{\mathrm{}}{}}\underset{k=1}{\overset{n}{}}_0^{t_{k1}}𝑑t_k(t_kt_{k+1})\frac{v(t_k)}{m}.$$
(51)
Here we defined in each term of order $`n`$ that $`t_{n+1}=0`$. It follows by direct computation that the first two moments of $`M_{22}`$ are identical to those of $`M_{11}`$,

$$M_{22}=M_{11},M_{22}^2=M_{11}^2.$$
(52)
These explicit results already suggest that the statistics of the two diagonal matrix elements is the same. Indeed, the transformation $`t_k=t\stackrel{~}{t}_{n+1k}`$, $`v(t\stackrel{~}{t})=\stackrel{~}{v}(\stackrel{~}{t})`$ brings Eq. (51) into the form of Eq. (48) and leaves the properties of the Gaussian noise (10) invariant. Hence even the transient behavior of the diagonal elements is completely identical, for arbitrary values of $`V_2^{}`$.

The results for the cumulant-generating function $`\eta (\xi )`$ (hence also the growth rates of the moments, but not the transient behavior) can also be transferred to the offdiagonal matrix elements of $`M`$: The element $`M_{12}`$ fulfills the same differential equation as $`M_{11}`$, see Eq. (13), while $`M_{21}`$ fulfills the same differential equation as $`M_{22}`$. The initial conditions of the offdiagonal matrix elements differ from those of the diagonal elements. However, according to Eq. (24) this only affects the function $`f(\xi )`$ in the subleading corrections of the cumulant-generating function \[which, for the example of the second moment, results in factors in front of the exponential functions which are different than in Eq. (50)\].

Let us add that from Eqs. (48) and (51) we find for $`V_2^{}=0`$ the cross-correlator

$$M_{11}M_{22}=\frac{1}{2}+\frac{1}{2}M_{11}^2.$$
(53)
As a consequence, for $`V_2^{}=0`$ the trace $`\mathrm{tr}M=M_{11}+M_{22}`$ of the stability matrix has the following first two moments

$`\mathrm{tr}M`$ $`=`$ $`2,`$ (54a)
$`(\mathrm{tr}M)^2`$ $`=`$ $`1+e^{\mu (2)t/t_c}+2\mathrm{Re}e^{(1)^{1/3}\mu (2)t/t_c}.`$
## V Discussion

In this work we presented a uniform approach to the asymptotic statistics of finite-time Lyapunov exponents, for the model (described in Sec. II) of a particle moving in a random time-dependent potential. The cumulant-generating function $`\eta (\xi )`$ was found to be directly proportional to the eigenvalue $`\mu (\xi )`$ of a parameterized differential equation, defined by Eqs. (21). This facilitated an effective analysis of the statistics, including the non-Gaussian deviations of the distribution function. These deviations are especially important for the positive moments of the elements of the stability matrix, since their growth rate cannot be predicted by the Gaussian approximation Eq. (36).

We limited our attention to the case of time-dependent Hamiltonian systems with a single degree of freedom and a Hamiltonian (5) which is of the special type of kinetic energy plus potential energy, with time-dependence only in the potential energy. This case is of particular interest because of its direct applicability to specific dynamical systems as in the random wave-propagation problem of Refs. wolfson ; brown , and because of its applicability to one-dimensional wave localization. For the Hamiltonian (5) the matrix $`K`$ in the differential equation (7) is purely off-diagonal, with fluctuations only in the lower-left element. For Hamiltonians which do not separate into kinetic and potential energy, the differential equation (5) for $`M`$ involves the matrix $`K`$ in the more general form

$$K=\left(\begin{array}{cc}K_{11}& K_{12}\\ K_{21}& K_{22}\end{array}\right)=\left(\begin{array}{cc}\frac{^2H}{xp}& \frac{^2H}{p^2}\\ \frac{^2H}{x^2}& \frac{^2H}{xp}\end{array}\right).$$
(55)
A generalized statistical model now arises by introducing noise into all of the matrix elements of $`K`$. (One may also allow for correlations between the different matrix elements or for finite correlation times by introducing auxiliary variables for the noise in the standard way.)

Let us point out two particular cases for which a statistical description promises to result in direct applications to physical situations of interest. One case is more relevant to wave-function localization while the other is more relevant for chaotic dynamics.

a) The diagonal elements $`K_{11}=K_{22}=0`$ still vanish identically, but both off-diagonal elements $`K_{12}`$ and $`K_{21}`$ fluctuate with a vanishing mean. This situation appears to be related to the band-center case of one-dimensional localization in the Anderson model kappus ; goldhirsch  (where space is discretized on the lattice), since at the band-center the effective mass of the particle diverges (and hence the mean of $`K_{12}`$ vanishes).

b) Chaotic dynamics with an isotropic phase space may be modeled by independent fluctuations of all four matrix elements $`K_{ij}`$ with identical amplitude and vanishing mean. Hamiltonian dynamics gives rise to the further constraint $`K_{11}=K_{22}`$. Isotropic dynamics arises in typical chaotic maps (some maps, like the Baker map or the cat map, however, are not isotropic—the directions of stable and unstable manifolds are known by construction). Good candidates are the Poincaré surface of section of autonomous systems with two degrees of freedom, in which the motion in four-dimensional phase space is restricted to three-dimensional manifolds of constant energy and the coordinate along the flow field is taken as a time.

It would be interesting to compare the outcome of an analysis of model b) with the findings in the literature Prasad ; Diakonos  which indicate a certain degree of robustness (if not universality) of the distribution of finite-time Lyapunov exponents for typical chaotic systems.

###### Acknowledgements.

We gratefully acknowledge useful discussions with Philippe Jacquod and Holger Kantz, and especially with Steven Tomsovic who motivated us to study this problem.

## Appendix: Integrals for the higher cumulants

The cumulants of order $`n`$ result from the recursion relations Eq. (32) in the form of $`2n`$-fold integrals. Usually, the numerical evaluation of such integrals for large $`n`$ is very time-consuming, since the number of points on a grid covering the integration domain with lattice constant $`(1/N)`$, $`N1`$, grows rapidly with $`n`$ as $`N^{2n}`$. However, presently the integrand factorizes and the expense of the integration can be reduced from exponential to algebraic $`n`$-dependence $`nN`$. The principle can be demonstrated for the example of the two-fold integral

$$I^{(1)}=_{z_0}^{z_1}𝑑zI^{(2)}(z),I^{(2)}(z)=g(z)_{z_0}^z𝑑yI^{(3)}(y),$$
(56)
where $`g`$ is an arbitrary function and $`I^{(3)}`$ may itself be a multi-dimensional integral.

We introduce an index $`m`$ which denotes that the argument of a function is taken at the $`m`$th lattice point on the appropriate axis of the grid. The initial values of $`I_m^{(n)}`$ at $`m=0`$ (the lower integration boundary) are zero. We now can write recursively, by incrementally increasing the integration variables,

$`I_{m+1}^{(2)}`$ $`=`$ $`{\displaystyle \frac{g_{m+1}}{g_m}}I_m^{(2)}+{\displaystyle \frac{1}{N}}g_mI_m^{(3)},`$ (57a)
$`I_{m+1}^{(1)}`$ $`=`$ $`I_m^{(1)}+{\displaystyle \frac{1}{N}}I_{m+1}^{(2)}.`$ (57b)
Moreover, when $`I^{(3)}`$ itself is a multi-dimensional integral of type $`I^{(1)}`$, its current value can be obtained recursively in the same way as the value of $`I^{(1)}`$. Since each additional integral will give rise to only one additional equation \[similar either to Eq. (57a) or to Eq. (57b)\], the number of operations grows linearly with $`n`$, as advertised above. \[The recursion relations (57) have the additional advantage for the present problem that they avoid over- and underflow in the evaluation of the kernel $`K(y,z)=\mathrm{exp}(y^3/3+V_2^{}yz^3/3V_2^{}z)`$.\]