text:
# GRADED CONTRACTIONS OF VIRASORO ALGEBRAS

## 1 Introduction

The Virasoro algebra takes an important part in strings theory, integrable models and two-dimensional conformal field theory. Resently it has been investigated the extensions of Virasoro algebra leading to $`W`$-algebras, which contain the Virasoro algebra as a subalgebra. $`W`$-algebras are appeared in integrable models  and then in two-dimensional conformal field theory . Thereafter the notions of $`W`$-strings, $`W`$-gravitation have arisen and $`W`$-models of conformal field theory have been created . Hence the extensions of symmetries on the base of the Virasoro algebra are very fruitful.

A different way of extension of Virasoro algebra is a contraction. In this case new nonsemisimple algebras containing the Virasoro algebra as a subalgebra are arisen. Contractions of affine Kac-Moody algebras were studied in  and similar effect was obtained . Nonsemisimple algebras are used as algebras of internal symmetries of different physical models: gauge theories , WZNW models . Sugavara construction was developed for nonsemisimple algebras . The purpose of this paper is the investigation of contractions of Virasoro algebra and its representations with the help of the graded contraction method .

This paper is organized as follows. In the second section we briefly present the main ideas of graded contraction method. In Sec.3 the grading of Virasoro algebra by cyclic groups $`ZZ_2`$ is suggested. In Sec.4 we recall the base facts of reducible theory of modules over Virasoro algebra and describe their $`ZZ_2`$-grading. Sec.5 contains the results of some contractions of Virasoro algebra. In Sec.6 contractions of modules and some problems of their reducibility are investigated.

## 2 Graded contractions of Lie algebras and their representations

A Lie groups (algebras) contraction makes possible to get new Lie groups (algebras) from some initial ones. Moreover one can construct a representation of contracted algebras from a representation of starting algebra. The graded contractions  are defined in such a way that the grading of both Lie algebra and its representation are preserved.

First we recall some definitions. Let a Lie algebra $`L`$ is graduated by an abelian group $`G`$ as follows

$$L=L_j,[L_j,L_i]L_{j+i},j,iG.$$
(1)
A Lie algebra $`L^\epsilon `$ is called a $`G`$-graded contraction of the algebra $`L`$ if $`L^\epsilon `$ is isomorphic to $`L`$ as vector space, $`L^\epsilon `$ has grading (1) and new commutation relations are

$$[L_j^\epsilon ,L_k^\epsilon ]_\epsilon :=\epsilon _{jk}[L_j,L_k]\epsilon _{jk}L_{j+k}^\epsilon ,$$
(2)
where the matrix $`\epsilon `$ is a solution of the equations

$$\epsilon _{jk}\epsilon _{m,j+k}=\epsilon _{km}\epsilon _{j,m+k}=\epsilon _{mj}\epsilon _{k,m+j},\epsilon _{jk}=\epsilon _{kj}.$$
(3)
The first set of equations in (3) is the consequence of Jacoby identity and the second one follows from antisimmetry of commutator.

$`G`$-grading of $`L`$-module $`V`$ is obtained from decomposition (1) and looks as follows

$$V=V_m,L_jV_mV_{j+m},$$
(4)
$$[L_j,L_k]V_m=L_jL_kV_mL_kL_jV_mV_{j+k+m}.$$
(5)
A $`G`$-graded contraction of algebra representation is defined by

$$L_j^\epsilon V_m^{\epsilon ,\psi }:=\psi _{jm}L_jV_m\psi _{jm}V_{j+m}^{\epsilon ,\psi }.$$
(6)
where the matrix $`\psi `$ satisfy the following equations

$$\epsilon _{jk}\psi _{j+k,m}=\psi _{km}\psi _{j,k+m}=\psi _{jm}\psi _{k,j+m}.$$
(7)
Hereafter we shall use for grading only cyclic group $`G=ZZ_2`$. It is convenient to write down the vector space $`V`$ and action of algebra $`L`$ for $`ZZ_2`$-grading in the form

$$V=\left(\genfrac{}{}{0pt}{}{V_0}{V_1}\right),LV=\left(\genfrac{}{}{0pt}{}{L_0L_1}{L_1L_0}\right)\left(\genfrac{}{}{0pt}{}{V_0}{V_1}\right)=\left(\genfrac{}{}{0pt}{}{L_0V_0+L_1V_1}{L_1V_0+L_0V_1}\right),$$
(8)
explicitly showing the structure of grading. Then we have for contracted module

$$L^\psi V=\left(\genfrac{}{}{0pt}{}{L_0L_1}{L_1L_0}\right)^\psi \left(\genfrac{}{}{0pt}{}{V_0}{V_1}\right)=\left(\genfrac{}{}{0pt}{}{\psi _{00}L_0\psi _{11}L_1}{\psi _{10}L_1\psi _{01}L_0}\right)\left(\genfrac{}{}{0pt}{}{V_0}{V_1}\right).$$
(9)
## 3 Grading of Virasoro algebra

The Virasoro algebra $``$(c) is the central extension of the algebra of vector fields on the circle with the basis $`c,l_n(nZZ)`$ and the commutation relations

$$[l_n,l_m]=(nm)l_{n+m}+\frac{c}{12}(n^3n)\delta _{n+m,0},$$
(10)
where $`c`$ is the central charge commuting with all generators. This algebra has the natural $`ZZ`$-grading: deg $`l_k=k`$, deg $`c=0`$ and may be graded also by cyclic groups $`ZZ_n`$. The method of $`ZZ_2`$-grading of the Virasoro algebra consists in follows. The set of the Virasoro algebra generators is divided on even $`L_0=\{A_n,c\}`$ and odd $`L_1=\{B_n\}`$ parts, where

$$A_n=\frac{1}{2}\left(l_{2n}+\frac{c}{8}\delta _{n,0}\right),B_n=\frac{1}{2}l_{2n+1},$$
(11)
then the $`Z_2`$-grading conditions are held

$$=L_0L_1,[L_0,L_0]L_0,[L_0,L_1]L_1,[L_1,L_1]L_0.$$

and the commutation relations of Virasoro algebra for new generators look as

$`[A_n,A_m]`$ $`=`$ $`(nm)A_{n+m}+{\displaystyle \frac{2c}{12}}(n^3n)\delta _{n+m,0},`$ (12)
$`[B_n,B_m]`$ $`=`$ $`(nm)A_{n+m+1}+{\displaystyle \frac{2c}{12}}(n{\displaystyle \frac{1}{2}})(n+{\displaystyle \frac{1}{2}})(n+{\displaystyle \frac{3}{2}})\delta _{n+m+1,0},`$ (13)
$`[A_n,B_m]`$ $`=`$ $`(nm{\displaystyle \frac{1}{2}})B_{n+m}.`$ (14)
The first line of equations means that the subalgebra $`\{A_n,c\}`$ is again Virasoro algebra but with central charge 2c. In the case of $`ZZ_n`$grading we have Virasoro subalgebra with central charge $`nc`$. This looks like the result of generating the orbifold algebras by $`ZZ_\lambda `$-orbifold induction procedure $`c\lambda c`$ .

## 4 Reducibility of representations. $`ZZ_2`$-grading

In this section we recall a basic notions and results of the theory of Virasoro algebra representations  and describe their $`ZZ_2`$-grading. Let $`M`$ be the module over $`(c)`$ and there is the vector $`|v`$ such that

$$l_0|v=h|v,\widehat{c}|v=c|v,l_k|v=0,k>0,$$
(15)
then $`|v`$ is a highest vector. The space spanned by vectors $`l_{i_1}\mathrm{}l_{i_m}|v,`$ where $`i_1\mathrm{}i_m>0`$, form a highest vector representation $`M(h,c)`$. All vectors of module are classified by levels. The level number is given by $`n=i_1+\mathrm{}+i_m`$. For example $`l_2l_1|v`$ belongs to the third level. Therefore the module has natural $`ZZ`$-grading $`M(h,c)=M^n(h,c)`$. The number of basis vectors on the $`n`$-th level is defined by the number $`p(n)`$ of partitions of $`n`$ on positive integer numbers: dim$`M^n(h,c)=p(n)`$. Some values of $`p(n)`$ are given below:

$$p(0)=1,p(1)=1,p(2)=2,p(3)=3,p(4)=5,p(5)=7.$$

The investigation of reducibility conditions of module is an important problem. A module is called degenerate if there is a null-vector $`|\chi _n`$, such that

$$l_0|\chi _n=(h+n)|\chi _n,l_k|\chi _n=0,k>0,$$

where $`n`$ is called the degeneracy of level. The degenerated module contains the submodule with highest vectors $`|\chi _n`$ on the $`n`$-th level and therefore is a reducible module. Null-vector is found as a linear combination of basis vectors with unknown coefficients. The conditions of null-vectors lead to the system of linear equations with two parameters $`h`$ and $`c`$. The number of equations and the number of unknowns are given by $`N_1(n)=p(n1)+p(n2)`$ and $`N_2(n)=p(n)1`$ respectively.

On the second level we have two equations and only one unknown. The relation between h and c on this level is given by

$$h=\frac{5c\pm \sqrt{(c25)(c1)}}{16}.$$

The number of equations on the next levels increases faster then the number of unknowns and therefore the existence of null-vectors becomes a rarity. The problem of searching of null-vectors is completely solved with the help of bilinear form on the $`M(h,c)`$. Let $`w`$ be an antiautomorphizm $`(c)`$ defined by $`w(l_i)=l_i`$, $`w(c)=c,`$ then there is a symmetrical bilinear form $`()`$, such that $`(a|lb)=(w(l)a|b),`$ where $`a,bM(h,c),l(c)`$. Let the $`p(n)\times p(n)`$ matrix $`K^n(h,c)`$ is the matrix of bilinear form of basis vectors on $`n`$-th level. The vanishing of determinant of matrix $`K^n(h,c)`$ is the condition of existence of null-vector . There are some possibilities depending on the values of parameters $`h,c`$: 1) the module is irreducible, then there are no null-vectors; 2) the submodules generated by null-vectors are embedded one into another; 3) there are two submodules containing all other submodules. In the last case the modules $`M(h,c)`$ are defined by

$$c=1\frac{6(qp)^2}{pq},h=\frac{(qrps)^2(qp)^2}{4pq},$$

where $`0<r<p,0<s<q`$, for integers $`r,s`$ and $`q=p+1=3,4,\mathrm{}`$. The last modules describe the space of fields of minimal models in conformal fields theory. Degeneracy of module leads to the differential equations on the fields correlators .

We are coming now to the describing of $`ZZ_2`$-grading of module $`M(h,c)`$. We shall consider that the Virasoro algebra is $`ZZ_2`$-graded by the above described method. The representation space may be divided on two subspaces

$$M(h,c)=M_0(h,c)+M_1(h,c),$$

$$M_0(h,c)=M^{2n}(h,c),M_1(h,c)=M^{2n+1}(h,c).$$

The grading of elements $`l_{i_1}\mathrm{}l_{i_m}`$ of enveloping algebra is realized by decomposition on two subset with odd and even value of sum $`i_1+\mathrm{}+i_m`$. The vectors from $`M_0(h,c)`$ and $`M_1(h,c)`$ in terms of generators $`A`$ and $`B`$ look as

$$AA\mathrm{}AA\underset{k}{\underset{}{BB\mathrm{}BB}}|v,$$
(16)
where $`k`$ is even for $`M_0(h,c)`$ and odd for $`M_1(h,c)`$. Then the following equations of $`ZZ_2`$-grading are realized

$$AM_0M_0,AM_1M_1,BM_0M_1,BM_1M_0.$$

The structure of the $`ZZ_2`$-grading module may be written as

$$LM=\left(\genfrac{}{}{0pt}{}{AB}{BA}\right)\left(\genfrac{}{}{0pt}{}{M_0}{M_1}\right)=\left(\genfrac{}{}{0pt}{}{AM_0+BM_1}{BM_0+AM_1}\right).$$
(17)
Let us note that antiautomorphism $`w`$ in $`ZZ_2`$-graded notation looks as

$`w(A_i)=A_i,w(B)=B_{i+1},w(c)=c.`$

## 5 $`ZZ_2`$-graded contraction of Virasoro algebra

Consider the contraction given by solution $`\epsilon ^\alpha =\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)`$ of equations (3). The commutation relations (14) take the form

$`[A_n,A_m]`$ $`=`$ $`(nm)A_{n+m}+{\displaystyle \frac{2c}{12}}(n^3n)\delta _{n+m,0},`$
$`[A_n,B_m]`$ $`=`$ $`(nm{\displaystyle \frac{1}{2}})B_{n+m},[B_n,B_m]=0.`$
The contracted algebra has the structure of semidirect sum of Virasoro algebra with double central charge $`(2c)=\{A_n,2c\}`$ and infinite dimensional abelian one $`\{B_n\}`$. It may be called as inhomogeneous Virasoro algebra. The second solution of equations (3) $`\epsilon ^\beta =\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)`$ leads to the contracted Virasoro algebra, which is a direct sum of $`(2c)`$ and abelian $`\{B_n\}`$

$`[A_n,A_m]`$ $`=`$ $`(nm)A_{n+m}+{\displaystyle \frac{2c}{12}}n(n^21)\delta _{n+m,0},`$
$`[A_n,B_m]`$ $`=`$ $`0,[B_n,B_m]=0.`$
## 6 $`ZZ_2`$-Grading contractions of Virasoro modules

It is easily to understand the structure of contracted $`ZZ_2`$-module using the relations (17) and (9). We take the solution of equations (7) $`\psi ^\alpha =\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right).`$ In this case

$$LM=\left(\genfrac{}{}{0pt}{}{AB}{BA}\right)\left(\genfrac{}{}{0pt}{}{M_0}{M_1}\right)=\left(\genfrac{}{}{0pt}{}{AM_0}{BM_0+AM_1}\right).$$
(18)
Therefore all properties of contracted module may be obtained by assuming

$$BM_1=0.$$
(19)
which means that all vectors of the form (16) with $`k2`$ are vanish. In other words the contracted module is obtained from the initial one by factorization with respect to condition $`BB|v=0`$.

It follows from (11) that the parameters of contracted $`M_\psi (h,c)`$ and initial $`M(h_0,c_0)`$ modules are connected by relations $`c=2c_0,h=\frac{1}{2}(h_0+\frac{c_0}{8}).`$ It is easily to understand that the structure of vectors on the even level $`2n`$ of contracted module and the structure of vectors on the level $`n`$ of initial one are identical. Therefore dim$`M_\psi ^{2n}(h,c)=\mathrm{dim}M^n(h,c)=p(n)`$. The vectors of odd level look as $`A_{i_1}\mathrm{}A_{i_n}B_q|v`$. The dimension of odd level is given by

$$\mathrm{dim}M_\psi ^{2n+1}(h,c)=\underset{k=0}{\overset{n}{}}p(k).$$
(20)
Phisical requirement that the spectrum must be bounded from below leads to the highest vector condition (15). Since any element $`l_k`$ is expressed as multiple commutators of $`l_1`$ and $`l_2`$, then in (14) instead of $`l_k|v=0,k>0`$ it is enough to impose two conditions: $`l_1|v=0,l_2|v=0`$. Then the highest vector conditions (15) are simplified

$$A_0|v=h|v,\widehat{c}|v=c|v,B_0|v=0,A_1|v=0.$$
(21)
For $`\epsilon ^\alpha `$-contracted Virasoro algebra these conditions needs to be modified by adding the equation $`A_2|v=0,`$ since otherwise the generators $`A_k,k>1`$ are not generated by commutators of $`B_0`$ and $`A_1`$. Thus the existence conditions of null-vector on the level $`n`$ are

$$A_0|\chi =(h+n)|\chi ,B_0|\chi =0,A_1|\chi =0,A_2|\chi =0.$$
(22)
The increase of number of equations on even levels with compared to (21) leads to vanishing of the null-vectors on these levels. It is possible to verify this fact by direct calculations. The relations $`B_0|\chi =0`$ for odd level are held due to (19). Thus the number of conditions remains the same therefore one would expect the presence of null-vectors on odd levels. The number $`N_1`$ of equations and the number $`N_2`$ of unknown coefficients, which defined the existence of the null-vectors on $`(2n+1)`$-th level are easily obtained from (20) and (22)

$$N_1=\underset{k=0}{\overset{n1}{}}p(k)+\underset{k=0}{\overset{n2}{}}p(k),N_2=\underset{k=0}{\overset{n}{}}p(k)1$$

It follows from this equations that $`N_1=N_2`$ up to eleventh level and $`N_1>N_2`$ for higher-order levels.

Let us investigate some first levels in details. One find by direct calculations that there are null-vectors on the first and on the third levels

$`|\chi _1=B_1|v,|\chi _3=\left(B_2\frac{5}{4h+2}A_1B_1\right)|v,`$

for any values of $`h,c`$. On the fifth level the null-vector looks as

$`|\chi _5=\left(B_3+\alpha A_2B_1+\beta A_1B_2+\gamma A_1^2B_1\right)|v,`$

where $`\alpha ,\beta ,\gamma `$ are solutions of system with parameters $`h`$ and $`c`$

$`(4h+6)\beta =7,`$

$`(4+8h+c)\alpha +15\beta +(12h+6)\gamma =9,`$

$`6\alpha +5\beta +(8h+8)\gamma =0.`$

It is immediately follows from the first equation that there are no null-vectors on the straight line $`h=3/2`$. The substitution of the first equation into the others reduces the analysis of system to possible disposition of a straight lines on the plane $`(\alpha ,\gamma )`$. Let $`D(h,c)`$ be the determinant of matrix composed from the coefficients of the left part of system and let $`D_i(h,c)`$, $`i=1,2`$ be the determinant of the matrix, which is obtained from the previous one by the replacement of $`i`$-th column with the column of the right part of the system. There are three possibility depending on the parameters $`h`$ and $`c`$: 1) the straight lines are crossed – in this case there are null-vectors; 2) the straight lines are parallel – in this case there are not null-vectors; 3) the straight lines are coincided – there are null-vectors, moreover the null-vectors become two-dimensional or ”null-plane”. The second case is realized on the straight line $`h=\frac{3}{2},`$ and two curves $`D(h,c)=0`$ on (h,c)-plane, i.e.

$$h=\frac{3c\pm \sqrt{(c25)(c1)}}{16},$$

The third case is a possibile at two points on (h,c)-plane: $`A(h_1=\frac{3}{2},c_1=26)`$ and $`B(h_2=\frac{11}{24},c_2=\frac{184}{105})`$ which are the intersection of the curve $`D(h,c)=0`$ and the straight line $`D_2(h,c)=0`$, i.e. $`c=(176496h)/35.`$ However the point $`A`$ is on the straight line $`h=\frac{3}{2}`$, where there are not null-vectors. So the null-plane on the fifth level appears at the unique values of parameters corresponding to the point $`B(h_2,c_2)`$. Similar result was obtained in paper  where it was shown that certain Verma modules over the $`N=1`$ Ramond algebra contain degenerate two-dimensional singular vector spaces. For the rest values of parameters $`h`$ and $`c`$ the first case is realized.

On the seventh and ninth levels the numbers of equations and unknown coefficients are coincided and are equal to 6 and 11 respectively. However the calculations become too cumbersome, although it is clear that null-vectors are absent only in a special cases. Starting with eleventh level the number of equations begins to exceed the number of unknowns $`N_1(11)=19,N_2(11)=18`$ and the situation becomes apparently similar to the case of noncontracted module. The attempt of using the ordinary bilinear form as in the case of noncontracted algebra does not leads to the solution of the problem of searching of null-vectors. The point is that all vectors on odd levels are orthogonal to each other due to (19). It is an ordinary thing for a finite Lie algebra that the Killing form is degenerated under contraction. At this point it may be useful the results of the work , where the way of construction of nondegenerate forms for the contracted algebras was proposed.


format: text

text:
# Computer simulations of interferometric imaging with the VLT interferometer and its AMBER instrument

## 1 INTRODUCTION

The Very Large Telescope Interferometer (VLTI) of the European Southern Observatory (ESO) with its four 8.2 m unit telescopes (UTs) and three 1.8 m auxiliary telescopes (ATs) will certainly establish a new era of optical and infrared interferometric imaging within the next few years. With a maximum baseline of up to more than 200 m, the VLTI will allow the study of astrophysical key objects with unprecendented resolution opening up new vistas to a better understanding of their physics.

The near-infrared focal plane instrument of the VLTI, the Astronomical MultiBEam Recombiner (AMBER), will operate between 1 and 2.5 $`\mu `$m and for up to three beams allowing the measurement of closure phases. In a second stage its wavelength coverage is planned to be extended to 0.6 $`\mu `$m. Objects as faint as $`K=20`$ mag are expected to be observable with AMBER when a bright reference star is available, and as faint as $`K=14`$ mag otherwise.

Among the astrophysical key issues are, for instance, young stellar objects, active galactic nuclei and stars in late stages of stellar evolution. The simulation of interferometric imaging of a stellar object consists in principle of two components:
(i)   the calculation of an astrophysical model of the object, typically based on radiative transfer calculations predicting, e.g., its intensity distribution. To obtain a robust and non-ambiguous model, it is of particular importance to take diverse observational constraints into account, for instance the spectral energy distribution and visibilities.
(ii)   the determination of the interferometer’s response to this intensity signal, i.e. the simulation of light propagation in the atmosphere and the interferometer.
Often, only one of the above parts is considered in full detail. The aim of this study is to combine both efforts and to present a computer simulation of the VLTI performance for observations of one object class, the dusty supergiants. For this purpose, we calculated a detailed radiative transfer model for one of its most outstanding representatives, the supergiant IRC +10 420, and carried out computer simulations of VLTI visibility measurements. The goal is to estimate how accurate visibilities can be measured with the VLTI in this particular but not untypical case, to discuss if the accuracy is sufficient to distinguish between different theoretical model predictions and to study on what the accuracy is dependent.

## 2 The supergiant IRC +10 420: Evolution on human timescales

The star IRC +10 420 (= V 1302 Aql) is an outstanding object for the study of stellar evolution. Its spectral type changed from F8 I$`{}_{}{}^{+}{}_{\mathrm{a}}{}^{}`$ in 1973 to mid-A today corresponding to an increase of its effective temperature of 1000-2000 K within only 25 yr. It is heavily obscured by circumstellar dust due to strong mass loss with rates typically of the order of several $`10^4`$ M/yr (M: solar mass) . IRC +10 420 is believed to be a massive luminous star of initially $`20`$ to 40 M currently being observed in its rapid transition from the red supergiant stage to the Wolf-Rayet ddphase. Wolf-Rayet stars, in turn, finally evolve into a supernova explosion. IRC +10 420 is the only object observed until now in its transition phase to the Wolf-Rayet stage.

Several infrared speckle and coronographic observations  were conducted to study the dust shell of IRC +10 420. The most recent study of Blöcker et al.  reports the first diffraction-limited 73 mas (mas: milli-arcsecond) bispectrum speckle interferometry of IRC +10 420 and presents the first radiative transfer calculations that model both the spectral energy distribution and the visibility of this key object. We will briefly describe the main results and conclusions of this study which will serve as astrophysical input for the VLTI computer simulations presented in the next section.

Fig. 1 shows the SED and the reconstructed 2.11 $`\mu `$m visibility function of IRC +10 420. The visibility 0.6 at frequencies $`>4`$ cycles/arcsec indicates, for instance, that the stellar contribution to the total flux is $``$ 60% and the dust shell contribution is $``$ 40%.

An extensive grid of radiative transfer models was calculated for the dust shell of IRC +10 420 assuming spherical symmetry and considering black bodies and model atmospheres as central sources of radiation, different silicates, grain-size and density distributions, various dust temperatures at the shell’s inner boundary (determining the radius of the shell’s inner boundary) and optical depths. We refer to Blöcker et al.  for a full description of the model grid. It turned out that the observed dust-shell properties cannot be matched by single-shell models but require multiple components with different density distributions. The best model was found for a dust shell with a dust temperature of $`T_1=1000`$ K at its inner radius of $`r_1=69R_{}`$ ($`R_{}`$: stellar radius). At a distance of $`r_{\mathrm{sw}}=308R_{}`$ ($`Y=r_{\mathrm{sw}}/r_1=4.5`$) the density was enhanced by a factor of $`S=40`$ and its slope within the shell changed from $`1/r^2`$ to $`1/r^{1.7}`$. The corresponding fits for SED and 2.11 $`\mu `$m visibility are shown in Fig. 1. The shell’s model intensity distribution is shown in Fig. 1 and was found to be ring-like due to a limb-brightened dust-condensation zone. The ring diameter is equal to the inner diameter of the hot shell ($`69`$ mas), and the diameter of the central star amounts to $`1`$ mas.

This two-component model can be interpreted in terms of a termination of an enhanced mass-loss phase roughly 90 yr ago. The assumption that IRC +10 420 had passed through a superwind phase in its recent history is in line with its evolutionary status of an object in transition from the Red-Supergiant to the Wolf-Rayet phase.

## 3 Interferometry with the VLTI and the AMBER instrument

In the previous section, the spatial intensity profile of the dusty supergiant IRC +10 420 (see Fig. 1) was derived by means of radiative transfer models and their comparison with photometric and interferometric observations. This $`2.11\mu `$m intensity profile will serve as object intensity profile in the simulation of monochromatic VLTI observations. The next steps are the simulation of light propagation from the object to the detector (through atmosphere, telescopes, and the AMBER wide-field mode instrument), simulation of photon noise and detector read-out noise, and finally data processing of the interferograms.

### 3.1 Computer simulation of interferometric imaging

Fig. 2 shows a flow chart of our simulation of interferometric imaging with the VLTI (ATs, or UTs with adaptive optics) and the AMBER camera in the wide-field mode (i.e. without fiber optics spatial filtering). The total number of detectable photons $`N_{\mathrm{tot}}`$ is given by the object brightness, the wavelength and band-width, the collecting area of the telescope, the total transmission of the optics, and the quantum efficiency of the detector. Since we assume that the observations are carried out with a spectral resolution of $`R=\lambda /\mathrm{\Delta }\lambda =70`$, the number of photons available in one spectral channel for observations in the $`K`$-band is given by $`N_{\mathrm{channel}}=N_{K,\mathrm{tot}}R_K/R`$, with $`R_K=2.2\mu \mathrm{m}/0.4\mu \mathrm{m}=5.5`$ being the spectral resolution of the $`K`$-band (i.e. $`R/R_K13`$ is the number of spectral channels within the $`K`$-band; $`N_{K,\mathrm{tot}}`$ is the total number of $`K`$-band photons). The simulations described below refer to observations in one spectral channel (see also Table 1). The signal-to-noise ratio of the reconstruction can be increased by averaging the reconstructions of all $`R/R_K`$ channels.

In a first step, an array of Gaussian distributed random numbers is generated and convolved with the correlation function of the atmospheric refraction index variations (see Roddier) in order to generate wavefronts degraded by atmospheric turbulence. The typical size of the atmospheric turbulence cells is given by the Fried parameter $`r_0`$. After the simulation of the entrance pupil, the next step incorporates the tip-tilt correction of the wavefront over each subpupil, but we allow for a residual tip-tilt error $`\delta _{\mathrm{tt}}`$. In the next step a typical Michelson output pupil is simulated by shifting both subpupils from the entrance pupil position, with a subpupil separation equal to the baseline, to the output pupil position where both subpupils are only separated by $``$ one subpupil diameter (pupil reconfiguration). The output pupil is chosen such that (i) in the optical transfer function the off-axis peaks are separated from the central peak, and (ii) the interferograms are sampled with the smallest number of pixels to assure the lowest influence of detector noise. The following step includes light propagation through the beam combiner lens to the focal plane. The squared modulus of the Fourier transform of the complex amplitude in front of the beam combiner lens yields the intensity distribution of a Michelson interferogram of a point source. In the next step the required object intensity distribution is simulated (given here by the intensity distribution of IRC +10 420) to obtain the Michelson interferogram of the object: The Fourier transformation of the object intensity distribution, calculated at those spatial frequencies covered by the the simulated interferometer baseline vector, is multiplied with the off-axis peaks of the transfer function of the generated Michelson interferogram. Finally, Poisson photon-noise and detector read-out noise is injected to the interferograms. The noise level depends, among other parameters, on the number of detectable photons, the total optical throughput of the interferometer (VLTI + AMBER in the wide-field mode) and the quantum efficiency of the detector. Details of the simulation parameters specific to VLTI-AT/AMBER, as optical throughput, detector read-out noise and quantum efficiency, are given in Table 1.

### 3.2 Computer simulations of visibility measurements

We performed simulations of IRC +10 420 visibility observations in the VLTI-AT/AMBER wide-field mode and studied the influence of various observational parameters on the visibility accuracy. Visibility error bars were obtained for the following observational parameters:

1. different seeing during the observation of object and reference star (Fried parameters $`r_{0,\mathrm{object}}`$ and $`r_{0,\mathrm{ref}.}`$ ranging between 0.9 m and 1.2 m),
2. different residual tip-tilt error ($`\delta _{\mathrm{tt},\mathrm{object}}`$ and $`\delta _{\mathrm{tt},\mathrm{ref}.}`$ ranging between 0.1% and 20% of the Airy disk diameter),
3. different object brightness ($`K_{\mathrm{object}}`$=0.7 mag to 10.2 mag, $`K_{\mathrm{ref}.}`$=0.7 mag).

All computer experiments are based on 2400 interferograms of the $`S=40`$ intensity profile (see Fig. 1) of IRC +10 420 and refer to observations in one of the $`R/R_K13`$ spectral channels. The error bars correspond to the reduction of six statistically independent data sets and refer to the standard deviation $`\sigma `$. The dependence of the error bar on the number of data sets has been verified (see below). In the computer experiments (1) and (2) mentioned above, object and reference star were assumed to have the same brightness. We chose $`K`$=0.7 mag (see Table 1) in order to minimize brightness effects. Experiment (3) simulates IRC +10 420 itself ($`K`$=3.5 mag) and fainter objects, assuming object and reference star to have the same seeing and residual tip-tilt error ($`r_{0,\mathrm{object}}`$=$`r_{0,\mathrm{ref}}`$=0.9 m; $`\delta _{\mathrm{tt},\mathrm{object}}`$=$`\delta _{\mathrm{tt},\mathrm{ref}}`$=0.1%). These simulations supersede the corresponding results of a first exploratory study  due to improved numerics and a broader statistical basis.

The first computer experiment (different seeing conditions) consists of the simulations A to F (see top panel of Fig. 3). Simulation A represents the case of typical seeing conditions ($`r_{0,\mathrm{object}}`$=0.9 m, $`r_{0,\mathrm{ref}.}`$=0.9 m) and an almost perfect tip-tilt correction (residual tip-tilt error $`\delta _{\mathrm{tt}}=0.1`$%). It will serve as a reference for the following simulations. The brightness of object and reference star was chosen to be $`K`$=0.7 mag. The visibility error $`\sigma _V`$ amounts only to $`\pm 0.0067`$. The remaining statistical uncertainty due to the limited number of data sets has been checked by a test calculation taking into account twice as many statistically independent data sets, i.e. 12 data sets instead of 6. The visibility error changed by 24% from $`\pm 0.0067`$ to $`\pm 0.0051`$. However, the general dependencies of the visibility accuracy on observational parameters (seeing, residual tip-tilt error, brightness) are only scarcely affected.

If one considers a better seeing for both object and reference star (A-C), $`\sigma _V`$ decreases in almost linear proportion with increasing diameter of the atmospheric turbulence cells (Fig 3). However, improved but different seeing conditions (D-F) lead to larger visibility errors. This emphasizes that equal seeing conditions for object and reference star are more crucial than excellent seeing conditions of, e.g., the object only in order to obtain accurate visibilities.

The second set of simulations (A, G-J) illustrates the influence of a larger residual object tip-tilt error (see middle panel of Fig. 3). Increasing $`\delta _{\mathrm{tt}}`$ from 0.1% to 20% for both object and reference star increases the visibility error from 0.0067 to 0.0084 (A,F,G). Like in the simulations of different seeing conditions, it is essential to have the same tip-tilt error for both object and reference star in order to achieve high accuracy. Differences in $`\delta _{\mathrm{tt}}`$ enhance the visibility error (H,I).

Finally, in the third experiment of this series an IRC +10 420-like intensity distribution is assumed but much much fainter (i.e. fainter than $`K`$=3.5 mag) objects are considered (see bottom panel of Fig. 3, simulations A, K-N). The visibility error stays almost constant up to an $`K`$ magnitude of 7 mag. If the simulated object $`K`$-magnitude is 9.2 mag, the visibility error has increased by a factor of $`2`$, i.e. to $`\sigma _V=0.0156`$, but is still acceptable although the photon number decreased to $`N_{\mathrm{channel}}=1.7610^3`$ (i.e. by a factor of 2500, see Table 1). Fig. 3 demonstrates that for even fainter objects noise becomes more and more important leading to an steep increase of the visibility error. However, the signal-to-noise ratio can be improved ($`\sqrt{R/R_K}3.6`$) if the visibility measurements of all $`R/R_K`$ spectral channels within the $`K`$-band are averaged.

For illustration, Fig. 4 shows the results of the simulations A, G, J and D for baselines of 50 and 100 m together with the model predictions for IRC +10 420. The asterisk indicates that within these simulations the object brightness is that of IRC +10 420, viz. $`K`$=3.5 mag (photon number: $`3.4710^5`$) As it is obvious from Fig. 3 (bottom panel), this brightness is well within the range where the visibility error is not very sensitive to brightness changes. Accordingly, the results are very close to those of the simulations A, G, J, and D. Again, different seeing conditions (simulation D) for object and reference star are more crucial than, e.g., residual tip-tilt errors (simulations G and J). For typical seeing conditions (Fried parameter differences $`\begin{array}{c}<\\ \end{array}`$ 10%), the visibility accuracy achievable with the VLTI using the ATs in the wide-field mode is clearly sufficient to distinguish between different radiative transfer models of IRC +10 420 and, thus, to prove (or disprove) theoretical predictions.

### 3.3 Computer simulations of phase-closure imaging

In addition to the above visibility studies, image reconstruction simulations were performed as well. For this purpose, modulus (=visibility) and phase of the Fourier transform of the object’s intensity distribution were determined utilizing the phase-closure method. The phase closure method requires the simultaneous measurement of the object with at least three telescopes. For a corresponding computer experiment we chose three telescopes covering five configurations. The baseline between telescope 1 and telescope 2 was 8 m whereas the baseline between telescope 2 and telescope 3 amounted to 8, 16, 24, 32 and 40 m in order to facilitate a simple recursive algorithm for phase reconstruction from the measured closure phases. A non-redundant arrangement of the output pupil ensures a sufficient separation of the object information in frequency space.

As model object served a binary star with an intensity distribution as shown in Fig. 5 (bottom; dashed line). The components’ intensity ratio is 1:2. As in the case of the visibility simulations of IRC +10 420 each simulation is based on 6 statistically independent data sets. The number of interferograms per data set is 1200.

Fig. 5 shows the simulated modulus, phase and image reconstruction of the binary for total object $`K`$-magnitudes of 0.7 mag and 9.2 mag, resp. Even for objects as faint as $`K=9.2`$ mag the agreement between test object and reconstructed object is very good. Beyond $`K=10`$ mag the errors increase considerably. This is also illustrated in Fig. 6 which shows the photometry error (i.e. the deviation from the components’ intensity ratio of 1:2) as a function of the object brightness. The signal-to-noise ratio of the reconstruction can be increased if the reconstructions of all $`R/R_K`$ spectral channels are averaged.

## 4 Conclusions

We have presented computer simulations of interferometric imaging with the VLT interferometer and the AMBER instrument in the wide-field mode. These simulations include both the astrophysical modelling of a stellar object by radiative transfer calculations and the simulation of light propagation from the object to the detector and simulation of photon noise and detector read-out noise. We focussed on stars in late stages of stellar evolution and examplarily studied one of its most outstanding representatives, the dusty supergiant IRC +10 420. The model intensity distribution of this key object, obtained by radiative transfer calculations, served as astrophysical input for the VLTI/AMBER simulations. The results of these simulations show the dependence of the visibility error bar on various observational parameters. With these simulations at hand one can immediately see under which conditions the visibility data quality would allow us to discriminate between different model assumptions (e.g. the size of the superwind amplitude $`S`$). Inspection of Fig. 4 shows that in all studied cases the observations will give clear preference to one particular model. Therefore, observations with VLTI will certainly be well suited to gain deeper insight into the physics of dusty supergiants.


text:
# Bianchi VIII Empty Futures

## 1 Introduction

Many studies have been made of the asymptotic states of spatially homogeneous cosmological models. Studies of the early-time behaviour have focussed upon the possibility of chaotic behaviour in the most general cases , , , while investigations of the late-time asymptotics have focussed upon the relative generality of cosmological models which evolve to become at least as isotropic as the observed universe . In this paper we are going to consider some aspects of the late-time evolution of one of the four most general classes of spatially homogeneous universe, see , , , , and  for related studies. The Bianchi type $`VIII`$ universe has not attracted the same interest from cosmologists as some other homogeneous models because it does not contain an isotropic Friedmann universe as a special case. However, it can become arbitrarily close to isotropy. It is primarily of interest because of its status as one of the most general homogeneous universes and hence in some sense, yet to be made precise, may offer a leading order approximation to some part of the general cosmological solution of Einstein’s equations at late times, far from the initial singularity. In fact, if space is topologically compact it becomes the most general of the ever-expanding Bianchi type universes .

In what follows we will study certain aspects of the late-time behaviour of vacuum Bianchi type $`VIII`$ models. In particular, we will try to determine whether the chaotic Mixmaster behaviour which takes place close to the initial singularity still takes place in the far future as $`t+\mathrm{}`$. We will show how a qualitative analysis based on the Hamiltonian formalism can reveal some basic features of the future asymptote of vacuum Bianchi $`VIII`$ models. These provide a basis for studying the nature of the late-time behaviour of the general solution to the vacuum Einstein field equations.

It is also of interest to note that late-time cosmological behaviour in general relativity corresponds to early times in the so-called Pre Big Bang scenario suggested by symmetries of superstring theory. It has been argued  that plane-wave space-times are the generic initial states for these models. So in this context the study of the asymptotic future behaviour of general vacuum space-times is also of interest. Also, recent work  has pointed out that inhomogeneous models of low-energy bosonic field equations for all superstring models and M-theory ( $`D=11`$ supergravity ) exhibit an oscillatory approach to the initial singularity, so it would be interesting as well to study how these oscillatory singularities transform in type $`VIII`$ universes at late times using a similar qualitative analysis, see also .

We will first introduce the field equations for the vacuum Bianchi type $`VIII`$ model and then discuss the orthonormal frame approach and the Hamiltonian formalism to analyse the possible future asymptotic states.

## 2 The field equations for Bianchi type $`VIII`$

In the metric approach, where the metric components are the basic variables of the gravitational field, we can introduce group invariant and time-independent frame vectors $`e_a`$ such that the line element is given by

$$ds^2=dt^2+g_{ab}W^aW^b$$

where $`W^a`$ are time-independent 1-forms dual to the frame vectors $`e_a`$. In the case of vacuum Bianchi type $`VIII`$ we can introduce three time-dependent expansion scale factors as

$$g_{ab}=diag(a^2,b^2,c^2)$$

and Einstein’s field equations become

$$\frac{(\dot{a}bc)^.}{abc}=\frac{1}{2a^2b^2c^2}[(b^2+c^2)^2a^4]$$
(1)
$$\frac{(a\dot{b}c)^.}{abc}=\frac{1}{2a^2b^2c^2}[(a^2+c^2)^2b^4]$$
(2)
$$\frac{(ab\dot{c})^.}{abc}=\frac{1}{2a^2b^2c^2}[(a^2b^2)^2c^4]$$
(3)
$$\frac{\ddot{a}}{a}+\frac{\ddot{b}}{b}+\frac{\ddot{c}}{c}=0$$
(4)
Here, the dot denotes differentiation with respect to the synchronous time $`t`$. If we replace the synchronous time $`t`$ by a variable $`\eta `$ defined by

$$dt=(abc)d\eta $$

then the field equations (1)-(4) become

$$2(\mathrm{ln}a)^{^{\prime \prime }}=(b^2+c^2)^2a^4$$
(5)
$$2(\mathrm{ln}b)^{^{\prime \prime }}=(a^2+c^2)^2b^4$$
(6)
$$2(\mathrm{ln}c)^{^{\prime \prime }}=(a^2b^2)^2c^4$$
(7)
$$(\mathrm{ln}a^2)^{^{\prime \prime }}+(\mathrm{ln}b^2)^{^{\prime \prime }}+(\mathrm{ln}c^2)^{^{\prime \prime }}=(\mathrm{ln}a^2)^{}(\mathrm{ln}b^2)^{}+(\mathrm{ln}b^2)^{}(\mathrm{ln}c^2)^{}+(\mathrm{ln}c^2)^{}(\mathrm{ln}a^2)^{}$$
(8)
where ’ denotes differentiation with respect to $`\eta `$.

These equations are invariant under the transformation$`\eta \eta `$. Now, as is widely exploited in the literature, , , to a good approximation we have $`abct`$ as $`t0,`$ which implies

$$\eta \mathrm{ln}(t)$$

and the invariance of the field equations (5)-(8) under$`\eta \eta `$ corresponds to invariance under

$$t\frac{1}{t}$$

in (1)-(4): an apparent duality transformation relating early and late-time behaviour. If this symmetry were exact it would imply that the presence of chaotic behaviour at early times ($`t0`$) means that there must also be chaotic behaviour as $`t+\mathrm{},`$ as suggested in ref. . A crucial point in the above argument is the accuracy of the assumption that

$$abct$$

One can use the field equations (1)-(3) to test this assumption, by looking at the second derivative of the metric volume element $`abc`$. Taking the sum of the equations (1)-(3) we get

$$2(abc)(abc)^{..}=(a^2b^2)^2+2a^2c^2+2c^2b^2+c^4.$$
(9)
Since the right-hand side is the sum of four positive terms, if $`abc`$ $`0`$ we have

$$(abc)^{..}>0$$
(10)
This shows that the relation $`abct`$ cannot be exact and so the duality symmetry argument cannot be applied to relate late and early time behaviour for the vacuum Bianchi type $`VIII`$ model. However, the inequality (10) will be useful in deducing the late-time behaviour as $`t+\mathrm{}`$. There is another consequence of (10) worth pointing out : if $`abc`$ is a power law

$$abct^n$$

then (10) implies $`n>1`$. Since $`dt=(abc)d\eta `$, we would have

$$\eta \frac{1}{1n}t^{1n},n1$$

Note that near $`n=1`$ there is a structural instability with respect to the variation of $`t`$ with respect to $`\eta `$. When $`n>1`$ we have $`t0`$ as $`\eta \mathrm{}`$ and $`t+\mathrm{}`$ when $`\eta 0`$. However, this differs form the case where $`n=1,`$ where $`\eta \mathrm{ln}(t)`$ (which is a good approximation close to the initial singularity) since here $`\eta +\mathrm{}`$ as $`t+\mathrm{}`$. If $`n`$ $`<1`$, then we would have $`\eta 0`$ as $`t0`$ and $`\eta +\mathrm{}`$ as $`t+\mathrm{}.`$ Tiny variations on either side of $`n=1`$ would make $`\eta `$ a bad choice of variable for an early-time analysis of Mixmaster behaviour. This may be of importance for numerical studies of Mixmaster behaviour , which are notoriously difficult to carry out because of the slow time variation of the oscillations compared with the growth in volume.

Let $`u^a`$ be the 4-velocity of the Hubble flow, the Hubble parameter is defined by

$$H=\frac{1}{3}u_{}^{a}{}_{;a}{}^{}$$

This parameter can be expressed as follows

$$H=\frac{\dot{l}}{l}$$

where $`l(t)`$ corresponds to a geometric-mean expansion scale factor. The deceleration parameter $`q`$ is defined by

$$q=\frac{l\ddot{l}}{\dot{l}^2}$$

The Hubble parameter then satisfies the following equation

$$\frac{dH}{dt}=(1+q)H^2$$
(11)
Defining the a dimensionless shear parameter by

$$\mathrm{\Sigma }^2=\frac{\sigma ^2}{3H^2}$$

with $`\sigma ^2=\frac{1}{2}\sigma _{ab}\sigma ^{ab}`$, where $`\sigma _{ab}`$ is the shear tensor, we have the following simple relation between the shear parameter and the deceleration parameter in vacuum:

$$q=2\mathrm{\Sigma }^2.$$
(12)
Equations (11) and (12) can now be used to express the Hubble parameter as

$$H(t)=\frac{1}{t+2\mathrm{\Sigma }^2(t)𝑑t}$$
(13)
For the isotropic ($`\mathrm{\Sigma }=0`$) Milne universe we have $`H(t)=t^1`$ so equation (13) expresses the departure from the Milne expansion rate for vacuum homogeneous cosmologies.

We specialise to the Bianchi $`VIII`$ universe and calculate the Hubble scalar defined above for $`u^a=(1,0,0,0)`$ using the Christoffel symbols for Bianchi type $`VIII`$ models , so:

$$H(t)=\frac{1}{3}\frac{d(abc)}{dt}\frac{1}{abc}$$
(14)
Equation (13) then gives the following expression for $`abc`$ :

$$abc=\mathrm{exp}(\frac{3}{t+2s(t)}𝑑t)$$
(15)
with $`s(t)=\mathrm{\Sigma }^2𝑑t`$

From the above equation we can get $`(abc)^{..}`$, as

$$(abc)^{..}=\mathrm{exp}(\frac{3}{t+2s(t)}𝑑t)\frac{1}{(t+2s(t))^2}[93(1+2\mathrm{\Sigma }^2)]$$
(16)
For this to agree with the inequality (10) we require

$$93(1+2\mathrm{\Sigma }^2)>0$$

which implies

$$\mathrm{\Sigma }^2<1$$
(17)
As we will see this is consistent with bounds on $`\mathrm{\Sigma }`$ given by equation (28) of the next section.

## 3 The orthonormal frame formalism

We use the orthonormal frame formalism of Ellis and MacCallum . The orthonormal frame $`e_a`$ is chosen to be invariant under the action of the isometry group with $`e_0=u,`$ so that the commutation functions $`\gamma _{}^{c}{}_{ab}{}^{}`$ defined by

$$[e_a,e_b]=\gamma _{}^{c}{}_{ab}{}^{}e_c$$

depend on the synchronous time $`t`$ and are the basic variables of the gravitational field. With this choice of frame, the line-element takes the form

$$ds^2=dt^2+\delta _{ab}\omega ^a(t)\omega ^b(t),$$

where $`\omega ^a`$ are time dependent 1-forms dual to the vector fields $`e_a`$.

The physical state of the models can be described using the components of the shear tensor $`\sigma _{ab}`$ and the curvature anisotropies $`n_{ab}`$ which occur in the expression of the commutation functions, see Wainwright and Ellis  and , . In the above frame we define

$$\sigma _{ab}=diag(\sigma _{11},\sigma _{22},\sigma _{33})$$

and

$$n_{ab}=diag(n_1,n_2,n_3)$$

and the linear combinations

$$\sigma _+=\frac{1}{2}(\sigma _{22}+\sigma _{33})$$

$$\sigma _{}=\frac{1}{2\sqrt{3}}(\sigma _{22}\sigma _{33})$$

$$n_+=\frac{1}{2}(n_2+n_3)$$

$$n_{}=\frac{1}{2\sqrt{3}}(n_2n_3)$$

Now scale out the overall expansion of the universe by introducing the expansion-normalised variables

$$\mathrm{\Sigma }_\pm =\frac{\sigma _\pm }{H}$$
(18)
$$N_a=\frac{n_a}{H}$$
(19)
The cosmological model is then fully described by the dimensionless variables $`(\mathrm{\Sigma }_+,\mathrm{\Sigma }_{},N_1,N_2,N_3)`$. The variables $`\mathrm{\Sigma }_\pm `$ describe the anisotropy of the Hubble flow and the $`N_a`$ describe the spatial curvature of the group orbits. If we define a dimensionless curvature parameter by

$$K=\frac{R^{(3)}}{6H^2}$$

then we can deduce from the generalised Friedmann equation the relation

$$1=\mathrm{\Sigma }^2+K$$
(20)
Note that the curvature parameter $`K`$ is related to the variables $`N_a`$ by

$$K=\frac{1}{12}[N_{1}^{}{}_{}{}^{2}+N_{2}^{}{}_{}{}^{2}+N_{3}^{}{}_{}{}^{2}2(N_1N_2+N_2N_3+N_3N_1)]$$

We can introduce a dimensionless time variable $`\tau `$ defined by

$$\frac{dt}{d\tau }=\frac{1}{H}$$
(21)
and related to the mean expansion scale factor, $`l(t)`$, by

$$l=l_0e^\tau $$
(22)
These variables yield the following simplified evolution equations for the vacuum type $`VIII`$ state variables :

$$\mathrm{\Sigma }_{\pm }^{}{}_{}{}^{^{}}=(2q)\mathrm{\Sigma }_\pm S_\pm $$
(23)
$$N_{1}^{}{}_{}{}^{^{}}=(q4\mathrm{\Sigma }_+)N_1$$
(24)
$$N_{2}^{}{}_{}{}^{^{}}=(q+2\mathrm{\Sigma }_++2\sqrt{3}\mathrm{\Sigma }_{})N_2$$
(25)
$$N_{3}^{}{}_{}{}^{^{}}=(q+2\mathrm{\Sigma }_+2\sqrt{3}\mathrm{\Sigma }_{})N_3$$
(26)
where  denotes derivatives with respect to $`\tau `$ and the functions $`S_\pm `$ are given by :

$$S_+=\frac{1}{6}[(N_2N_3)^2N_1(2N_1N_2N_3)]$$

$$S_{}=\frac{1}{2\sqrt{3}}(N_2N_3)(N_1N_2N_3)$$

This system, together with equation (20), is an autonomous system of first-order ordinary differential equations for which a monotone function can be found with respect to the time variable $`\tau `$. Different Bianchi types lead to different restrictions on the ranges of the variables $`N_a`$. For Bianchi type $`VIII`$ we require $`N_1<0`$, $`N_2>0`$, $`N_3>0`$ and the set defined by these constraints is an invariant set for the evolution equations. On this invariant Bianchi type $`VIII`$ set there exists the following monotone increasing function :

$$Z=(N_1N_2N_3)^2$$
(27)
As a consequence of eq. (20) the following bounds can placed on the state variables  :

$$0\mathrm{\Sigma }^21$$
(28)
$$1\mathrm{\Sigma }_\pm +1$$
(29)
$$0N_{1}^{}{}_{}{}^{2}12$$
(30)
$$0N_1(N_2N_3)12$$
(31)
$$02N_1(N_2+N_3)12$$
(32)
$$0(N_2N_3)^212$$
(33)
Therefore, the invariant set $`S`$ for Bianchi type $`VIII`$ is unbounded in such way that $`N_2+N_3`$ can assume arbitrarily large values. Together with the existence of the monotone increasing function $`Z,`$ given by equation (27), this means we do not expect Bianchi type $`VIII`$ orbits to approach an equilibrium point as $`\tau +\mathrm{}`$.

## 4 Hamiltonian Cosmology

As discussed in section 2, in the metric approach, we can introduce a set of group-invariant and time-independent one-forms $`W^a`$, so that with an arbitrary time variable $`t^{}`$ the metric has the form

$$ds^2=N(t^{})dt_{}^{}{}_{}{}^{2}+g_{ab}(t^{})W^aW^b$$
(34)
In vacuum, the Bianchi type $`VIII`$ metric is diagonal and we write the metric components in the form

$$g_{ab}=diag(e^{2\beta _1},e^{2\beta _2},e^{2\beta _3})$$
(35)
with

$$\beta _1=\beta ^02\beta ^+$$

$$\beta _2=\beta ^0+\beta ^++\sqrt{3}\beta ^{}$$

$$\beta _3=\beta ^0+\beta ^+\sqrt{3}\beta ^{}$$

In the Hamiltonian approach the spatial Ricci curvature defines a potential function $`V(\beta _+,\beta _{})`$,

$$V(\beta _+,\beta _{})=e^{2\beta ^0}R^{(3)}$$

The variable $`\beta ^0`$ is related to the overall length scale for the metric (34) and to the Hubble parameter by the relations

$$l=l_0e^{\beta ^0}$$
(36)
$$H=N^1\dot{\beta ^0}$$
(37)
The evolution of the Bianchi type $`VIII`$ universe is then described by the motion of a point in the mini-superspace $`(\beta ^0,\beta ^+,\beta ^{})`$ evolving in the potential $`V(\beta ^+,\beta ^{})`$. The variable $`\beta ^0`$ can be chosen to be the time variable so that $`\dot{\beta ^0}=1`$ and $`N=H^1`$. Equations (22) and (36) imply that we have $`\beta ^0=\tau `$.

This Hamiltonian approach allows us to study the early and late-time asymptotic behaviour of Class A Bianchi-type universes qualitatively. For vacuum Bianchi type $`I`$, the potential function vanishes, $`V(\beta ^\pm )=0`$, so that the evolution is a straight line in $`\beta ^\pm `$-space. For Bianchi type $`VIII`$ the equipotential curves of $`V(\beta ^\pm )`$ are approximately triangular in shape, , with an infinite open channel which starts at the vertex on the $`\beta ^+`$ axis and becomes increasingly narrow as $`\beta _++\mathrm{}`$. As $`\tau +\mathrm{}`$ (evolution away from the initial singularity) the walls of the triangle contract. If the universe point moves initially in a straight line (corresponding to Kasner-like behaviour) it will bounce off a potential wall and then move along another straight line (a transition to another Kasner-like behaviour). The triangular shape of the potential allows the evolution to be chaotic, which corresponds to the Mixmaster behaviour and is familiar from studies of the evolution near the singularity ($`\tau \mathrm{}).`$ But, what happens as $`\tau +\mathrm{}`$ ? The Hamiltonian picture implies that the triangular potential contracts, so this would enable the universe point to bounce off the walls in a chaotic sequence, ad infinitum. We shall argue however from the inequality (10), or the bound (28), and the fact that $`N_2+N_3`$ can take arbitrarily large values, that this chaotic future evolution cannot occur as $`\tau +\mathrm{}`$, contrary to the claim of Halpern, .

In order to show this, we introduce two further equations, relating the $`\mathrm{\Sigma }_\pm `$ to the $`\beta ^\pm `$ variables. Since

$$\mathrm{\Sigma }_\pm =\frac{d\beta ^\pm }{d\tau },$$
(38)
where $`\mathrm{\Sigma }_\pm `$ coincides with expression calculated in the orthonormal frame formalism, using equation (12), we have the relation

$$\frac{q}{2}=\mathrm{\Sigma }^2=\mathrm{\Sigma }_{+}^{}{}_{}{}^{2}+\mathrm{\Sigma }_{}^{}{}_{}{}^{2}=(\frac{d\beta ^+}{d\tau })^2+(\frac{d\beta ^{}}{d\tau })^2$$
(39)
The Bianchi type $`VIII`$ model exhibits Mixmaster behaviour towards the past and this chaotic evolution is reflected in the infinite sequence of Kasner era’s and bounces off the triangular potential walls in the Hamiltonian picture , , . When bounces occur, the function

$$(\frac{d\beta ^+}{d\tau })^2+(\frac{d\beta ^{}}{d\tau })^2$$
(40)
exhibits a series of peaks (a bounce produces an abrupt change in $`\beta ^+`$ or $`\beta ^{}`$). By equation (39) this implies that $`\mathrm{\Sigma }^2`$ should exhibit a series of peaks as well. A graph of the function $`\mathrm{\Sigma }`$ is given in Wainwright and Ellis . Now, if the triangular potential contracts as $`\tau +\mathrm{}`$, the interval between bounces decreases as well. This means that time derivatives of $`\beta _\pm `$ could diverge (the peaks of $`\beta _\pm (\tau )`$ would become increasingly narrow) and so the function (40) would be unbounded. This would imply, by equation (39), that $`\mathrm{\Sigma }^2`$ would be unbounded, which contradicts the inequality (17) or (28). In this picture of evolution in an eternally contracting potential, one can imagine two solutions which avoid this problem.

First, it could be that $`\beta _\pm (\tau )`$ is of the form

$$\beta _\pm (\tau )=e^{a\tau }osc(be^{a\tau })+c,$$

with $`a`$, $`b`$, $`c`$ constants. In the above the $`osc`$ function could have a periodic or chaotic oscillatory behaviour. Clearly this could give increasingly rapid chaotic oscillations of decreasing amplitude and using equation (38) one finds a bounded, eternally chaotic, oscillating $`\mathrm{\Sigma }_\pm `$. However, defining the function $`p_0(\tau )`$ by

$$\frac{dp_0}{d\tau }=(2q)p_0$$

one can show the following relations in the Hamiltonian approach

$$\stackrel{~}{N}_a=\frac{12e^{2\beta _a}}{(p_0)}$$
(41)
where $`\stackrel{~}{N}_a=n_aN_a,`$ ( no sum ) with $`n_1=C_{}^{1}{}_{23}{}^{}`$, $`n_2=C_{}^{2}{}_{31}{}^{}`$ and $`n_3=C_{}^{3}{}_{12}{}^{}.`$ Using equations (41), (35) and (12), one finds,

$$|N_1N_2N_3|\mathrm{exp}(6\tau 6(1\mathrm{\Sigma }^2)𝑑\tau )$$
(42)
We see that if $`\mathrm{\Sigma }_\pm `$ oscillates eternally and satisfies (17) or (28), then (42) implies that $`|N_1N_2N_3|`$ increases as $`\tau +\mathrm{}`$ (see end of section 3). Now, equation (30) implies that $`N_1`$ is bounded. Since $`\beta _\pm `$ tend to constant values, the expression for $`\beta _1`$ in (35) and equation (41) with $`a=1`$ implies that

$$\frac{e^{2\beta _0}}{(p_0)}$$

is bounded, and this means that equations (35) and (41) with $`a=2`$ and $`a=3`$ will imply that $`N_2`$ and $`N_3`$ are each bounded. Hence the Bianchi type $`VIII`$ invariant set would be bounded and this clearly contradicts the fact that $`N_2+N_3`$ diverges as $`\tau +\mathrm{}`$ because of equation (42).

Second, the universe point could slow down, but then we could end up reaching a point where $`\beta _\pm `$ constant and $`\mathrm{\Sigma }^20`$. One can easily see that in this case equation (23) will imply that equation (20) will not be satisfied. Indeed, $`\mathrm{\Sigma }_{}=0`$ will imply that $`N_2=N_{3\text{ }}`$ for Bianchi $`VIII`$ , so equation (23) requires $`\mathrm{\Sigma }_{}^{}{}_{}{}^{^{}}=0`$. Equations (24),…, (26), (41) and (35) imply that the $`N_a`$ become constants, but then (23) requires $`S_+=0`$, otherwise $`\mathrm{\Sigma }_+`$ does not satisfy the bound (17) or (28). In order to have $`S_+=0`$ we must satisfy either $`N_1=0`$ or $`N_1=N_2=N_3`$. The first case implies that the curvature parameter $`K=0`$, and this is inconsistent with the Friedmann equation (20), and the second possibility implies that $`N_1=N_2=N_3=0`$ because of the sign restrictions mentioned in the previous section for Bianchi type $`VIII`$. This also leads to $`K=0`$ and contradicts (20) as well. Note that this argument is not necessarily true in the non-vacuum situation (see appendix).

So what are the alternative possibilities consistent with the Hamiltonian picture? The infinitely long channel along the $`\beta ^+`$ axis offers a solution. The universe point could evolve eternally between the two walls of the infinite channel, for instance there could be an infinity of oscillations between the two walls of the channel. This evolution would correspond to a non-chaotic oscillation of the universe point about the axisymmetric solution (see below) characterised by

$$\beta ^{}=0$$

In this case, because the channel becomes increasingly narrow as $`\beta _++\mathrm{}`$, we expect $`\beta _{}`$ to undergo increasingly rapid oscillations of decreasing amplitude as $`\tau +\mathrm{}`$. Similar behaviour of the universe point when interacting with a channel of the Bianchi IX potential can be found in Ryan and Shepley  for a so-called ’mixing bounce’. Equation (42) also implies that $`N_2+N_3+\mathrm{}`$ as $`\tau +\mathrm{}`$ and so the orbits in the state space escape to infinity.

Non-chaotic oscillatory behaviour for Bianchi type $`VIII`$ seems also to be observed in the behaviour of the Weyl curvature of other Bianchi type universes according to Coley  and in Wainwright et al. . In this interesting work the late-time behaviour of Bianchi type $`VII_0`$ models with perfect fluid source was investigated, (see also refs. , , . If the universe contains matter with pressure $`p`$ and density $`\rho ,`$ obeying an equation of state

$$p=(\gamma 1)\rho $$

with $`\gamma `$ constant, then if $`\gamma <\frac{4}{3}`$ the universe isotropizes with respect to shear but not with respect to the Weyl curvature. So even if the shear decreases at late times, the anisotropic Weyl curvature dominates the dynamics, producing a form of self-similarity breaking at late times. If $`\frac{4}{3}<\gamma 2`$, then $`\mathrm{\Sigma }^2`$ is oscillatory as $`\tau +\mathrm{},`$ with $`\beta _{}`$ a decreasing oscillatory function, while the Weyl scalar diverges. It was further argued that since an open subset of a more general class of models (i.e. Bianchi type $`VIII`$) will display features of Bianchi $`VII_0`$ models, the non-self-similar behaviour of the Weyl curvature could occur at some stage during the evolution of Bianchi $`VIII`$ models.

Let us introduce the dimensionless electric and magnetic parts of the Weyl curvature

$$\stackrel{~}{E}_{ab}=E_{ab}\frac{1}{H^2}$$

$$\stackrel{~}{H}_{ab}=H_{ab}\frac{1}{H^2}$$

For all Bianchi class A models, in the standard frame $`\stackrel{~}{E}_{ab}`$ and $`\stackrel{~}{H}_{ab}`$ are diagonal and we define:

$$\stackrel{~}{E}_+=\frac{1}{2}(\stackrel{~}{E}_{22}+\stackrel{~}{E}_{33})$$
(43)
$$\stackrel{~}{E}_{}=\frac{1}{2\sqrt{3}}(\stackrel{~}{E}_{22}\stackrel{~}{E}_{33})$$
(44)
and likewise for $`\stackrel{~}{H}_+`$ and $`\stackrel{~}{H}_{}`$. By defining

$$N_+=\frac{1}{2}(N_2+N_3)$$
(45)
$$N_{}=\frac{1}{2}(N_2N_3).$$
(46)
We have the following useful relations between the state variables introduced in the orthonormal frame approach and the electric and magnetic parts of the Weyl curvature :

$$\stackrel{~}{E}_+=\mathrm{\Sigma }_++(\mathrm{\Sigma }_{+}^{}{}_{}{}^{2}\mathrm{\Sigma }_{}^{}{}_{}{}^{2})+S_+$$
(47)
$$\stackrel{~}{E}_{}=\mathrm{\Sigma }_{}2\mathrm{\Sigma }_+\mathrm{\Sigma }_{}+S_{}$$
(48)
$$\stackrel{~}{H}_+=\frac{3}{2}N_1\mathrm{\Sigma }_+3N_{}\mathrm{\Sigma }_{}$$
(49)
$$\stackrel{~}{H}_{}=3N_{}\mathrm{\Sigma }_++\frac{1}{2}(N_14N_+)\mathrm{\Sigma }_{}.$$
(50)
If $`N_2+N_3`$ becomes arbitrarily large, then by equation (50), we can have

$$\underset{\tau +\mathrm{}}{lim}\stackrel{~}{H}_{}=+\mathrm{}$$
(51)
and the Weyl scalar, given by

$$W=\stackrel{~}{E}_{+}^{}{}_{}{}^{2}+\stackrel{~}{E}_{}^{}{}_{}{}^{2}+\stackrel{~}{H}_{+}^{}{}_{}{}^{2}+\stackrel{~}{H}_{}^{}{}_{}{}^{2},$$

can become arbitrarily large as $`\tau +\mathrm{}`$, see Wainwright  for a discussion of the self-similarity breaking. Similar behaviour for the Weyl curvature was found in perfect fluid Bianchi type $`VII_0`$ universes with equation of state parameter satisfying $`1<\gamma 2`$ .

## 5  The axisymmetric (LRS) vacuum Bianchi type $`VIII`$ model

Our qualitative analysis shows that the late-time behaviour of vacuum Bianchi type $`VIII`$ models will be non-chaotic. Evolution could occur which exhibits oscillations about the axisymmetric (LRS) vacuum Bianchi type $`VIII`$ model with $`\mathrm{\Sigma }_{}=0`$. The asymptotic late-time behaviour of this solution is known . The Bianchi type $`VIII`$ LRS model is described by the following restrictions on the state variables, and these define an invariant set for the evolution equations (23)-(26) :

$$N_1<0$$

$$N_2=N_3>0$$

$$\mathrm{\Sigma }_{}=0$$

This state space is unbounded, since $`N_2`$ can take arbitrarily large values. However, we can use the following bounds on $`N_1N_2`$ (see inequalities (28)-(33)):

$$0N_1N_23$$

We introduce a new variable $`Z`$ defined by Wainwright, ,

$$Z^2=N_1N_2$$

It can then be shown that the system of evolution equations (23)-(26) admits the following equilibrium point :

$$N_1=0$$

$$\mathrm{\Sigma }_+=\frac{1}{2}$$

$$Z=\frac{3}{2}$$

This corresponds to a flat plane-wave asymptote, which has $`W=0`$. More precisely, it is the Bianchi type $`III`$ form of flat space-time, with metric

$$ds^2=dt^2+t^2(dx^2+e^{2x}dy^2)+dz^2.$$
(52)
The Bianchi type $`VI_h`$ diagonal plane-wave metric with $`h=1`$ reduces to this line element for a particular choice of parameters.

Since the above equilibrium point solution, (52), is of LRS Bianchi type $`III`$ with $`n_{}^{\alpha }{}_{\alpha }{}^{}=0`$, it is a particular case of the following general metric which includes Bianchi types $`I`$, $`VI_0`$, $`VIII`$, $`V`$ and $`VI_h`$ with $`n_{}^{\alpha }{}_{\alpha }{}^{}=0`$

$`ds^2`$ $`=`$ $`dt^2+(X^2+Y^2w^2)dx^22Y^2we^{(a_0+q_0)x}dxdy`$ (53)
$`+Y^2e^{2(a_0+q_0)x}dy^2+Z^2e^{2(a_0q_0)x}dz^2`$
where $`a_0`$ and $`q_0`$ constants and $`X(t)`$, $`Y(t)`$, $`Z(t)`$, $`w(t)`$ are functions of the synchronous time $`t,`$ and

$$w(t)=2bXY^3Z^1𝑑t$$

where $`b`$ is constant. The line element (52) has,

$$b=0,X=Y=t,Z=1,a_0=q_0=\frac{1}{2}.$$

Now since this line element is a future attractor for the LRS vacuum Bianchi type $`VIII,`$ and since (53) contains non-LRS solutions, it is interesting to perform a stability analysis of this line element for non-LRS perturbation modes. The perturbation technique for special cases of (53) has already been developed by Barrow and Sonoda . In this approach one defines a new time variable as

$$\frac{d}{du}=X(t)\frac{d}{dt}=^{}$$

as well as the new variables $`F=\frac{X^{^{}}}{X}`$, $`G=\frac{Y^{^{}}}{Y}`$ and $`H=\frac{Z^{^{}}}{Z}`$ so that the field equations for (53) in vacuum become

$$F^{^{}}=2(a_{0}^{}{}_{}{}^{2}+q_{0}^{}{}_{}{}^{2})\frac{2b^2X^2}{Y^4Z^2}F(G+H),$$
(54)
$$G^{^{}}=2(a_{0}^{}{}_{}{}^{2}+a_0q_0)+\frac{2b^2X^2}{Y^4Z^2}G(G+H),$$
(55)
$$H^{^{}}=2(a_{0}^{}{}_{}{}^{2}a_0q_0)H(G+H),$$
(56)
with the first integral constraints

$$\frac{2b^2X^2}{Y^4Z^2}=2F(G+H)+2GH(6a_{0}^{}{}_{}{}^{2}+2q_{0}^{}{}_{}{}^{2})$$
(57)
and

$$2Fa_0=q_0(GH)+G+H.$$
(58)
The unperturbed solution is the line element (52) with $`F_0=1`$, $`G_0=1`$ and $`H_0=0`$. We require equations for the following small perturbations, $`f,g,h`$ :

$$F=F_0+f,$$

$$G=G_0+g,$$

$$H=H_0+h.$$

Using the equations (54)-(58), one gets the following (the equation for $`F`$ is eliminated using the constraint (58)) :

$$g^{^{}}=4g4h$$
(59)
$$h^{^{}}=h$$
(60)
The eigenvalues associated with this problem are $`\lambda _1=4`$ and $`\lambda _2=1`$ and so the solution (52) is asymptotically stable under small non-LRS perturbations consistent with (53), thus supporting the prediction that the universe point will evolve about the axisymmetric solution with $`\beta _{}=0`$ in the narrow infinite channel along the $`\beta _+`$-axis. Taking into account the results of the previous section, this means we expect the future asymptotic state of the full non-LRS vacuum Bianchi type $`VIII`$ to have

$$\mathrm{\Sigma }=\frac{1}{2}.$$

The fact that $`\mathrm{\Sigma }_{}`$ vanishes is a consequence of the following : since $`\beta _{}0`$, the line element (34) becomes axisymmetric and so the field equations (5), (6) and (7) lead to the LRS solution for the scale factors, thus $`\beta _+\frac{\tau }{2}`$, however from equation (41) one finds that if $`\mathrm{\Sigma }_{}`$ oscillates and does not decrease, then the curvature parameter $`K`$ can be unbounded, which contradicts the generalised Friedmann equation (20). Also, although $`\mathrm{\Sigma }_{}`$ decreases, it can be shown to undergo increasingly rapid oscillations : if $`\beta _{}`$ is bounded, oscillatory and decreasing, then equations (23), (38) and (41) imply that $`\beta _{}`$ ( and thus $`\mathrm{\Sigma }_{}`$ ) undergoes increasingly rapid oscillations. If the universe point oscillates about the axisymmetric solution, then the metric shares a feature of some gravitational-wave space-times, namely that the solution oscillates about a given form of flat space-time as $`\tau +\mathrm{}`$.

## 6 Conclusion

A qualitative analysis based on the Hamiltonian formalism and the orthonormal frame approach gives the following information about the future asymptotic state of vacuum Bianchi type $`VIII`$ models :

$`a`$. The chaotic Mixmaster solution close to the initial singularity evolves at late times into a simpler non-chaotic solution, contrary to the claims of ref. , and corresponding to the fact that the universe point evolves along the infinite open channel of the Bianchi $`VIII`$ potential and oscillates about the LRS Bianchi type $`VIII`$ solution in such a way that the shear parameter exhibits increasingly rapid oscillations as $`\tau +\mathrm{}`$ .

$`b.`$The form of the line element (34) tends to the Bianchi $`III`$ form of flat space-time because $`\beta _{}0`$. It is interesting to compare this with asymptotic states of other vacuum Bianchi types : Bianchi types $`IV`$, $`VI_h`$, and $`VII_h`$ tend to plane-wave space-times; Bianchi type $`VII_0`$ tends to the special equilibrium point $`N_1=0`$, $`N_2=N_3=`$ constant, with $`\mathrm{\Sigma }_+=1`$ and $`\mathrm{\Sigma }_{}=0`$, which in the state space is a line emanating from a point corresponding to the Taub form of flat space; finally, Bianchi type $`V`$ tends to the Milne form of flat space. So it appears that different forms of flat and plane-wave space-times play an important role in the future state of these eternally expanding vacuum Bianchi models and one might expect that they occur in a leading-order approximation to part of the general inhomogeneous vacuum solution to the Einstein field equations at late times.

## 7 Acknowledgements

We would like to thank John Wainwright for sending us helpful information about the late-time behaviour of LRS Bianchi type $`VIII`$ models.

## 8 Appendix

For a non-vacuum homogenous cosmology the Friedmann equation (20) becomes

$$1=\mathrm{\Sigma }^2+K+\mathrm{\Omega }$$
(61)
where $`\mathrm{\Omega }`$ is the density parameter defined by ($`8\pi G=1`$)

$$\mathrm{\Omega }=\frac{\rho }{3H^2}$$

with $`\rho `$ the energy density of a perfect fluid with equation of state $`p=(\gamma 1)\rho .`$ It appears that if $`0<\gamma \frac{2}{3}`$, the LRS Bianchi $`VIII`$ models admits the flat Friedmann model as future attractor i.e. $`\mathrm{\Sigma }=N_1=N_2=N_3=0`$ and $`\mathrm{\Omega }=1`$, so that the argument developed in section 4 does not hold. This is simply the well known cosmic no hair attractor of power-law inflationary universes whose matter content violates the strong energy condition. In fact, for non-vacuum Bianchi type $`VIII`$, the analysis of section 4 can be adapted as follows when the strong energy condition holds. Note first that the deceleration parameter is now given by

$$q=2\mathrm{\Sigma }^2+\frac{1}{2}(3\gamma 2)\mathrm{\Omega }$$

and so we have $`q0`$ if $`3\gamma 20`$. For all $`\gamma `$ satisfying this inequality, one can apply the analysis of section 4 and the universe point will again evolve in the infinite channel as $`\tau +\mathrm{}`$.

We end this appendix by a remark about the stiff fluid case with $`\gamma =2`$. This case is of interest namely because the stiff fluid is equivalent to the presence of a massless scalar field. The evolution equation for $`\mathrm{\Omega }`$ is

$$\frac{d\mathrm{\Omega }}{d\tau }=[2q(3\gamma 2)]\mathrm{\Omega }$$
(62)
Now when $`\gamma =2`$, the field equations (1),…, (3) are identical for vacuum and non-vacuum cases. So, the inequality (17) is valid and this corresponds to

$$q<2$$

and so equation (62) implies that $`\mathrm{\Omega }0`$ as $`\tau +\mathrm{}`$, so that the evolution tends to a vacuum solution.


text:
# On wave propagation in inhomogeneous systems

## 1 Introduction

The Korringa-Kohn-Rostoker (KKR) method  has been used extensively in the study of the electronic structure and related properties of materials. It has been very successful in calculations of the electronic structure of ordered elemental solids , impurities , and disordered alloys . Based on the multiple-scattering theory , the KKR method calculates from first principles quantities such as the electron density and the ground-state energy of the system. A KKR method for layered systems, the so-called layer KKR (LKKR) method, has also been developed . This method is well adapted to calculations of low-energy electron diffraction (LEED)  and electron emission  spectra.

In recent years, propagation of classical waves \[electromagnetic (EM) or elastic waves\] in composite materials with dielectric or, respectively, elastic properties which are periodic functions of the position, with a period comparable to the wavelength of the corresponding field, has been the object of considerable attention . These materials, photonic and phononic crystals, respectively, whether they exist naturally or are artificially fabricated, exhibit a rich variety of physical properties of interest to fundamental and applied research. There are striking analogies between the propagation of electrons in ordinary crystals and EM/elastic waves in photonic/phononic crystals, so that a great variety of multiple-scattering methods originally developed for electronic-structure calculations have been transferred to the field of photonic and phononic crystals . The aim of the present article is to present the KKR formalism for electron, EM and elastic waves from a unifying point of view.

## 2 Multipole expansion of a wave field

We first consider, as an example of a field with one degree of freedom, the case of a spinless particle of mass $`m`$ and energy $``$, in a potential $`V(𝐫)`$. According to quantum mechanics, the particle is described by a scalar field: the single-particle wave function

$$\mathrm{\Psi }(𝐫;t)=\mathrm{\Psi }(𝐫)\mathrm{exp}\left(it/\mathrm{}\right),$$
(1)
where $`\mathrm{}`$ is the Planck constant and $`\mathrm{\Psi }(𝐫)`$ satisfies the time-independent Schrödinger equation

$$\left[\frac{\mathrm{}^2}{2m}^2+V(𝐫)\right]\mathrm{\Psi }=0.$$
(2)
In a region of constant potential, $`V_0`$, Eq. (2) reduces to the Helmholtz wave equation

$$\left[^2+q^2\right]\mathrm{\Psi }=0,$$
(3)
where $`q=\sqrt{2m(V_0)}/\mathrm{}`$. A complete set of spherical-wave solutions of Eq. (3) is given by

$$\mathrm{\Psi }_\mathrm{}m(𝐫)=f_{\mathrm{}}(qr)Y_{\mathrm{}}^m(\widehat{𝐫}),$$
(4)
where $`Y_{\mathrm{}}^m(\widehat{𝐫})`$ are the usual spherical harmonics and $`f_{\mathrm{}}`$ may be any linear combination of the spherical Bessel function, $`j_{\mathrm{}}`$, and the spherical Hankel function, $`h_{\mathrm{}}^+`$. The most general wave function in a constant potential field can be written as linear sum of the spherical waves given by Eq. (4), as follows:

$$\mathrm{\Psi }(𝐫)=\underset{\mathrm{}m}{}\mathrm{a}_\mathrm{}mf_{\mathrm{}}(qr)Y_{\mathrm{}}^m(\widehat{𝐫}),$$
(5)
where $`\mathrm{a}_\mathrm{}m`$ are coefficients to be determined.

A vector field has in general, not one, but three degrees of freedom. Let us, for example, consider a harmonic elastic wave, of angular frequency $`\omega `$, propagating in a homogeneous medium of mass density $`\rho `$ and Lamé coefficients $`\lambda `$, $`\mu `$. This is described by a displacement vector field:

$$𝐔(𝐫,t)=\mathrm{Re}\left[𝐔(𝐫)\mathrm{exp}(\mathrm{i}\omega t)\right],$$
(6)
where $`𝐔(𝐫)`$ satisfies the equation

$$\left(\lambda +2\mu \right)\left(𝐔\right)\mu \times \left(\times 𝐔\right)+\rho \omega ^2𝐔=0.$$
(7)
The most general solution of Eq. (7) consists of two types of elastic waves which propagate independently. These are purely longitudinal (irrotational) waves ($`\times 𝐔=\mathrm{𝟎}`$), which satisfy the vector Helmholtz equation

$$\left[^2+q_l^2\right]𝐔=0,$$
(8)
where $`q_l=\omega /c_l`$, with $`c_l=\sqrt{(\lambda +2\mu )/\rho }`$ being the longitudinal-wave speed; and purely transverse (divergence-free) waves ($`𝐔=0`$), which satisfy the vector Helmholtz equation

$$\left[^2+q_t^2\right]𝐔=0,$$
(9)
where $`q_t=\omega /c_t`$, with $`c_t=\sqrt{\mu /\rho }`$ being the transverse-wave speed. A complete set of longitudinal ($`L`$) spherical-wave solutions of Eq. (8) is given by

$$𝐔_\mathrm{}m^L(𝐫)=\frac{1}{q_l}\left[f_{\mathrm{}}(q_lr)Y_{\mathrm{}}^m(\widehat{𝐫})\right].$$
(10)
A complete set of transverse ($`M,N`$) spherical-wave solutions of Eq. (9) is given by

$$𝐔_\mathrm{}m^M(𝐫)=f_{\mathrm{}}(q_tr)𝐗_\mathrm{}m(\widehat{𝐫})$$
(11)
and

$$𝐔_\mathrm{}m^N(𝐫)=\frac{\mathrm{i}}{q_t}\times f_{\mathrm{}}(q_tr)𝐗_\mathrm{}m(\widehat{𝐫}),$$
(12)
where $`𝐗_\mathrm{}m(\widehat{𝐫})`$ are the so-called vector spherical harmonics, which are defined by

$$\sqrt{\mathrm{}(\mathrm{}+1)}𝐗_\mathrm{}m(\widehat{𝐫})=𝐋Y_{\mathrm{}}^m(\widehat{𝐫})\mathrm{i}𝐫\times Y_{\mathrm{}}^m(\widehat{𝐫}).$$
(13)
The most general displacement field in a homogeneous medium can be written as a linear sum of the spherical waves given by Eqs. (10)-(12), as follows

$`𝐔(𝐫)=`$ $`{\displaystyle \underset{\mathrm{}m}{}}\{\mathrm{a}_\mathrm{}m^Mf_{\mathrm{}}(q_tr)𝐗_\mathrm{}m(\widehat{𝐫})+\mathrm{a}_\mathrm{}m^N{\displaystyle \frac{\mathrm{i}}{q_t}}\times f_{\mathrm{}}(q_tr)𝐗_\mathrm{}m(\widehat{𝐫})`$ (14)
$`+\mathrm{a}_\mathrm{}m^L{\displaystyle \frac{1}{q_l}}\left[f_{\mathrm{}}(q_lr)Y_{\mathrm{}}^m(\widehat{𝐫})\right]\},`$
where $`\mathrm{a}_\mathrm{}m^P`$, $`P=M,N,L`$, are coefficients to be determined.

Another example of a vector field is the EM field. A harmonic EM wave, of angular frequency $`\omega `$, is described by its electric-field component

$$𝐄(𝐫,t)=\mathrm{Re}\left[𝐄(𝐫)\mathrm{exp}(\mathrm{i}\omega t)\right].$$
(15)
We need not write down explicitly the magnetic-field component of the wave since this can be readily obtained from $`𝐄(𝐫,t)`$ using one of Maxwell equations. In a homogeneous medium characterized by a dielectric function $`ϵ(\omega )ϵ_0`$ and a magnetic permeability $`\mu (\omega )\mu _0`$, where $`ϵ_0`$, $`\mu _0`$ are the electric permittivity and magnetic permeability of vacuum, Maxwell equations imply that $`𝐄(𝐫)`$ satisfies a vector Helmholtz equation, subject to the condition $`𝐄=0`$, with a wave number $`q=\omega /c`$, where $`c=1/\sqrt{\mu ϵ\mu _0ϵ_0}`$ is the velocity of light in the medium. The spherical-wave expansion of $`𝐄(𝐫)`$ is given by the terms of Eq. (14) which satisfy $`𝐄=0`$, i.e.,

$`𝐄(𝐫)={\displaystyle \underset{\mathrm{}m}{}}\left\{\mathrm{a}_\mathrm{}m^Hf_{\mathrm{}}(qr)𝐗_\mathrm{}m(\widehat{𝐫})+\mathrm{a}_\mathrm{}m^E{\displaystyle \frac{\mathrm{i}}{q}}\times \left[f_{\mathrm{}}(qr)𝐗_\mathrm{}m(\widehat{𝐫})\right]\right\},`$ (16)
where $`\mathrm{a}_\mathrm{}m^P`$, $`P=H,E`$, are coefficients to be determined.

## 3 Scattering by a single scatterer

We consider a scatterer of finite range, $`S`$, with its center at the origin of coordinates, and assume that the appropriate in each case characteristics (potential, Lamé coefficients and mass density, electric permittivity and magnetic permeability) are different from those of the surrounding homogeneous medium. A De Broglie, EM or an elastic plane wave incident on this scatterer is described, respectively, by Eq. (4), (16) or (14) with $`f_{\mathrm{}}=j_{\mathrm{}}`$ (since the plane wave is finite everywhere) and appropriate coefficients $`\mathrm{a}_L^0`$, where $`L`$ denotes collectively the indices $`\mathrm{}mP`$. Similarly, the scattered wave is described by Eq. (4), (16) or (14) with $`f_{\mathrm{}}=h_{\mathrm{}}^+`$, which has the asymptotic form appropriate to an outgoing spherical wave: $`h_{\mathrm{}}^+(\mathrm{i})^{\mathrm{}}\mathrm{exp}(\mathrm{i}qr)/\mathrm{i}qr`$ as $`r\mathrm{}`$, and appropriate expansion coefficients $`\mathrm{a}_L^+`$. The wavefield for $`r>S`$ is the sum of the incident and scattered waves. The spherical-wave expansion of the field for $`r<S`$ is obtained in similar manner by the requirement that it be finite at the origin ($`𝐫=\mathrm{𝟎}`$). By applying the proper in each case boundary conditions at the surface of the scatterer, we obtain a relation between the expansion coefficients of the incident and the scattered field, as follows:

$$\mathrm{a}_L^+=\underset{L^{}}{}T_{LL^{}}\mathrm{a}_L^{}^0,$$
(17)
where $`T_{LL^{}}`$ are the elements of the so-called scattering transition matrix. Eq. (17) is valid for any shape of scatterer; for spherically symmetric scatterers each spherical wave scatters independently of all others, which leads to a transition matrix diagonal in angular momentum. Explicit forms for $`T_{LL^{}}`$ for various cases can be found elsewhere .

## 4 Multiple scattering theory

We consider an assembly of non-overlapping scatterers centered on sites $`𝐑_i`$ in a homogeneous medium. In general: an outgoing wave about $`𝐑_i^{}`$ \[described by Eq. (4), (14) or (16) with $`f_{\mathrm{}}=h_{\mathrm{}}^+`$ and expansion coefficients $`\mathrm{b}_L^{+i^{}}`$\] can be expanded into spherical waves about $`𝐑_i`$, incident on $`𝐑_i`$. This expansion has the form of Eq. (4), (14) or (16) with $`f_{\mathrm{}}=j_{\mathrm{}}`$ and expansion coefficients $`\mathrm{b}_{}^{}{}_{L}{}^{i}(i^{})`$, given by

$$\mathrm{b}_{}^{}{}_{L}{}^{i}(i^{})=\underset{L^{}}{}\mathrm{\Omega }_{LL^{}}^{ii^{}}\mathrm{b}_L^{}^{+i^{}}.$$
(18)
Explicit expressions for the so-called free-space structural Green functions, $`\mathrm{\Omega }_{LL^{}}^{ii^{}}`$, for the scalar and transverse vector fields can be found elsewhere . In the case of the elastic field, since a longitudinal (transverse) spherical wave about $`𝐑_i^{}`$ remains a longitudinal (transverse) wave when expanded about another center $`𝐑_i𝐑_i^{}`$, the matrix elements $`\mathrm{\Omega }_{LL^{}}^{ii^{}}`$ are obtained independently for longitudinal and transverse waves. Therefore, $`\mathrm{\Omega }_{LL^{}}^{ii^{}}`$ has a block-diagonal form in the mode index $`P`$ (we recall that $`L\mathrm{}mP`$): the matrix elements of the block corresponding to transverse waves are identical to the corresponding matrix elements for the EM field, while the matrix elements of the block corresponding to longitudinal waves are identical to the corresponding matrix elements for the scalar field . We note that, for any field, by definition, $`\mathrm{\Omega }_{LL^{}}^{ii^{}}`$ equals zero for $`i=i^{}`$.

The wave scattered from the scatterer at $`𝐑_i`$ is determined from the total wave incident on this scatterer; therefore

$$\mathrm{b}_L^{+i}=\underset{L^{}}{}T_{LL^{}}^i(\mathrm{a}_L^{}^{0i}+\mathrm{b}_{}^{}{}_{L^{}}{}^{i}),$$
(19)
where $`\mathrm{a}_L^{0i}`$ are the coefficients in the multipole expansion \[given by Eq. (4), (14) or (16)\] about $`𝐑_i`$ of an externally incident wave; and $`\mathrm{b}_{}^{}{}_{L}{}^{i}`$ are the coefficients in the multipole expansion about $`𝐑_i`$ of the outgoing waves scattered from all the other scatterers at $`𝐑_i^{}𝐑_i`$. From Eq. (18) we have

$$\mathrm{b}_{}^{}{}_{L}{}^{i}=\underset{i^{},L^{}}{}\mathrm{\Omega }_{LL^{}}^{ii^{}}\mathrm{b}_L^{}^{+i^{}}.$$
(20)
Substituting Eq. (20) into Eq. (19) we obtain

$$\underset{i^{},L^{}}{}\left(\delta _{ii^{}}\delta _{LL^{}}\underset{L^{\prime \prime }}{}T_{LL^{\prime \prime }}^i\mathrm{\Omega }_{L^{\prime \prime }L^{}}^{ii^{}}\right)\mathrm{b}_L^{}^{+i^{}}=\underset{L^{}}{}T_{LL^{}}^i\mathrm{a}_L^{}^{0i}.$$
(21)
We now introduce the structural Green functions $`D_{LL^{}}^{ii^{}}`$ which give the coefficients \[in an expansion such as (4), (14) or (16)\] of the wave incident on the scatterer at $`𝐑_i`$, due to an outgoing wave from the scatterer at $`𝐑_i^{}`$. An outgoing wave from the $`i^{}`$th scatterer can reach the $`i`$th scatterer directly, or indirectly after scattering any number of times by any number of scatterers (including those at $`𝐑_i`$ and $`𝐑_i^{}`$). Let $`D_{LL^{}}^{ii^{}}`$ express the sum of the contribution to the coefficients of the incident wave on $`𝐑_i`$ from all possible scattering paths originating from an outgoing wave from the $`i^{}`$th scatterer. One can easily prove by iteration the following equation

$$D_{LL^{}}^{ii^{}}=\mathrm{\Omega }_{LL^{}}^{ii^{}}+\underset{i^{\prime \prime },L^{\prime \prime },L^{\prime \prime \prime }}{}\mathrm{\Omega }_{LL^{\prime \prime }}^{ii^{\prime \prime }}T_{L^{\prime \prime }L^{\prime \prime \prime }}^{i^{\prime \prime }}D_{L^{\prime \prime \prime }L^{}}^{i^{\prime \prime }i^{}},$$
(22)
where $`T_{L^{\prime \prime }L^{\prime \prime \prime }}^{i^{\prime \prime }}`$ are the elements of the transition matrix which describes the scattering by the single scatterer at $`𝐑_{i^{\prime \prime }}`$. We generalize Eq. (22) by treating the scattering at $`𝐑_i`$ in two stages. The first-stage scattering is described by $`T_{LL^{}}^{ri}`$, which correspond to arbitrarily defined scatterers (reference scatterers), and the second-stage scattering by $`T_{LL^{}}^iT_{LL^{}}^{ri}`$. We obtain:

$$D_{LL^{}}^{ii^{}}=D_{LL^{}}^{rii^{}}+\underset{i^{\prime \prime },L^{\prime \prime },L^{\prime \prime \prime }}{}D_{LL^{\prime \prime }}^{rii^{\prime \prime }}\left(T_{L^{\prime \prime }L^{\prime \prime \prime }}^{i^{\prime \prime }}T_{L^{\prime \prime }L^{\prime \prime \prime }}^{ri^{\prime \prime }}\right)D_{L^{\prime \prime \prime }L^{}}^{i^{\prime \prime }i^{}},$$
(23)
where $`D_{LL^{}}^{rii^{}}`$ are the solution of Eq. (22) when $`T_{LL^{}}^i=T_{LL^{}}^{ri}`$.

## 5 Bulk systems

We consider the case in which the scatterers constitute a three dimensional (3D) crystal structure specified by Bravais lattice vectors $`𝐑_n`$ and non-primitive translation vectors, indicating the positions of the scatterers within the unit cell, $`𝐭_\alpha `$; in this case the site index $`i`$ stands for the composite index $`n\alpha `$. Bloch’s theorem implies that $`\mathrm{b}_L^{}^{+n^{}\alpha ^{}}=\mathrm{exp}\left[\mathrm{i}𝐤(𝐑_n^{}𝐑_n)\right]\mathrm{b}_L^{}^{+n\alpha ^{}}`$.

The normal modes of the crystal are obtained, by putting the external incident wave equal to zero in Eq. (21), which leads to the following secular equation:

$$det\left[\delta _{\alpha \alpha ^{}}\delta _{LL^{}}\underset{L^{\prime \prime }}{}T_{LL^{\prime \prime }}^\alpha \mathrm{\Omega }_{L^{\prime \prime }L^{}}^{\alpha \alpha ^{}}(𝐤)\right]=0,$$
(24)
where

$$\mathrm{\Omega }_{LL^{}}^{\alpha \alpha ^{}}(𝐤)=\underset{n^{}}{}\mathrm{\Omega }_{LL^{}}^{n\alpha ;n^{}\alpha ^{}}\mathrm{exp}\left[\mathrm{i}𝐤(𝐑_n^{}𝐑_n)\right],$$
(25)
which does not depend on $`n`$. Both $`T_{LL^{}}^\alpha `$ and $`\mathrm{\Omega }_{LL^{}}^{\alpha \alpha ^{}}(𝐤)`$ in Eq. (24) are functions of the frequency of the wave, but the $`T_{LL^{}}^\alpha `$ depend only on the properties of a single scatterer, whereas $`\mathrm{\Omega }_{LL^{}}^{\alpha \alpha ^{}}(𝐤)`$ depend only on the geometry. Thus Eq. (24) reflects a complete separation of the individual scatterer and structural aspects of the problem. The matrix elements $`\mathrm{\Omega }_{LL^{}}^{\alpha \alpha ^{}}(𝐤)`$ for the scalar field were introduced by Korringa, Kohn, and Rostoker  in relation to the calculation of the electronic band structure of periodic solids and are commonly referred to as the KKR structure constants. The calculation of the structure constants, which needs to be done only once for a given lattice, usually requires Ewald-summation techniques . We note that Eq. (24) involves in principle infinite-dimensional matrices. In actual calculations, however, it is sufficient to truncate the angular momentum index, $`\mathrm{}`$, to some relatively small number, $`\mathrm{}_{max}`$. The KKR method has been successfully used for many years in the calculation of the electron band structure of solids , and more recently it has also been applied to the calculation of the frequency band structure of EM and elastic fields in relation to photonic and phononic crystals .

One of the advantages of the KKR theory compared to other approaches is its capability of dealing with defects and disorder . For the description of point defects at a finite number of sites, one needs to calculate the structural Green functions, $`D_{LL^{}}^{ii^{}}`$, of the defect system. This can be done in real space, using Eq. (23) and considering the periodic crystal (without any defect) as the reference system. In this case, the sum over $`i^{\prime \prime }`$ in Eq. (23) is restricted to those sites at which there are defects and induce a perturbation $`\delta T_{L^{\prime \prime }L^{\prime \prime \prime }}^{i^{\prime \prime }}=T_{L^{\prime \prime }L^{\prime \prime \prime }}^{i^{\prime \prime }}T_{L^{\prime \prime }L^{\prime \prime \prime }}^{ri^{\prime \prime }}`$ on the $`T`$ matrix. For the periodic arrangement of scatterers of the reference system, the evaluation of $`D_{LL^{}}^{rii^{}}`$ through Eq. (22) involves an infinite number of sites $`i^{\prime \prime }`$ and it is achieved by a lattice Fourier transform of Eq. (22)

$$D_{LL^{}}^{r\alpha \alpha ^{}}(𝐤)=\mathrm{\Omega }_{LL^{}}^{\alpha \alpha ^{}}(𝐤)+\underset{\alpha ^{\prime \prime },L^{\prime \prime },L^{\prime \prime \prime }}{}\mathrm{\Omega }_{LL^{\prime \prime }}^{\alpha \alpha ^{\prime \prime }}(𝐤)T_{L^{\prime \prime }L^{\prime \prime \prime }}^{r\alpha ^{\prime \prime }}D_{L^{\prime \prime \prime }L^{}}^{r\alpha ^{\prime \prime }\alpha ^{}}(𝐤)$$
(26)
and subsequent integration over the volume $`v`$ of the Brillouin zone (BZ)

$$D_{LL^{}}^{rii^{}}=\frac{1}{v}_{\mathrm{BZ}}\mathrm{d}^3k\mathrm{exp}\left[\mathrm{i}𝐤\left(𝐑_n𝐑_n^{}\right)\right]D_{LL^{}}^{r\alpha \alpha ^{}}(𝐤).$$
(27)
If a number of $`N`$ different scatterers, described by $`T`$ matrices $`T_{LL^{}}^j`$, $`j=1,2,\mathrm{},N`$, are randomly distributed on the sites of a given crystal, the coherent potential approximation (CPA)  can be used to treat the disordered system within a mean-field context. The effective CPA medium consists of identical scatterers at all sites of sublattice $`\alpha `$ of the crystal, characterized by a $`T`$ matrix $`T_{LL^{}}^{c\alpha }`$. The structural Green functions $`D_{LL^{}}^{cn\alpha ;n^{}\alpha ^{}}`$ of the CPA medium are calculated by equations analogous to (26) and (27). The CPA self-consistency condition demands that the correction to the scattering due to the difference of the actual scatterers from the CPA scatterer vanishes on the average

$$\underset{j=1}{\overset{N}{}}c_{j;\alpha }\left[𝐈\left(𝐓^j𝐓^{c\alpha }\right)𝐃^{c\mathrm{\hspace{0.17em}\hspace{0.17em}0}\alpha ;0\alpha }\right]^1\left(𝐓^j𝐓^{c\alpha }\right)=0\mathrm{for}\mathrm{every}\alpha ,$$
(28)
where $`c_{j;\alpha }`$ denotes the concentration of scatterer $`j`$ on sublattice $`\alpha `$.

## 6 Layered systems

A layered system consists of a number of planes (layers) of scatterers with the same two dimensional (2D) periodicity. To begin with, we consider just one layer, at $`z=0`$, in which case the scatterers are centered on the sites $`𝐑_n+𝐭_\alpha `$, where $`\{𝐑_n\}`$ constitutes a 2D Bravais lattice and $`𝐭_\alpha `$ denote the positions of the scatterers within the 2D unit cell. The 2D reciprocal vectors $`𝐠`$, and the surface Brillouin zone (SBZ) corresponding to this lattice are defined in the usual manner . The wave vector parallel to the plane of atoms can always be written as follows:

$$𝐪_{}=𝐤_{}+𝐠^{},$$
(29)
where the reduced wave vector, $`𝐤_{}`$, lies in the SBZ and $`𝐠^{}`$ is a certain reciprocal vector. In what follows we write the wave vector of a plane wave of given $`𝐪_{}=𝐤_{}+𝐠`$ and given wave number $`q_\nu `$, where $`\nu `$ specifies the polarization mode of the wave, if there exist longitudinal and transverse modes, as follows:

$$𝐊_{𝐠\nu }^\pm =(𝐤_{}+𝐠,\pm \left[q_\nu ^2(𝐤_{}+𝐠)^2\right]^{1/2}).$$
(30)
We note that when $`q_\nu ^2<(𝐤_{}+𝐠)^2`$, the above defines a decaying wave; the positive (negative) sign in Eq. (30) corresponds to a wave propagating or decaying to the right (left).

The coefficients $`\mathrm{b}_L^{+n\alpha }`$ in the multipole expansion of the scattered field from the layer are obtained from Eq. (21). Using the Bloch condition: $`\mathrm{b}_L^{}^{+n^{}\alpha ^{}}=\mathrm{exp}\left[\mathrm{i}𝐤_{}(𝐑_n^{}𝐑_n)\right]\mathrm{b}_L^{}^{+n\alpha ^{}}`$, we obtain

$$\underset{\alpha ^{},L^{}}{}\left[\delta _{\alpha \alpha ^{}}\delta _{LL^{}}\underset{L^{\prime \prime }}{}T_{LL^{\prime \prime }}^\alpha \mathrm{\Omega }_{L^{\prime \prime }L^{}}^{\alpha \alpha ^{}}(𝐤_{})\right]\mathrm{b}_L^{}^{+n\alpha ^{}}=\underset{L^{}}{}T_{LL^{}}^\alpha \mathrm{a}_L^{}^{0n\alpha },$$
(31)
where $`\mathrm{a}_L^{0n\alpha }`$ are the coefficients in the multipole expansion of the incident plane wave. The 2D KKR structure constants

$$\mathrm{\Omega }_{LL^{}}^{\alpha \alpha ^{}}(𝐤_{})=\underset{n^{}}{}\mathrm{\Omega }_{LL^{}}^{n\alpha ;n^{}\alpha ^{}}\mathrm{exp}\left[\mathrm{i}𝐤_{}(𝐑_n^{}𝐑_n)\right],$$
(32)
like their counterparts in 3D given by Eq. (25), can be evaluated using Ewald-summation techniques .

Writing the incident, reflected and transmitted plane waves with respect to an origin (a scattering center) in this plane, we obtain the amplitudes, $`\left[\mathrm{\Phi }_{\mathrm{rf}}\right]_{𝐠i}^\pm `$ and $`\left[\mathrm{\Phi }_{\mathrm{tr}}\right]_{𝐠i}^\pm `$ of the reflected and transmitted beams, respectively, in terms of the amplitudes, $`\left[\mathrm{\Phi }_{\mathrm{in}}\right]_{𝐠^{}i}^\pm `$, of the incident wave as follows:

$`\left[\mathrm{\Phi }_{\mathrm{rf}}\right]_{𝐠i}^s`$ $`={\displaystyle \underset{i^{}}{}}M_{𝐠i;𝐠^{}i^{}}^{ss}\left[\mathrm{\Phi }_{\mathrm{in}}\right]_{𝐠^{}i^{}}^s`$
$`\left[\mathrm{\Phi }_{\mathrm{tr}}\right]_{𝐠i}^s`$ $`={\displaystyle \underset{i^{}}{}}M_{𝐠i;𝐠^{}i^{}}^{ss}\left[\mathrm{\Phi }_{\mathrm{in}}\right]_{𝐠^{}i^{}}^s,`$ (33)
where $`s=+()`$ specifies a wave traveling or decaying to the positive (negative) $`z`$ direction, and $`i`$ denotes the components of the wave field. A scalar field ($`𝚽\mathrm{\Psi }`$) has just one component ($`i=1`$), the displacement vector of the elastic field ($`𝚽𝐔`$) has three components ($`i=1,2,3`$), while the EM field ($`𝚽𝐄`$), because of its transverse nature, has only two independent components, $`i=1,2`$. Explicit expressions for the transmission/reflection matrix elements $`M_{𝐠i;𝐠^{}i^{}}^{ss^{}}`$ for the scalar, elastic, and EM fields can be found elsewhere .

The transmission/reflection matrices for a slab which consists of a stack of layers with the same 2D periodicity parallel to a given plane are obtained from the transmission/reflection matrices of the individual layers in the manner described in Refs. . Knowing the transmission/reflection matrices for the slab we can readily obtain the transmission, reflection, and absorption coefficients of the plane wave of Eq. (30) incident on the slab. The LKKR technique summarized in this section is particularly well adapted to layered systems because the computation time scales linearly to the number of layers in the slab. Moreover, because this method does not require periodicity in the direction perpendicular to the layers, it can easily treat heterostructures, slabs with impurity planes  or even a random succession of different layers, as long as these have the same 2D periodicity.

The LKKR technique provides also the complex band structure of an infinite crystal, viewed as a succession of layers parallel to a given crystallographic plane .

Point defects or a random distribution of point defects in the layers of a slab can also be treated along the lines described in Section 5 .

I. E. Psarobas was partly supported by the University of Athens, and partly by the Institute of Communication and Computer Systems (ICCS) of the National Technical University of Athens. V. Yannopapas was supported by the State Scholarship Foundation of Greece (I.K.Y.).