Quantum Field Theory of Exciton Correlations and Entanglement in Semiconductor Structures

The usual ideology when dealing with many-body systems using second quantization is to treat elementary excitations, say electrons, as the main entities determining the physical properties of the system: the states are enumerated by population numbers, the response is described using Green’s functions defined in terms of the elementary excitation, and so on. The situation changes drastically when interaction allows for the formation of bound states. The whole view of the system has to be revisited. Perhaps the most famous example is given by superconductors. In the normal state the material fits the canonical universal description when, roughly speaking, most of the properties are explained (at least qualitatively) by the position of the Fermi level. In the superconducting state, however, one has complete reconstruction of the ground state. Now properties of thematerial are defined by Cooper pairs with behavior qualitatively different from that of individual electrons. For one the exclusion principle has significantly diminished effect, so that pairs can even condense. Such change in the character of elementary excitations leads to significant consequences: resistivity drops practically to zero.

approach. Also, we discuss some nonlinear effects in the excitonic systems obtained by means of the TDDFT approach.

Semiconductor optical response 2.1 Dynamics of semiconductor excitations
In the most general setup the problem of combined time evolution of semiconductor excitations and electromagnetic field requires considering very complex Hamiltonians of electrons in a periodic lattice and coupled to quantized field, whose dynamics, in turn, may be affected by the spatial variation of the dielectric function. However, the physics of phenomena of our main interest, as will be seen, is quite rich on its own. Therefore, in order to demonstrate the major effects we consider the simplest case of a semiconductor with well separated spherical bands excited by a classical electromagnetic field. An example of such a material is GaAs, where the approximation of spherical bands proved to be good. Keeping in mind the applications of techniques discussed below for GaAs we consider the case when the states in the valence band are characterized by the projections σ = ±3/2 of angular momentum (heavy holes).
We will be interested in interband optical transitions, when the external electromagnetic field has its frequency tuned close to the value of the gap separating the conduction and valence bands. Due to the high main frequency the vector potential can be presented in a form convenient for adopting the rotating wave approximation A(x, t)=A Ω (x, t)e −iΩt + A * Ω (x, t)e iΩt with relatively slowly changing amplitude A Ω (x, t). The dynamics of excitations is described by the semiconductor-light Hamiltonian H = H SC + H F + H exc , which is composed of the Hamiltonians of the nonperturbed semiconductor H SC , the free electromagnetic field H E and the light-matter interaction H exc . We will concentrate on interband optical transitions so that the frequency of the relevant photon modes is close to the value of the gap (throughout the chapter we use units withh = 1). If the external field is not too strong such an interaction is well described by the rotating wave approximation that takes into account only resonant transitions. Adopting this approximation and using the coordinate representation we have in the rotating frame Here we have introduced c † s (x) and v σ (x), operators creating electron with spin s and hole in spin state σ, respectively, at point x. These are fermion operators satisfying the canonical anticommutation relations {c s 1 (x 1 ), c s 2 ( The first line in Eq. (1) represents the Hamiltonian of non-interacting electrons and holes with m e,h being the respective masses and ǫ c and ǫ v denoting the positions of the bottom of the conduction band and of the top of the valence band, respectively, so that ǫ c − ǫ v is the gap. The last three lines describe the electrostatic interaction with the Coulomb potential V(x).
The Hamiltonian of interaction with external electromagnetic field is where d s,σ = d * σ,s = −i s| ∇ |σ e/m 0 ,w i t hm 0 the electron mass in vacuum, quantifies coupling between the respective states in conduction and valence bands. Interaction of light with the semiconductor occurs through absorption and emission of electron-hole pairs thus specifying the quantity of main interest.
Because the main technical tool used in this part is perturbation theory, it is convenient to incorporate the time dependence into the Heisenberg representation of quantum operators O = e iH SC t Oe −iH SC t ,sothati∂O/∂t =[H SC , O]. For electron-hole pair operator the equation of motion has the form where operator L s,σ (x 1 , x 2 ) describes the one-pair dynamics and U (x) accounts for interaction of the pair with surrounding charges In the case when its contribution vanishes (for instance, when operators in Eq. (4) act on vacuum state, i.e. state with empty conduction band and filled valence band) the semiconductor dynamics becomes very simple and the structure of excitations is determined by the spectral decomposition where the formal summation over μ implies summing over the discrete quantum numbers and integrating over continuous ones. In Eq. (7) E μ and φ μ are, respectively, the eigenvalues and the eigenfunctions, Lφ μ = E μ φ μ . The operator L is self-adjoint and, hence, its eigenvalues are real and the eigenfunctions form a complete orthonormal set. Convoluting both sides of Eq. (4) (taken with U = 0) with φ μ we find that operators Such electron-hole states are called excitons and, respectively, B † μ are called exciton operators. As follows from Eq. (7) excitons are characterized by a set of quantum numbers describing the solutions of the respective Schrödinger equation. For convenience we also include the spin variables into this set, so that μ = {s μ , σ μ , n μ , ℓ μ , K μ },w h e r eK μ is the momentum of the exciton center of mass, n μ is the principal quantum number, and ℓ μ is orbital momentum.
Among the full variety of exciton states not all of them play an equally important role in the dynamics of a semiconductor interacting with external electromagnetic field. In order to see this we employ the fact that wave-functions φ μ form a complete set and obtain the relation between electron-hole and exciton operators where the sum is taken for fixed values of electron and hole spins. This relation allows us to express H exc in terms of the exciton operators where A μ (t) are the projections of the external field onto the respective exciton mode As a result only states with ℓ = 0 are directly coupled with the electromagnetic field. Among four possible heavy-hole exciton states two of them are dark because the respective transitions between valence and conduction bands are dipole-forbidden (Ivchenko, 2005). This results in a life-time of dark excitons that is significantly greater than the life time of bright states. Additionally, the interaction of bright excitons with light raises their energy compared to dark excitons (Combescot & Leuenberger, 2009). These circumstances are very important from the perspective of exciton Bose-Einstein condensation as will be discussed below.

Dark excitons
A simple way to see why the dark exciton states have lower energies than the bright exciton states (Combescot & Leuenberger, 2009) is to rearrange the electron-hole exchange scattering diagram as shown in Fig. 1. Then it becomes obvious that the electron-hole exchange corresponds to the exchange of a virtual photon between electron-hole pairs. Since only bright excitons interact with photons, it is only possible for bright excitons to exchange virtual photons among each other. This interaction pushes the energies of the bright excitons above the energies of the dark excitons. The reason for adding energy can be intuitively understood by comparing direct and exchange Coulomb interactions for electrons; i.e. the direct interaction among electrons leads to repulsion, corresponding to adding energy, while the exchange interaction leads to attraction, corresponding to reducing energy. In the case of the Coulomb interaction between electrons and holes, this situation is completely reversed; i.e. the direct interaction between electrons and holes leads to attraction, corresponding to reducing energy, while the exchange interaction leads to repulsion, corresponding to adding energy. It is this exchange-based repulsion which lets the bright exciton energies be higher than the dark exciton energies. Interband Coulomb exchange interaction which shifts the bright exciton energy above the dark exciton energy. Transition process either in terms of valence-conduction electrons (b) and (c), or in terms of electron-hole (d). This interband Coulomb process is nothing but an exchange of virtual photon between electron-hole pairs.

Perturbation theory
Coupled operator equations of motion of semiconductor excitations and electromagnetic field, turn out to be too complex for a detailed analysis and, therefore, a reliable approximation scheme must be applied. While such scheme can be worked out at the level of the equations of motion it is convenient to start from the Hamiltonian formulation of the problem and to use the standard approach developed in the quantum field theory on the ground of interaction representation.
The equation of motion of the external electromagnetic field driven by semiconductor excitations has the form of classical wave equation with a source This shows that the effect of material excitations is described by the exciton polarization where ρ(t) is the density matrix at instant t. Specifying the semiconductor state by the density matrix allows one to address the most general situation. We, however, start from the special case when the system at all times is in some pure state |Ψ(t) ,sothatP μ (t)= Ψ(t)| B μ |Ψ(t) . Next, we treat the interaction with the electromagnetic field, H exc as a perturbation and follow the standard prescription. First, we account the nonperturbed dynamics by introducing Iterating the equation of motion we find where In what follows we will need the explicit form of such expansion. For a compact notation, however, it is convenient to introduce the time ordering operator T + , so that one has S (t)=T + exp −i t 0 dt ′ H exc (t ′ ) . Following the same line of arguments we can also derive Ψ(t) = Ψ(0) S † (t),w h e r eS † (t)=T − exp i t 0 dt ′ H exc (t ′ ) . Thus we obtain for the exciton polarization (and, actually, for any observable) Finally, the perturbational series for P μ (t) is obtained substituting instead of S (t) and S † (t) their expansions following Eq. (17).
This was the general consideration, which is applicable for arbitrary system and perturbation. For the problem of our main interest, however, the fact of great importance is that exciton operators entering H exc change the total number of electrons and holes. In particular, if the initial state |Ψ(0) is characterized by a definite number of particles this form of perturbation implies that only terms with matching numbers of exciton creation and annihilation operators would make nonzero contribution. At the same time a lot of terms enter even low orders of the perturbation series [as illustrated by Eq. (21) below]. Therefore, for analysis of the series it is convenient to represent terms graphically using as building blocks The raising (lowering) line corresponds to the exciton annihilation (creation) operator in the interaction representation at particular instant and a hollow vertex attached to the line corresponds to the external field taken at the same instant. Thus the order of the diagram is determined by the number of vertices. Integration over time is shown by filled vertex and taking the diagonal matrix element is indicated by a horizontal line, for example While writing down an expression corresponding to particular term in the perturbational series for P μ (t) one should take into account that elements to the left and to the right from B μ originate from expansions for S † and S, respectively, which have factors i and −i in correspondence rule (19). The diagrammatical representation of the first few terms of the perturbational series for P μ (t) has the form

Linear and nonlinear responses of initially unperturbed semiconductor
We illustrate the application of this analysis by considering the optical response of a non-perturbed semiconductor (Erementchouk & Leuenberger, 2010b;Ostreich et al., 1998), that is when the initial state is vacuum |0 (empty conduction band and filled valence band). We will emphasize this fact by using the dashed horizontal (vacuum) line in diagrams. It is seen that only diagrams starting and ending at the vacuum line and above it provide non-zero contributions into the series. The total number of lines in this case is even and, hence, the total number of vertices is odd implying that we have only odd orders of the perturbation theory in this case. In particular among the diagrams shown in Eq. (21) only the second diagram survives yielding the polarization of linear response For the following consideration it is constructive to analyze this expression in somewhat excessive details. We introduce the vacuum exciton propagator Φ and differentiating with respect to τ we find that it satisfies the dynamical equatioṅ with the initial value Φ μ,ν (0)=δ μ,ν , where we have taken into account that B μ B † ν |0 = δ μ,ν |0 following from the commutation relation for the exciton operators Operators C μ,ν describe the deviation of excitons from bosons Thus we find Simple but important feature of the linear response is that wave vectors containing in the external excitation directly transferred to the linear exciton polarization. For example, if In the next (third) order the exciton polarization is given by the following diagrams In order to present this expression in more familiar form we differentiate Eq. (27) with respect to t remembering that we also need to differentiate the upper limits of the respective integrals. This yieldṡ where the elements with the hollow vertices are taken at the instant t and the respective diagrams describe the modification of the instantaneous effect of the electromagnetic field and thus account for the phase-space filling effect. It can be seen that fifth and fourth diagram cancel each other by virtue of [B μ , B ν ]=0. The first and the second diagrams This term is canceled by the commutator appearing after combining the third and the sixth diagrams. Thus the phase-space filling effect is described by The effect of the Coulomb interaction is described by the last term in Eq. (28), which has the form where we have introduced D λ,μ = B λ , B μ , H SC . This operator can be presented as . Thus operator D λ,μ describes the Coulomb interaction between excitons λ and μ. Taking into account the symmetry with respect to ν ↔ κ we can rewrite M μ (t) in terms of polarizations of linear response resulting in where . These coefficients contain typical average of the form B λ B μ B † κ B † ν , which determines the spin selection rules governing different optical processes. In order to derive them we present this average in terms of the electron and hole operators and rearrange operators to have vvv † v † 0 ccc † c † 0 and then expand each term using the Wick theorem for fermions. It produces four terms, which are represented by the diagrams in Fig. 2. The points on the upper line of a diagram represent the spin states of electrons in the conduction band and the points on the lower line stand for the spin states of holes in the valence band. For example, the anticommutator {v † σ κ , v σ μ } ∝ δ σ κ ,σ μ requires the equality of the respective hole spins in the valence band. We denote this equality by connecting the vertices κ and μ on the lower line by the arc.
Upper and lower lines correspond to electron and hole spins, respectively. Arcs connecting two vertices denote equal spins. The diagrams (a) and (b) lead to helicity selection rules ∝ δ κ,σ ν δ σ λ ,σ μ and ∝ δ σ κ ,σ μ δ σ λ ,σ ν , respectively. Diagrams (c) and (d) enter the average with minus sign and for bright excitons require the spin states in the conduction band to be the same.
Also, diagrams in Fig. 2 show which coordinates are identified by delta-functions appearing after anti-commutation and thus demonstrate how exciton wave functions are convoluted in such averages. Thus diagrams in Fig. 2a and 2b describe direct exciton scattering, while those shown in Fig. 2c and 2d take into account scattering with exchange by electron or hole.
The memory term in Eq. (32) accounts for the effect of exciton-exciton interaction. Unfortunately, an exact evolution of this integral is impossible (it is related to four-particle propagator) and one has to rely on approximation schemes. It should be noted, however, that the most significant effect of this interaction is when excitons form a bound state (possibly metastable). When the contribution of such a state is small or non-existent (e.g. there are no bound states in the spectrum of two copolarized excitons) one can employ a short memory approximation, which accounts for the effect of the nonlocal term as a modification of β ν,κ λ,μ .
The essential difference between linear and nonlinear responses is that the latter is a combination of several linear polarizations. As a result if the external excitation has components with different wave vectors the nonlinear polarization contains not only all of them but also their combinations. More precisely one can easily show that P . The four momenta have to add up to zero and therefore this is called four-wave mixing response. Because the direction of the response is different from the direction of excitation this allows observing the effect of interaction that is not blurred by linear response or non-absorbed light. This makes four-wave mixing spectroscopy possible.
We would like to emphasize the generality of the derivation of the exciton optical response. For example, the specific form of the exciton states has not been used for the derivation of main formulas and therefore the same arguments can be repeated for excitons in arbitrary confinement potential. This allows using this approach for description of nonlinear optical response of disordered quantum wells (Erementchouk et al., 2011).
Another important feature of the diagrammatic representation is that it establishes the connection between nonlinear optical response and other phenomena which involve exciton dynamics. As an example we would like to discuss the problem of entanglement of photons interacting with a semiconductor quantum well. For simplicity we neglect the possible effect of variation of the dielectric function and assume that the states of the quantized electromagnetic field are plane waves so that the vector potential is presented as where k = {σ, k} combines all photon quantum numbers, polarization and wave vector. Then the two-photon states are described by the density matrix, which in the interaction picture has the form where |Ψ(t) is the state of the semiconductor-photon system. As the first approximation it suffices to consider vacuum as the initial state of the semiconductor |Ψ(0) = |0 and to neglect the processes of photon re-absorption, which is justified if the photon lifetime within the quantum well is short (it should be noted that the situation may change in a cavity). In the lowest approximation the photon annihilation operators in Eq. (34) act on a two-photon state where a k (t) are the photon operators in the Heisenberg representation. The interaction Hamitonian in this case has the form H Thus the same diagrams as before can be drawn with the only difference that B line is accompanied with A † . Quick analysis shows that only two diagrams contribute into Ψ It is immediately seen that the first diagrams describes the emission of photons along the direction of external excitation and the emitted photons are disentangled: the polarization two-photon state is the direct product of the external excitation polarization. The two-exciton process represented by the second diagram, however, admits oblique emission with non-trivial dependence of entanglement of emitted photons on the direction of observation and the polarization state of the excitation field (Erementchouk & Leuenberger, 2010a) as is summarized in Fig. 3. The polarization of external excitation is described using Poincare sphere, that is the excitation is presented as a combination of left-and right-polarized components with amplitudes A − = e −iχ/2 sin(β/2) and A + = e iχ/2 cos(β/2), respectively, where β is the polar angle on the Poincare sphere and χ is azimuthal angle (χ/2 is the angle between the axis of the ellipse of polarization and the plane spanned by the wave vectors of emitted photons). Near the frequency of the heavy-hole exciton resonance entanglement may reach maximum E N = 1 only in the case of linear polarization of the pump field and entanglement demonstrates interesting dependence on the orientation of the plane of polarization. Near the light-hole exciton resonance the most advantageous orientation of the ellipse of polarization is χ = 0, however, the direction along which the most emitted photons are emitted strongly depends on the ellipticity of the external excitation.

Optical response of Bose-Einstein condensate
Bose-Einstein condensation (BEC) is a phenomenon when at non-zero temperature the majority of particles occupy only a few states. This is in striking contrast to a distribution prescribed by the classical theory, where it is governed by the Boltzmann exponent exp(−E/k B T) and significant difference in occupations may be expected only when the energy levels are sufficiently far away from each other ΔE/k B T ≫ 1. The Bose-Einstein condensation, as is well known from the standard textbook consideration of ideal Bose gases (see e.g. Chapter 12 in (Huang, 1987), where it is clearly shown how a condensate emerges during the transition to the thermodynamic limit), does not require such level separation and may as well appear in a system with continuous spectrum.
The effect of BEC is tightly connected to such highly unusual from the classical point of view phenomena as superconductivity and superfluidity, which enjoy detailed developed theories (Lifshitz & Pitaevskii, 2002) and still are inexhaustible sources of new questions. In contrast to superfluidity and superconductivity, for which experiments have taken the lead over theoretical considerations, experimental studies of BEC fall well behind the theory (Moskalenko & Snoke, 2000): being predicted in 1925 BEC was obtained in a laboratory only in 1995. The difficulty of observing BEC motivates the constant search for more optimal systems. With this regard the significant attention has been paid to excitons in semiconductors. The physics of transition of semiconductors into the condensate state is similar to the superconducting transition (in conventional superconductors). While the elementary excitations are electrons, at sufficiently low temperatures Cooper pairs of electrons play the important role, which are formed by phonon-mediated attraction. Pairs are no longer subject to the exclusion principle and, moreover, at low densities obey the boson commutation relation. Thus, they may undergo the transition into the BEC-state. With excitons in semiconductors a similar scenario may take place and may even be more favorable because the Coulomb interaction binds electron and hole instead of preventing them from forming a bound state as it is in superconductors.
The finite life time of excitons, however, makes it difficult to reach the condensate stateexcitons decay before the equilibrium is established. Therefore, recently indirect excitons in coupled quantum wells became the object of special interest Butov, Lai, Ivanov, Gossard & Chemla, 2002;Snoke et al., 2002). The electrons and holes are spatially separated in such a structure that leads to increased life-time. Recently another possibility, BE condensate formed by dark excitons, started to attract attention (Combescot & Leuenberger, 2009).
The obvious difficulty related to dark excitons is how to observe them. One of possible ways to test properties of dark excitons is to use indirect interband spectroscopy, which relies on dynamics of bright states modified by the presence of dark excitons.
The problem of the optical response can be approached along the same lines as in the previous sections. The Hamiltonian of light-matter interaction is treated as perturbation and using the interaction picture the exciton polarization of bright excitons P = B is found in terms of S operator, which produces the perturbation theory. Before we apply this ideology we need to revisit the notion of averaging in formulas containing ... . In Section 2.4 the initial state of semiconductor was taken to be vacuum. Here, however, we need to take into account that initially the system is in thermal equilibrium and therefore its state is given by density matrix ρ rather than by a vector of state. Thus for an operator O we need to consider In the case when the system has BE condensate this expression significantly simplifies because the main contribution results from the condensate states, which contain the macroscopic number of particles. This observation is formally expressed by the spectral decomposition for the density matrix ρ = ∑ n w n |ψ n ψ n |,w h e r e|ψ n are some orthogonal states and w n are their weights. The ratio of the weights of non-condensate and condensate states is small, w nc /w c ≪ 1, and thus the contribution of the respective terms in Eq. (38) can be neglected leaving us with We model the condensate state by a coherent state, which is obtained from vacuum by .H e r eλ and function χ are parameters of the state and will be determined later. It should be noted, however, that translational symmetry of the system implies that χ(x 1 , x 2 ) should posses the symmetry with respect to translations of both arguments χ(x 1 + a, x 2 + a)=e iP·a χ(x 1 , x 2 ) with some P. Clearly P can be eliminated in a moving frame and, hence, among states created by B † σ the smallest energy would be of those with P = 0. Thus we need to consider The spin states of hole and electron entering the pair creation operator B † σ are denoted by σ = {σ, s} and are such that the pair does not interact with the electromagnetic field. As has been shown above, there are two such states with σ − = {3/2, −1/2} and σ + = {−3/2, 1/2}. In the absence of external static magnetic field the energy of these states is the same thus we naturally have fragmented condensate approximately described by the density matrix ρ = (|ψ + ψ + | + |ψ − ψ − |)/2, where |ψ ± = D † σ ± |0 . Using this approximation in Eq. (39) we which reduces the problem of finding exciton polarization to the problem with pure initial state similar to analyzed in Section 2.4. The first and the second terms in Eq. (42)  Applying the diagrammatic representation of the perturbation series (as illustrated in Eq. (21) one can immediately see that the series will contain only the same diagrams as in Section 2.4. In contrast, if the condensate was made of bright excitons then all diagrams shown in Eq. (21) would contribute. For example, the first diagram, without the external field, would describe the radiative decay of excitons in the condensate. For dark condensate, however, this diagram turns to zero because cs ψ ± = cσ ψ ± = 0. Thus one can see that only diagrams with matching number of creation and annihilation operators of electrons or holes not containing in the condensate are not vanishing and one need to keep only diagrams with the vacuum line, where vacuum is understood for electrons with spins or for holes with spinσ.
Unitarity of the shift operator D + (λ) allows one to present the average as taken over vacuum Only electron and hole operators with spins s and σ are affected by this transformation. In order to find how they transform it is convenient to use the momentum representation, e.g. v σ (x)=( 2π) −d/2 dkv σ (k)e ik·x with d being the dimensionality of the problem. In this representation we have and thus

Solving the system of equations we find the transformation induced by the shift operator is equivalent to the Bogoliubov transformation
where In order to clarify the physical meaning of the transformation we present Eq. (46) in coordinate representation where E(k)=ǫ c − ǫ v − Ω + k 2 /2m e + k 2 /2m h . As the zeroth approximation we obtain α(k)=1a n dβ(k)=cφ(k),w h e r eφ(k)=2 2r 3 B /[π(r 2 B k 2 + 1) 2 ] is the Fourier transform of 1s-exciton states in 3d and r B is the exciton Bohr radius. The constant c is found from the "normalization condition": the electron (or hole) density is equal to the condensate density that is c 2 = n + (2π) 3 . In the dilute regime, when the condition holds we have β(k) ≪ 1 and hence λχ(k) ≈ cφ(k) thus completely defining the condensate state.
The major effect of the dark condensate on bright excitons is that the condensate changes the structure of excitons. In order to see this we consider the polarization of linear response given by Eq. (22) with transformed exciton operators. Let us consider for definiteness the case of right polarized external excitation, which is coupled to bright exciton with σ = {σ,s} = {−3/2, −1/2}. The time dependence of propagator Φ μ,ν (τ) is determined by where A(p, q)=α(q)β(q − p) − α(q − p)β(q) and The last term in Eq. (52) yields the modification of the energy through perturbation of the condensate. This term can be estimated to lead to relatively long beatings in the linear response and will be neglected in the following.
Similarly to the case of the dynamics of initially unperturbed semiconductor Eq. (52) provides the structure of elementary excitations given by the spectral decomposition of the kernel : its eigenfunctions define the modified exciton operators, which then should be used in expansion (9). We present results for the lowest energy state in Fig. 4. The exchange by hole with condensate leads to renormalizations of the mass of the hole and its interaction with electron within the bright exciton. As a result the observed exciton states exhibit a blue shift for the binding energy and increasing radius (see Fig. 4a and inset in Fig. 4b) as functions of the condensate density. In Fig. 4 we have plotted the absorption spectrum for different densities taking into account the variation of the exciton oscillator strength, which is determined by |φ(x = 0)| 2 ∝ 1/r 3 B (η).

Ab initio approach -time-dependent density functional theory
In this Section, we present an alternative, ab initio, approach to describe ultrafast processes, based on time-dependent density functional theory (TDDFT). (Runge & Gross, 1984) The main advantage of this approach compared to the Green's function many-body theory methods is its simplicity. Being an effective single-electron theory, TDDFT uses time-dependent single-electron density to describe the nonequilibrium response. The corresponding single-particle Kohn-Sham wave equation depends on one vector coordinate and one time variable, which gives a big advantage compared to the Green's function approaches, where in strongly nonequilibrium case one needs to take into account many-time Green's functions, and the truncation used to have a finite system of equations is often not well justified. Moreover, technically to solve this system of many-time Green's functions is a very complicated task since in general case there is no time translation invariance, and one cannot use Fourier transforms and needs to consider time variables on the complex Keldysh type time-contour (see for example (Kadanoff & Baym, 1962)). This analysis can consume a very large amount of time (see, e.g., Refs. (Freericks et al., 2006;Onida et al., 2002)).
Another difficulty in the Green's function case comes from the correct treating of the Coulomb interaction. In the case of pulses shorter than the Coulomb scattering time, one cannot take these effects into account within a phenomenological scattering time parameter, like in the Boltzmann equation approach. Thus, one needs to take into account the electron correlation effects more accurately in this case, which is a very complicated problem in the majority of cases. In TDDFT, provided one has the correct exchange-correlation (XC) potential, these effects are taken into account exactly in the framework of a simple single-particle Kohn-Sham equation with relatively simple XC potential responsible for the electron correlation effects. Even though the form of such a potential in the case of materials containing transition metal atoms or atoms with f-electrons in the valence band is not solved yet, in the case of familiar semiconductor and molecular systems standard local-density approximation (LDA) and generalized gradient approximation (GGA) are often proved to be good approximations, like in description of single-electron excitations in molecular systems. (Elliott et al., 2009) One of the most important questions in TDDFT is its ability as an effective single-particle theory to describe multi-particle effects, including multiple-electron excitations and bound electron-hole states, like excitons and biexcitons. To be more specific, one needs to find the form of the corresponding XC potential to describe these effects. At the moment our knowledge about the structure of such potentials is rather limited. Though, recently significant progress in this direction has been made.
In the case of a weak and slow perturbation, the algorithm based on the many-body Bethe-Salpeter equation (BSE) to construct the XC potential able to describe excitonic effects was proposed in Refs. (Botti et al., 2004;Marini et al., 2003;Reining et al., 2002). It was shown also that in the linear response regime in frequency representation one can construct a pure TDDFT potential, the time-dependent optimized effective potential in exact-exchange approximation (XX-TDOEP), to describe excitonic effects in the optical absorption spectra. (Kim & Görling, 2002a;b) Despite the progress of the approaches mentioned above, these methods are rather complicated, which makes them difficult to apply in general nonlinear case of strong ultrafast perturbations. Moreover, the question whether they can describe higher order effects, like biexcitons or other nonlinear collective effects, including Bose-Einstein condensation of excitons remains open.
Recently, we proposed formally simple and rather general TDDFT approach based on the density-matrix representation, which is able to describe some of the phenomena mentioned above. (Turkowski et al., 2009;Turkowski & Ullrich, 2008;Turkowski et al., 2010) In this Section, we present the general formalism used in this approach. In particular, we present the system of the TDDFT equations for the density matrix elements, which is the TDDFT version of semiconductor Bloch equations, that includes nonlinear excitonic effects in all orders of magnitude. As we demonstrate, the solution of these equations in the case of excitons, and their generalization on the case of biexcitons, gives the corresponding energies of the bound states and wave functions, allows to analyze nonlinear effects like exciton-exciton interaction, which can be applied to the description of four-wave mixing experiments. Finally, we discuss the possibility to describe with the approach a highly nonlinear, collective effect of the Bose condensation of excitons.

General formalism
In TDDFT the time evolution of the system is studied by solving time-dependent Kohn-Sham equation: where the time-dependent Hamiltonian consists of the kinetic energy part (first term), the nuclear potential V nucl (r) for the electrons, the Hartree potential V H [n](r)= dr ′ n(r ′ )/|r − r ′ |, the XC potential V xc [n](r, t) and the external field term V(r, t). Generally speaking, in the case of an external electric field E(r, t), this field has to be included into the Kohn-Sham Hamiltonian through the standard substitution ∇→∇ − (i/c)A ext (r, t), and the scalar potential term ϕ ext (r, t),w h e r et h e external vector potential A ext (r, t) and ϕ ext (r, t) define the electric field: In the case of extended systems one needs to use the vector potential in order to describe the periodic system. Moreover, one should use the current-TDDFT with the macroscopic current that leads to the periodicity of the system. However, in the case when the field frequency is larger than the level spacing, the scalar potential V(r, t)=ϕ ext (r, t)= −E(t)r, which significantly simplifies the solution and which we shall use below, is a good approximation. (Schäfer & Wegener, 2002) The XC potential takes into account all electron exchange and correlation effects beyond the Hartree approximation. The Hartree and the XC potentials depend on the time-dependent electron charge density, which can be expressed in terms of the wave-functions Ψ i k (r, t), solutions of Equation (54), where k is the momentum and i is the band (with the energy ε v i k ) and all other quantum number indices, and ε F is the Fermi energy in the case of extended bulk system. The bandstructure and the Fermi energy can be obtained from the stationary solution of the Kohn-Sham problem (standard DFT theory (Kohn, 1999)). It is easy to generalize the corresponding problem on the case of finite systems. The simplest well-known approximation for the XC potential (the exchange only part) is the adiabatic LDA approximation: V xLDA (n(r)) = − (3/π) 1/3 n 1/3 (r). In the case of highly-correlated systems and strongly nonequilibrium processes more complicated forms for V xc [n](r, t) have to be used. In particular, in the case of the ultrafast processes the memory effects (dependence of the XC potential on the charge densities in all previous times) can play an important role.
In order to study the optical transitions in the system, it is convenient to express the time-dependent wave functions as linear combinations of the ground-state wave functions: where c v i l j k (t) are momentum and time-dependent complex coefficients, and l i = v i , c i are the valence and conduction band indices. One can solve the problem by studying the time dependence of these coefficients, however from the physical point of view it is more convenient to analyze the solution for the density matrix: . The elements of this matrix correspond are related to the probability of the optical transitions between different bands and their occupation. The density matrix satisfies the Liouville-von Neumann equation: the last equation the total polarization is the sum of polarizations for all possible interband transitions: is the corresponding bandgap and Ω cell is the volume of the unit cell.
For simplicity, in this Subsection we shall concentrate on the two-band case with 2 × 2density matrix with two independent matrix elements: the conduction band occupancy ρ cc k (t) and the polarization ρ cv k (t). The other two elements, the valence band occupancy ρ vv k (t) and the polarization de-excitation matrix element ρ vc k (t), can be found from the first two elements by using the conservation of particles equation ρ vv k (t)+ρ cc k (t)=1 and the definition of the matrix elements, which gives ρ vc k (t)=ρ cv * k (t). The independent TDDFT matrix equations have the following form: They correspond to the many-body theory semiconductor Bloch equations, (Haug & Koch, 2004) but in the TDDFT case the correlation effects are taken into account exactly without making the Hartree-Fock truncation.
In order to get a better feeling of the correspondence between both theories, it is instructive to compare both systems of the equations. Applying the same expansion of the wave function in terms of the stationary wave functions and writing down the corresponding Liouville-von-Neumann equation for the density matrix in the case of the Hartree-Fock equation one can obtain familiar set of the semiconductor Bloch equations (in linear expansion in the polarization, see below): Comparizon of the systems of equations (63), (64) and (66), (67) suggests that the electron-hole interaction term dq (2π) 3 V(k − q)ρ cv q (t) resposible for the Rydberg series of the bound states in the Bloch equation case is contained in the nonlinear functional of the polarization ρ cv q (t) of the matrix element V cv xck (t) (through the polarization that enters in the charge density, Eq. (61)). We have analyzed the optical absorption spectrum in the case of several XC potentials and found that in some cases the spectrum contains an excitonic peak below the conduction band edge. (Turkowski & Ullrich, 2008) In particular, we have found that the optical absorption spectrum demonstrates a pronounceable excitonic peak when the XC kernel contains the Coulomb singularity 1/q 2 , like in the case of the KLI and Slater potentials. (Krieger et al., 1992) This result is in agreement with Kim and Görling (Kim & Görling, 2002a;b) who showed that in the translational-invariant systems in order to have the excitonic peaks one needs to have such a singularity. On the other hand, it was found that the standard LDA and GGA potentials are "too weak" to produce the peaks.
In order to get a deeper understanding of the structure of the XC kernels necessary to produce the excitonic bound state one can analyze the linearized TDDFT equation for the polarization which corresponds to the Wannier equation −(∇ 2 /2m r ) − (1/ǫr) φ(r)=Eφ(r) for the exciton eigenenergies and eigenfunctions (m r is the reduced electron-hole effective mass and ǫ is the static dielectric constant of the material). (Wannier, 1937) The solution of the last equation demonstrates a Rydberg series of the excitonic binding energies qualitatively described by the Elliott formula (Haug & Koch, 2004).
The corresponding TDDFT Wannier equation can be obtained by linearizing equation (64): with the effective electron-hole interaction (in the momentum representation). The corresponding real-space equation can be obtained after the Fourier transforms ρ(R, ω)=∑ k e −ikR ρ cv k (ω) and V eh (R, R ′ , ω)= ∑ k,q e −ikR F kq (ω)e iqR ′ ,w h e r eR is a direct lattice vector. Since the excitonic wave function extends over many lattice sites, one can consider R as a continuous variable. In this case the TDDFT Wannier equation takes the following form is the reduced mass and E KS g is the KS band gap). The solution of the last equation gives the exciton eigenfunctions ρ(R, ω) and eigenenergies, which are defined by a nonlocal, frequency-dependent electron-hole interaction V eh (r, r ′ , ω). This interaction is defined by the XC kernel. The analysis of the solution in the case of the of LDA kernel, shows again that it is too weak to produce bound states. On the other hand, it was found that a phenomenological local kernel f contact xc (r, r ′ )=−Aδ(r − r ′ ) (A is a positive constant), and a long-range kernel with the Coulomb singularity ∼ 1/q 2 in the momentum space, f LRC xc (r, r ′ )= −α/4π|r − r ′ | (α is an adjustable parameter, which might be interpreted as an effective inverse screening), give correct lowest excitonic energy at proper choice of the parameters A and α. Thesevaluesareofthesameorderofmagnitudeforonegroupofsemiconductors-zincblende or wurtzite, but can be one or two orders of magnitude different for the semiconductors from different groups. (Turkowski et al., 2009) This suggests that these kernels, similarly to the LDA case, might be defined by system parameters, in particular the electron density and the volume of the unit cell. The generalization of this formalism, based on two-particle density matrix, was used to study the biexcitonic binding energies. (Turkowski et al., 2010)

Biexcitons
The above formalism can be generalized to the case of multiple excitations, in particular on the case of biexcitons, correlated double electronic excitations. (Turkowski et al., 2010) In principle, one can obtain double excitations and possibly coupled (biexcitonic) states within single-particle TDDFT in the case of non-adiabatic XC kernel. In this case, nonlinear Casida equation for the eigenenergies will have extra solutions in addition to single-particle excitations. In this Subsection, we analyze how one can obtain biexcitonic states within adiabatic approximation since this case corresponds to a transparent biexciton eigenproblem. In order to find such an approach, one may use the natural orbital (NO) representation for the stationary electron eigenfunctions. (Giesbertz et al., 2008;Pernal et al., 2007) In this case multi-particle excited states can be described by elements of one-and two-electron density matrices, defined as where Ψ is the N-particle wave function and x i =( r i , s i ) denotes the space coordinate and spin index. (Giesbertz et al., 2008;Pernal et al., 2007) In principle, all ground state properties can be obtained from the the single-particle matrix γ(x 1 , x ′ 1 ), due to one-to-one correspondence between the matrix and the ground state many-body wave function Ψ (the density matrix functional theory generalization of the Hohenberg-Kohn theorem (Gilbert, 1975)). Though to study the excited states the two-electron density matrix is necessary. We shall concentrate on an effective two-electron theory described by the Hamiltonian H(r 1 , r 2 , t)=ĥ ad (r 1 , t)+ĥ ad (r 2 , t)+w[n 2 ](r 1 , r 2 , t), whereĥ ad is the single-particle TDDFT Hamiltonian (55). In order to have biexcitonic states in the adiabatic approximation, one can introduce an effective two-particle interaction w[n 2 ](r 1 , r 2 , t) which depends on two-particle density n 2 (r 1 , r 2 , t)=Ψ * (r 1 , r 2 , t)Ψ(r 1 , r 2 , t). Similar to the excitonic case, one can expand the two-electron wave-function in terms of the NOs χ k (r), which in the singlet case gives Ψ(r, r ′ , t)=∑ k,l C kl (t)χ k (r)χ l (r ′ ),whereC kl (t) is a symmetric matrix with respect to quantum number indices k abd l. In this case, it is easy to show that where γ kl (t)=2 ∑ m C km (t)C * T lm (t) and Γ klmn (t)=2C kl (t)C * mn (t). Some of the elements of the last two matrices are proportional to the excitonic and the biexcitonic wave functions (with indices cv and ccvv in the notations of the previous Subsection, correspondingly). Matrix elements C kl (t) satisfy the following equation of motion: with the initial condition C kl (t = 0) ∼ δ kl . The matrix elements h kr (t) are defined above and w klmn (t)= dr 1 dr 2 χ * k (r 1 )χ * l (r 2 )w[n 2 ](r 1 , r 2 , t)χ m (r 1 )χ n (r 2 ).
From equation (76) one can obtain the equations for γ kl (t) and Γ klmn (t): The solution of the eigenproblem for b ± nm,q gives one the biexcitonic eigenvectors and eigenenergies. In the last equation H ± nm,n ′ m ′ ,qq ′ are functionals of the excitonic functions and interaction elements F kk ′ and G kk ′ pp ′ (see Ref. (Turkowski et al., 2010)).
We tested the formalism in the case of several semiconductors by using phenomenological local one-particle kernel f contact xc (r, r ′ ) and two-particle "contact biexciton" kernels g local 1 (r, r 1 , r 2 )=−C 0 A 1 δ(r − r 1 )δ(r − r 2 ) and g local 2 (r, r ′ , r 1 , r 2 )=−A 2 δ(r − r ′ )δ(r − r 1 )δ(r − r 2 ). It was found that indeed the TDDFT can describe the biexcitonic states in the adiabatic approximation.

Nonlinear effects
It is straightforward to extend the formalism developed above to the nonlinear case, including dynamical exciton-exciton interaction and memory effects. These processes play an important role in the case of ultrafast processes, including four-wave mixing experiments. So far, these nonlinear effects were studied only in the framework of many-body effective models. In most cases, the problem was analyzed by solving a third-order polarization equation. Beyond the importance of developing the TDDFT approach to describe the ultrafast processes, there is another important reason for this. Namely, from the experimental data one can learn about the non-adiabatic structure of the XC kernels, since our knowledge on the non-adiabatic kernels is much more limited comparing to the static adiabatic case. Below, we analyze some possible types of the response of the system by taking into account the memory effects and by using the known asymptotic limits of the XC kernels at low and high frequencies, and compare qualitatively the TDDFT results with the corresponding phenomenological many-body solution.
From equation (64) one can obtain the system of equations for the first and the third order polarizations: where and β kqpp ′ (t, t ′ , t ′′ , t ′′′ ) is a sum of matrix elements of f XC and its first two derivatives with respect to the particle density, similar to Eqs. (86) and (87). In the case of four-wave mixing experiments the system of equations (84), (85) can be solved by solving first the linear equation, and then the nonlinear effects can be found by solving Eq. (85). To study nonlinear effects one can also analyze the approximate effective third order equation for the total polarization which corresponds to the system (84), (85): similar to the many-body equation analyzed in Refs. (Ostreich et al., 1995;1998). In the last equation, Ω is the Rabi frequency and n c is the maximum density corresponding to the Pauli blocking term (we neglect the momentum variable below).
As it follows from Eq. (89), the nonlinear time-dependent effects are defined by the memory function F(t − t ′ ), which depends on the non-adiabatic part of the XC kernel.
In the many-body approach, the memory function usually depends on the exciton-exciton correlation function, which is difficult to find, so in this case a phenomenological approach has to be used. For example, as it was proposed in Ref. (Ostreich & Sham, 1999), in the case of slowly varying polarization the memory term can be approximated by so the time-correlation effects are defined by the function g(t)= t 0 F(t ′ )dt ′ .T h i sf u n c t i o n can be expressed in terms of the spectral density ρ(ω): g(t) ∼ ∞ 0 dωρ(ω)ω −1 e −iωt .I nt h e low-frequency limit, which defines the long-time asymptotic behavior of the system, the spectral density can be approximated by a power low-function ρ(ω) ∼ ω α .T h i sf u n c t i o n defines the dissipation processes in the system, i.e. the role of the environment (other excitons) on the behavior of given exciton. In the cases when α is smaller, equal or larger than 1, the dissipation is called "sub-ohmic", "ohmic" and "super-ohmic" (Caldeira & Leggett, 1983). Since the spectral function must decay at large frequencies, the general form of the spectral density was approximated by where ω F is the frequency scale and A is the normalization constant. Thus, the memory function can be approximated by (Ostreich & Sham, 1999) In the case of an one-dimensional model for excitons, the authors of Ref. (Ostreich & Sham, 1999) found α = 1.
One can in principle construct an f XC such that the functions on the right hand sides of the Eqs. (89) and (92) are equal. In this case, the equation for polarization will coincide with the many-body equation, and the solution in both cases will be the same. However, as we show below due to some constraints on the frequency-dependent f XC , the last equation should be corrected. Namely, it is known that the exact asymptotic of the XC kernel at large frequencies is f XC ∼ a + bω −2 (van Leeuwen, 2001). In the case of low frequencies, the information about the exact behavior of f XC is more limited. In particular, it is known that it can have poles in the case of finite system in the discrete part of the spectrum. To get an idea about possible frequency-dependence of the XC kernel for all ranges of frequencies one can consider the case of the homogeneous electron gas, when f XC (ω → 0) → 0, f XC (ω → ∞) → ω −3/2 (Marques & Gross, 2003). From these results one can suggest the following rather general form for the non-adiabatic part of the XC kernel: where α is of order of 1, and β = 2, though the case β = 3/2 is also worth of special attention. Below we solve Eq. (88) in the case of different values of α and β = 1.5 and 2. We approximate the XC kernel in the following way: n(r, t)] is the adiabatic part and f XC (t − t ′ ) the last term is defined in Eq. (93). Substitution of this expression into Eq. (89) leads to the following form of the memory function: F(t)=A ∞ 0 dωω α [1 +(ω/ω F ) α+β ] −1 e −iωt , where A is the integral over the derivative of the adiabatic part with respect to the particle density multiplied by the static wave functions (see Eq. (89)).
We analyze qualitatively possible solutions of Eq. (88) by considering two characteristic cases: an approximate time-evolution of the excitonic density and collective excitations in the case of two different memory functions: the exponentially decaying kernel (91) and the TDDFT-type algebraically decaying kernel (93). The time-dependence of the excitonic density n(t) can be obtained from the equation for polarization by using the ansatz P 0 = n(t)exp(iφ), which gives n(t) n(0)[1 + n(0)Re t 0 g(t ′ )dt ′ ] −1 . One can show that the equilibration takes much longer time in a more realistic case of the TDDFT spectral function, comparing to the exponential one (see Fig. 5).
One can analyze the collective excitations in the excitonic system by separating the slow and fast components of polarization, P = P 0 + P 1 , so the fast component satisfies the following approximate equation: ∂P 1 /∂t = −|P 0 | 2 iβP 1 + t 0 F(t − t ′ )P 1 (t ′ )dt ′ . The eigenvalues of this equations can be found numerically from the corresponding equation in frequency representation: (for details, in particular for the normalization of the function ρ and the spectral sum rule for F(t), see Ref. (Ostreich & Sham, 1999)). It is possible to show that the number of possible collective modes increases from 0 to 2 with the exciton density increasing (the zero-frequency  5. The time-dependence of the exciton density in the case of β = 2" (a) and exponential spectral function (b). We have used the parameters for the 3D model of semiconductor from (Ostreich & Sham, 1999): β = 52ω x /3, F(0)=14ω 2 x , ω x = 6.7meV.Thetimeisgiveninunits of 1/ω x and n(0)=0.1. mode corresponds to the Goldstone mode). Since at small frequencies the change of the right hand side of Eq. (94) with frequency is faster in the case of power spectral function, the critical value for the density above which there are collective oscillations is lower in this case. Also, the corresponding energies for these oscillations are lower in the case of power spectral function.
Finally, similar to the many-body case, one can analyze the two-dimensional Fourier spectrum of the system by taking into account memory effects. It is possible to show that the presence of the memory function in Eq. (88) can not only result in a shift of the excitonic peak in the spectrum but also lead to coupled exciton-exciton states. The detailed results of these studies will be published in the nearest future.
To summarize, we have shown that our TDDFT approach for excitons, despite being at the early stage of the development, shows to be a very promising and powerful method that can be used in many applications, in particular in studies of ultrafast processes.