text:
# Low-temperature kinetics of exciton-exciton annihilation of weakly localized one-dimensional Frenkel excitons

## I Introduction

Exciton-exciton annihilation is an important process that strongly influences the optical and opto-electronic properties of materials at high excitation densities. In particular, exciton-exciton annihilation affects the nonlinear and lasing properties of organic systems, such as J-aggregates and polymer films.  Though the importance of this process is well-recognized, the microscopic understanding of the annihilation kinetics is still rather poor. This holds especially under the conditions of strong exciton delocalization and (or) low temperature, where the usual bi-molecular theory of exciton-exciton annihilation is expected to break down. The aim of this paper is to study the annihilation kinetics in weakly disordered one-dimensional Frenkel exciton systems, where the exciton coherence size can be considerable (tens of lattice units). This study is of relevance to the optical properties and exciton dynamics in J-aggregates and molecular antenna systems.

The standard approach to describe the kinetics of exciton-exciton annihilation relies on the bi-molecular rate equation, in which it is assumed that the effective annihilation rate is proportional to the exciton density. This equation reads

$$\dot{n}=\gamma n\alpha n^2,$$
(1)
where $`n`$ represents the average exciton density, understood here as the number of excitations per molecule, $`\gamma `$ is the single-excitation (radiative and non-radiative) relaxation rate, and $`\alpha `$ is the co-called annihilation constant having here the dimension of 1/time. The effective rate of exciton-exciton annihilation in the system is then indeed proportional to the average exciton density $`n`$ and is given by $`\alpha n`$. The solution to Eq. (1) reads

$$n=\frac{\gamma n_0}{\gamma e^{\gamma t}+\alpha n_0\left(e^{\gamma t}1\right)},$$
(2)
where $`n_0`$ is the initial population of excitations. If the annihilation dominates the single-exciton relaxation ($`\alpha n_0\gamma `$) the excitation population decreases according to a hyperbolic (non-exponential) law

$$n=\frac{n_0}{1+\alpha n_0t}.$$
(3)
The typical picture that one commonly has in mind when modeling the annihilation process as is done in Eq. (1), is as follows. First, it is usually understood that, as a result of strong disorder and (or) high temperature, the Frenkel excitons represent, in fact, molecular excitations.  Next, it is assumed that the excitations (diffusively) move over the system. If the diffusion rate is large compared to the rate of nearest-neighbor annihilation, denoted by $`w_0`$, two excitations annihilate each other (by fusing into one high-lying molecular excitation that quickly loses its energy by vibrational relaxation) when they have reached neighboring molecules. In this case $`\alpha =w_0`$. On the other hand, if $`w_0`$ dominates the diffusion rate, the annihilation event may occur at a distance large compared to the the nearest-neighbor separation. In this case, the annihilation constant is determined by a convolution of the annihilation rate with the pair correlation function of two excitations,  and may in principle depend on time.

The bi-molecular rate description has limitations, which become rather important at low temperatures. Obviously, the bi-molecular approach does not account for the fact that the excitons in J-aggregates are generally quite strongly spread, with typical coherence lengths of several tens of molecules.  Even at room temperature, this length is of the order of 10 molecules,  owing to the large inter-molecular excitation transfer interaction in these systems. It is an open and interesting question whether this finite length may be accounted for by interpreting $`\alpha `$ in Eq. (1) as an effective annihilation constant. In this paper, we will address this question for low-temperature exciton systems.

We will study the kinetics of exciton-exciton annihilation of one-dimensional excitons that are weakly localized by static disorder. We will use the same framework as was done in Refs.  to calculate annihilation rates. The new element of the present paper is to use these rates to follow the kinetics of annihilation. An important step in describing the annihilation of extended excitons, is to distinguish between inter- and intra-segment annihilation.  The rationale for this distinction is as follows. As appears from numerical simulations of disordered exciton chains, the exciton states residing close to the bottom of the exciton band (the region that dominates the optical response) can be classified into groups of a few (two or three) states. The states within each separate group are all localized on the same segment of the aggregate, with a typical size $`N^{}`$ (often referred to as the number of coherently bound molecules), while the segments corresponding to different groups do not overlap. In fact, it turns out that the two or three exciton states within each such group are very similar in structure and energy to the lowest two or three states that exists on an ordered chain of length $`N^{}`$ In particular, the lowest state of such a group has a wave function spread over the segment without nodes and can be interpreted as the local ground state. The next higher lying state of the group has a well-defined node and looks like a first local excited state, etc. The energy difference of the local ground and first excited states agrees well with that of a perfect chain of length $`N^{}`$.

Obviously, to describe exciton-exciton annihilation, one should consider at least the two-exciton states. As is well known, one-dimensional Frenkel excitons are weakly interacting fermions (see Refs. ). Thus, the wave functions of states with two excitons can be composed of Slater determinants of two one-exciton wave functions. Under the condition of weak localization, two different types of two-exciton states then appear: (i) those with two excitons belonging to the same localization segment, and (ii) those with the two excitons localized on different segments. This immediately leads to the distinction of intra-segment and inter-segment annihilation as fundamentally different annihilation channel.

This paper is organized as follows. In Sec. II, we present our microscopic model of annihilation, express the annihilation rate in terms of the basic interactions and wave functions, and make the formal step towards the annihilation kinetics at low temperatures, where the diffusive motion of excitons towards each other may be neglected. In Sec. III we use the distinction between inter- and intra-segment annihilation to derive qualitatively analytical expressions for the low-temperature annihilation kinetics. A more detailed study is presented in Sec. IV, where we basically exactly solve the kinetics, formally defined in Sec. II, through numerical simulations. In Sec. V we summarize our findings and discuss the relevance to experimental low-temperature annihilation data.

## II Two-exciton annihilation model

### A Motivation

Under usual experimental conditions, only a small part of the localization segments on molecular aggregates are excited. For example, the authors of Ref.  estimated that in their experiments, one hundred molecules per aggregate were produced at the highest excitation power applied ($`0.98GW/cm^2`$). As a physical aggregate normally consists of $`10^4`$ molecules,  while the typical localization segment in their particular case counted 20 molecules, these authors concluded that less than one exciton was created per segment of localization (on average, one exciton per five segments). A simple consideration based on the Poisson distribution for the probability of finding an integer number of excitons per segment, shows that on the physical aggregate about 80 out of 500 segments are expected to be singly excited, while only 8 are doubly excited. Triply (and more) excited segments are almost absent. Bearing in mind that excitons, created on the same segment of an aggregate or on closely spaced separate segments, will annihilate first, we conclude that a two-exciton model of annihilation seems to be quite reasonable as a first step.

### B Rate of two-exciton annihilation

As a working model, we adopt a linear chain of $`𝒩`$ three-level molecules as depicted in Fig. 1 The two lower molecular states, denoted ”0” and ”1”, are assumed to form an exciton band, as a result of a sufficiently strong resonant dipole-dipole inter-molecular coupling. The corresponding exciton Hamiltonian, taken in the nearest-neighbor approximation, reads

$$H_{ex}=U\underset{n=1}{\overset{𝒩1}{}}(b_{1n}^+b_{1,n+1}+b_{1,n+1}^+b_{1n})+\underset{n=1}{\overset{𝒩}{}}E_{1n}b_{1n}^+b_{1n},$$
(4)
where $`U<0`$ is the nearest-neighbor hopping integral chosen to be negative, as is the case for J-aggregates, and $`b_{1n}^+(b_{1n})`$ denotes the Pauli creation (annihilation) operator of the first excited state of molecule $`n`$ (the state with all the molecules in their ground states serves as the vacuum state $`|0`$ and has zero energy). The second term in Eq. (4) represents the Hamiltonian of non-interacting molecules, in which $`E_{1n}=E_1+\mathrm{\Delta }_n`$ is the energy of the first excited state of molecule $`n`$ with $`E_1`$ and $`\mathrm{\Delta }_n`$ being, respectively, the mean value of the energy and a static random offset. The latter simulates on-site (diagonal) disorder and results in localization of the excitonic states. Fluctuations of the nearest-neighbor coupling are neglected. It will be assumed that $`\mathrm{\Delta }_n`$ is distributed uniformly within the interval $`[\mathrm{\Delta },\mathrm{\Delta }]`$, so that the typical magnitude of the disorder is given by the standard deviation $`\sigma =\mathrm{\Delta }/\sqrt{3}`$. For $`\sigma U`$ the exciton eigenfunctions are localized within rather large segments of the chain: $`1N^{}𝒩`$ Throughout this paper, we will assume that this condition holds.

The high-lying electronic-vibrational molecular term, depicted as ”2” in Fig. 1, serves as the intermediate state through which annihilation occurs.  We consider one of the electronic-vibrational levels to be resonant with the two-exciton optical states and to undergo an efficient phonon-assisted relaxation to the ground vibronic state (see Fig. 1). The annihilation process itself consists of transferring the energy of two excitons to the high-lying molecular term. Assuming this step to occur due to the resonant dipole-dipole inter-molecular interaction, we may write the corresponding Hamiltonian as follows

$$H_a=\frac{1}{2}\underset{m,n=1}{\overset{𝒩}{}}\frac{V}{|nm|^3}b_{1n}b_{1m}(b_{2n}^++b_{2m}^+)+h.c.,$$
(5)
where $`V`$ is the matrix element of the annihilation operator for nearest neighbors and $`b_{2n}^+(b_{2n})`$ denotes the Pauli creation (annihilation) operator of the high-lying state of molecule $`n`$. The operator (5) annihilates the two excitations occupying molecules $`m`$ and $`n`$ and excites one of these molecules in the high-lying state. The implication of $`H_a`$ for third-order nonlinear optics of J-aggregates has been studied in Ref. .

In accordance with our arguments in Sec. II A, we will assume that not more than two excitons are created by the pump per linear chain. Moreover, we will assume that $`|V|`$ is small compared to the rate, $`\mathrm{\Gamma }`$, of phonon-assisted relaxation in the high-lying molecular state. We may then use perturbation theory to calculate the rate of annihilation. Moreover, the back process (exciton fission) can then be neglected. The resulting expression for the rate of exciton-exciton annihilation starting from the two exciton eigenstate $`|\mu \nu `$ ($`1\mu <\nu 𝒩`$) is simply given by the ”golden rule”

$$w_a^{\mu \nu }=\frac{2\pi }{\mathrm{}}\rho (E_f)\underset{n=1}{\overset{𝒩}{}}|2n|H_a|\mu \nu |^2,$$
(6)
where $`\rho (E_f)`$ is the density of final states (hereafter replaced by $`1/\mathrm{\Gamma }`$). Because of the fermionic nature of one-dimensional Frenkel excitons, one can compose the two-exciton eigenfunctions as Slater determinants of the one-exciton eigenfunctions:

$$|\mu \nu =\underset{m=1}{\overset{𝒩}{}}\underset{n<m}{\overset{𝒩}{}}\psi _{\mu m;\nu n}|1m,1n,$$
(8)
$$\psi _{\mu m;\nu n}=\phi _{\mu m}\phi _{\nu n}\phi _{\mu n}\phi _{\nu m},$$
(9)
where $`\{\phi _{\nu n}\}`$ are the eigenfunctions of the one-exciton problem

$$\underset{m=1}{\overset{𝒩}{}}H_{ex}^{nm}\phi _{\nu m}=E_\nu \phi _{\nu n}.$$
(10)
Here, $`H_{ex}^{nm}=1n|H_{ex}|1m`$ and $`E_\nu `$ is the eigenenergy of the one-exciton state $`\nu `$. Substituting Eq. (9) into Eq. (6) one obtains

$$w_a^{\mu \nu }=\frac{2\pi V^2}{\mathrm{}\mathrm{\Gamma }}\underset{m=1}{\overset{𝒩}{}}\left[\underset{n=1}{\overset{𝒩}{}}^{}\frac{\psi _{\mu m;\nu n}}{(mn)^3}\right]^2,$$
(11)
where the prime denotes that $`nm`$. In particular, for a dimer $`(𝒩=2)`$ only one two-exciton state exists and its annihilation rate is given by

$$w_a^{12}=w_0=\frac{4\pi V^2}{\mathrm{}\mathrm{\Gamma }}.$$
(12)
In order to arrive at Eq. (11) we used the fact that Frenkel excitons are non-interacting fermions whenever the nearest-neighbor approximation is used for the hopping integrals. They become interacting quasi-particles when including the coupling to far neighbors. The importance of the latter can be estimated through the changes which the long-range terms produce in the density of exciton states. It is known that in one-dimensional aggregates, the long-range dipole-dipole interactions shift the exciton band bottom by approximately 20% compared to the nearest-neighbor model.  The smallness of this shift suggests the corrections due to long-range interactions to be of a perturbative nature. Indeed, for few-particle states, Frenkel excitons are weakly interacting (well-defined) fermions, despite the long-range coupling.

Furthermore, as follows from the results of both numerical simulations  and theoretical estimates  of the linear optical properties of disordered Frenkel chains, the oscillator strengths of the optical transitions near the lower band edge grow by approximately a factor of 2.5 due to the long-range dipole-dipole interactions. This results from a larger extension of the optically active exciton states in the exact dipole-dipole model compared to the nearest-neighbor model (at a fixed disorder strength). In principle, this effect is not of a perturbative nature. It can, however, be included into the final formulas for the annihilation rates (see below) by rescaling the number of coherently bound molecules $`N^{}`$ . The above arguments justify the nearest-neighbor framework as a reasonable approach to describe Frenkel excitons.

It is worth stressing, though, that the dipole-dipole interaction in the annihilation channel ($`H_a`$) can generally not be taken in the nearest-neighbor approximation, because the annihilation of two excitons localized on separate localization segments is determined by the coupling to far neighbors.

### C Low-temperature annihilation kinetics

As in the present paper we are mainly interested in the annihilation kinetics itself, we will neglect any other possible channel of population relaxation, such as the radiative transitions from two-exciton to one-exciton states (also acting towards lowering the exciton population), as well as a possible multi-phonon relaxation from the high-lying term to the one-exciton states (acting, on the contrary, towards raising again the exciton population). To calculate the kinetics of the exciton-exciton annihilation, we will assume that excitons are created by a resonant laser pulse, that is short compared to the inverse of the J-band width. Under these conditions, the initial populations of the two-exciton state $`|\mu \nu `$ is proportional to the corresponding oscillator strength $`F_{\mu \nu }`$ given by

$$F_{\mu \nu }=\left|\mu \nu |D^2|0\right|^2=\left(\underset{m,n=1}{\overset{𝒩}{}}\psi _{\mu m;\nu n}\right)^2,$$
(13)
where $`D=_{n=1}^N(b_{1n}^++b_{1n})`$ is the chain’s dipole operator (the chain length is assumed to be smaller than the emission wavelength).

After the initial creation process, excitons may in principle move over the chain, At low temperature, however, the possibility to move is very restricted and the optically excited localized Frenkel excitons are practically immobile.  The reason is that at low temperature ($`T`$ width of J-band), an exciton created in one of the local ground states may move to an other similar state only when the latter has an energy lower than the former. The typical energy offset between the local ground states is of the order of the width of their energy distribution (i.e., the width of the J-band). Therefore, after one jump the exciton typically resides in the tail of this distribution. The number of states with still lower energy then drastically reduces, giving rise to a strong increase of the mean distance to such lower energy states. In fact, already after one jump the exciton has a strongly suppressed chance to jump further, i.e., such a type of the spatio-energetic diffusion (towards lowering the energy) is stopped rapidly and does not yield a sufficient possibility for two excitons to approach each other and annihilate. It is worth noting that experiments also indicate the absence of such a diffusion, which would manifest itself in a red shift of the exciton emission spectrum relative to the absorption spectrum. The experimental data show that such a Stokes shift is either absent or has a small magnitude.  A similar situation occurs in glasses doped with rare-earth ions.

Following the above arguments, we will assume that two excitons annihilate from the positions where they have been created. Then, the time dependence of the population of the excited two-exciton states is given by

$$P_a(t)=2\underset{\mu \nu }{}f_{\mu \nu }\mathrm{exp}(w_a^{\mu \nu }t),$$
(14)
where $`f_{\mu \nu }=F_{\mu \nu }/_{\mu \nu }F_{\mu \nu }`$ and angular brackets denote an average over the disorder realizations. Note that we have normalized the population such that it equals 2 at $`t=0`$. The formula (14) will be the basis of our further analysis of the low-temperature annihilation kinetics.

## III Qualitative picture

Before carrying out numerical simulations, we first provide a qualitative analysis of Eq. (14). Following the arguments concerning the nature of the low-energy weakly localized one-dimensional states, we separate the summation in Eq. (14) into two parts: $`P_a(t)=P_a^{intra}(t)+P_a^{inter}(t)`$. The first part, $`P_a^{intra}(t)`$, includes all those terms $`\{\mu \nu \}`$ where the one-exciton states $`\mu `$ and $`\nu `$ are localized on the same chain segment (doubly excited segments). The second part, $`P_a^{inter}(t)`$, contains those terms where $`\mu `$ and $`\nu `$ reside on different segments. The fact that this distinction can only be made for low-energy states is no restriction, as anyhow these states are the ones that dominate the ground state to one-exciton and the one-to-two-exciton absorption spectrum. Using the picture of exciton states on a chain of effective length $`N^{}`$, one arrives at the intra-segment annihilation rate:

$$w_a^{intra}=\frac{5\pi ^6}{18(N^{}+1)^3}w_0,$$
(15)
where $`w_0`$ is given by Eq. (12) and the factor $`5\pi ^6/18270`$.

The second term, $`P_a^{inter}(t)`$, governs the annihilation of two excitons created on different localization segments and is characterized by the rate

$$w_a^{inter}=\frac{N^{}+1}{R^6}w_0,$$
(16)
where $`R`$ is the distance between the two excited segments. Note that the rate of the inter-segment annihilation scales linearly with $`N^{}+1`$. As the wave functions of both segments enter the expression for $`w_a^{inter}`$, one might intuitively expect a quadratic dependence on $`N^{}+1`$. However, only one of the two excited segments, namely the one that passes to the ground state $`|0`$ in the annihilation process, coherently contributes to $`w_a^{inter}`$ giving the factor $`N^{}+1`$ (so-called superradiant transition). The transition within the other excited segment occurs to the high-lying molecular state of each molecule. It is important to note that the latter events are summed incoherently, as is evident from Eq. (11), thus preventing the appearance of an extra power of $`N^{}+1`$.

Keeping in mind that $`RN^{}`$ as well as that the numerical factor in Eq. (15) is fairly large, one immediately deduces from Eqs. (15) and (16) that $`w_a^{intra}w_a^{inter}`$ provided that $`N^{}`$ is of the order of or larger than several units, which is the condition we will focus on in the simulations. From this, one expects that the kinetics of the exciton-exciton annihilation will consist of two distinct parts: a very fast short-time decay, described by $`P_a^{intra}(t)`$, and a very slow long-time tail for which $`P_a^{inter}(t)`$ is responsible. In order to estimate the relative weight of these two components, we will take into account the fact that the oscillator strengths of double excitation of a particular localization segment and excitation of two different segments are of the same order. Then, statistical arguments based on the Poisson distribution seem to be sufficient for making estimates. If the initial density of excitations equals $`n_0`$ (in our case, $`n_0=2/𝒩`$), the probability of finding a typical localization segment (of length $`N^{}`$) to be $`k`$-fold excited reads

$$p(k)=\frac{(n_0N^{})^k}{k!}e^{n_0N^{}}.$$
(17)
Therefore, the probabilities of double excitation of a typical localization segment and of excitation of two different segments are given by $`0.5\left(n_0N^{}\right)^2e^{n_0N^{}}`$ and $`2n_0N^{}e^{n_0N^{}}`$, respectively. Thus, the relative contribution of the shorter component in $`P_a(t)`$ is of the order of $`0.25n_0N^{}`$, which drops upon increasing the disorder strength $`\sigma `$ (or upon a corresponding decrease of $`N^{}`$) and grows with increasing the density of excitations $`n_0`$.

Now, we turn to a discussion of the character of the annihilation kinetics. We recall that in the bi-molecular model, the decay is non-exponential, as a result of the nonlinearity of the driving equation (1). We will argue that in our case, the kinetics is also non-exponential, which, however, does not result from a nonlinearity. Let us consider first the intra-segment channel of annihilation. Here, according to Eq. (14), the nonexponentiality is expected even in the absence of disorder (i.e., for a regular chain), because the annihilation rate $`w_a^{\mu \nu }`$ obviously depends on which exciton states are involved in the annihilation process (see the next Section). In the presence of disorder, there is an additional source for the non-exponential behavior of annihilation through the intra-segment channel. It originates from the fluctuations in the sizes of the localization segments, $`N`$. Equation (15) gives the typical magnitude of the intra-segment annihilation rate. In reality, $`N^{}`$ in Eq. (15) should be replaced by a fluctuating value $`N`$. Consequently, the sum over the overlapping states in Eq. (14) can be approximately substituted by an average over a distribution of $`N`$, $`G(N)`$. One then obtains

$$P_a^{intra}(t)2𝑑NG(N)\mathrm{exp}\left[w_a^{intra}(N)t\right].$$
(18)
As follows from numerical simulations, the standard deviation of N, $`[𝑑NG(N)(N\overline{N})^2]^{1/2}`$, is of the order of the mean, $`\overline{N}=𝑑NG(N)N`$, i.e., the distribution $`G(N)`$ is rather broad. Due to this fact, the resulting nonexponentiality is expected to be considerable. The numerical simulations presented in Sec. IV confirm this picture.

The origin of the non-exponential behavior of the inter-segment annihilation is two-fold. First, the corresponding annihilation rate, as in the previous case, depends on the size of the localization segment $`N`$ (see Eq. (16) ). Thus, fluctuations of the latter will affect the annihilation kinetics even for a fixed distance between excited segments. However, the character of the inter-segment annihilation is determined mostly by the strong dependence of the annihilation rate (16) on the distance between two excited segments, $`R`$. The annihilation kinetics caused by the fluctuations of $`R`$ is simply given by an average of the pair kinetics $`2\mathrm{exp}[w_a^{inter}(R)t]`$ over all realizations of R. In order to obtain an analytical estimate, let us assume that the density of excited segments is low (as is in our case), so that $`R`$ can be treated as a continuous stochastic variable. Then the annihilation kinetics is given by an integral similar to that in Eq. (18) with $`G(N)`$ replaced by a suitable $`R`$-distribution function, $`G(R)`$. We will adopt a uniform distribution for $`R`$, $`G(R)=N^{}/𝒩`$, assuming that the probability of finding a segment to be excited is equal to the inverse of the number of segments in the chain, $`𝒩/N^{}`$. In evaluating the integral, we will extend the integration over the entire positive axes, neglecting thus the minimal distance between two adjacent segments as well as the finiteness of the chain. Both approximations are justified at a low density of excitations. One thus arrives at

$$P_a^{intra}(t)2\left[1\frac{N^{}}{𝒩}_0^{\mathrm{}}𝑑R\left(1e^{w_a^{inter}(R)t}\right)\right]2\left[1\mathrm{\Gamma }\left(\frac{5}{6}\right)\frac{N^{}}{𝒩}(N+1)^{1/6}(w_0t)^{1/6}\right],$$
(19)
where $`\mathrm{\Gamma }(x)`$ is the Gamma-function. Equation (19) is correct provided that the second term on the right-hand site is less than unity. This holds in a very large time interval, determined by the inequality $`\mathrm{\Gamma }\left(5/6\right)N^{}(N+1)^{1/6}(w_0t)^{1/6}<𝒩`$. It follows from the stretched exponential behavior of Eq. (19) that further averaging of $`P_a^{intra}(t)`$ over the $`N`$-distribution will not change the character of the kinetics and results, in fact, in replacing $`N`$ by $`N^{}`$.

To conclude this section, we note that the (”artificial”) quantity $`P_a(t)`$ can be simply rescaled to the measurable magnitude - the density of excitons $`n(t)=P_a(t)/𝒩`$. Introducing the initial exciton density $`n_0=2/𝒩`$, we obtain

$$n^{intra}(t)n_0\left[1\frac{1}{2}\mathrm{\Gamma }\left(\frac{5}{6}\right)n_0N^{}(N^{}+1)^{1/6}(w_0t)^{1/6}\right],$$
(20)
## IV Numerical simulations and discussion

To study the annihilation kinetics in more detail, we have carried out numerical simulations for a chain of length $`𝒩=200`$. For such a length, the mean initial density of excitations is $`n_0=0.01`$. We calculated the one-exciton eigenfunctions $`\phi _{\nu n}`$ by diagonalizing numerically the Frenkel Hamiltonian (4) for a particular realization of disorder, and then composed two-exciton eigenfunctions according to Eq. (6). Using further Eq. (11) and Eq. (13), we computed the rate of annihilation, $`w_a^{\mu \nu }`$, and the oscillator strength, $`F_{\mu \nu }`$, for any two-exciton state $`|\mu \nu `$. Then Eq. (14) was used to evaluate the annihilation kinetics. The resulting kinetics was obtained by averaging over 20 realizations of disorder. An increase of this number did not lead to considerable changes in the calculated curves. As a time unit, we used $`w_0^1`$. The results of the simulations for different values of $`\mathrm{\Delta }/U`$ are depicted in Figs. 2-5 by thick solid lines.

Figure 2 represents the annihilation kinetics for a perfect chain ($`\mathrm{\Delta }=0`$), as well as a least-square fit by means of an exponential $`2\mathrm{exp}(w_0t/a)`$ (thin solid line), achieved at $`a=1.6310^4`$. The fit clearly demonstrates that the calculated curve can not be matched by a single exponential. This unambiguously means that not only the states with $`\mu =1`$ and $`\nu =2`$, having the largest oscillator strengths, contribute to the sum in Eq. (14), but the other states contribute a comparable amount. The time scale of the kinetics depicted in Fig. 2 qualitatively corresponds to that calculated by using Eq. (15) with $`N^{}`$ replaced by $`𝒩`$: at $`𝒩=200`$, one obtains $`w_a^{intra}3\times 10^5w_0`$.

We also tried to fit the numerical curve in Fig. 2 by means of the bi-molecular equation (3) taken in the form $`2/(1+bw_0t)`$. The best fit was achieved at $`b=1.0710^4`$ and is plotted in Fig. 2 by the dashed line. At first glance, it seems that the latter almost matches the numerical data except, maybe, at the initial stage. However, the fitting constant $`b`$, carrying, in fact, the meaning of the density of excitations (see the discussion presented in the Introduction), underestimates the real value $`n_0=0.01`$ by two orders of magnitude.

Figures 3 and 4 show the numerical results obtained for disordered chains with different degrees of disorder, $`\mathrm{\Delta }/U`$. In Fig. 3, we plotted the annihilation decay curves in a wide time interval, while Fig. 4 presents the initial stages of the annihilation process. From the numerical results presented in Figs. 3 and 4 several conclusions can be deduced. First of all, it is clearly seen that the entire kinetics indeed consists of two well-distinguished components: a fast short-time decay, becoming faster as the disorder is increased, and a very slow long-time tail. The weight of the faster component is smaller than that for the slower one and drops upon increasing the disorder strength $`\mathrm{\Delta }/U`$.

It is reasonable to relate these two components to the intra- and inter-segment channels of exciton-exciton annihilation, respectively, in accordance with the qualitative picture discussed in the previous Section. Indeed, let us use for the typical size $`N^{}`$ of a localization segment the well-known estimate

$$N^{}+1=\left(3\pi ^2\frac{U}{\sigma }\right)^{2/3}.$$
(21)
Recall that in our case $`\sigma =\mathrm{\Delta }/\sqrt{3}`$. Substituting Eq. (21) into (15), one gets

$$w_a^{intra}=\frac{5\pi ^2}{506}\left(\frac{\mathrm{\Delta }}{U}\right)^2w_00.1\left(\frac{\mathrm{\Delta }}{U}\right)^2w_0,$$
(22)
Equation (22) gives us the disorder scaling of the intra-segment annihilation rate. Accordingly, we arrive at the following estimates: $`w_a^{intra}4\times 10^3w_0,2\times 10^2w_0`$ and $`6\times 10^2w_0`$ respectively for $`\mathrm{\Delta }/U=0.2,0.4`$ and $`0.8`$. Indeed, these numbers qualitatively match the time scales of the fast components (see Fig. 4).

We also plotted in Figs. 3 and 4 the least square fits of the numerical data by means of the function $`2c(w_0t)^{1/6}`$ (thin solid line). One observes that at higher degree of disorder $`(\mathrm{\Delta }/U=0.8)`$, when the weight of the faster component is smallest, the fitting function fairly well follows the numerical curve over almost the entire time interval of decay. The value of the fitting constant $`c=0.140.15`$ is of the same order as the one deduced from the theory, Eq. (19), according to which it must be $`2\mathrm{\Gamma }(5/6)(N^{})^{7/6}/𝒩=2\mathrm{\Gamma }(5/6)(3\sqrt{3}\pi ^2U/\mathrm{\Delta })^{7/9}/𝒩0.24`$ (The discrepancy probably stems from neglecting $`N^{}`$ as a minimal separation in Eq. (19).) This unambiguously means that the inter-segment channel dominates the long-time part of the annihilation kinetics. It should be especially stressed that the bi-molecular fits, shown in Figs. 3 and 4 by the dashed lines, fail absolutely in the description of the numerical data.

In order to show the character of the decay (exponential or non-exponential) in the case of disordered chains, we depicted in Fig. 5 the log-plot of the calculated annihilation kinetics for $`\mathrm{\Delta }/U=0.2`$. As can be seen, neither of the two components shows an exponential behavior.

## V Summary and concluding remarks

In this paper we have studied the low-temperature kinetics of exciton-exciton annihilation of weakly localized one-dimensional Frenkel excitons using a two-exciton static model (immobile quasi-particles) with diagonal disorder. Our analysis leads to three main conclusions:

i - The entire kinetics consists of two well-distinguished components: a very slow long-time decay and a much faster short-time drop. The latter component becomes faster with higher degree of disorder. The weight of the faster component is much smaller than that of the slower one, and decreases with increasing disorder strength.

ii - Neither of these two components shows an exponential behavior.

iii - The usual bi-molecular theory fails in the description of the behavior found.

These findings are well-understood from the existence of two competing options for two excitons to annihilate. The slower component is driven by the annihilation of exciton states localized on different segments of the chain, while the faster one originates from the annihilation of doubly excited segments. Fluctuations of distances between two excitons and sizes of the localization segments explain the non-exponential nature of the slower and faster components, respectively.

It is worthwhile to estimate the typical rates of both annihilation channels for existing J-aggregates. In order to do this, we need information concerning the parameters $`U`$, $`V`$, and $`\mathrm{\Gamma }`$. For J-aggregates $`U1000`$ cm<sup>-1</sup> is quite typical.  Less information consists concerning the annihilation interaction $`V`$, but as we assumed it to be of dipolar origin it seems not unreasonable to take a value similar to $`U`$. This is, in fact, supported by semi-empirical calculations of higher molecular singlet states of pseudo-isocyanine (PIC) molecules.  These calculations indicate a molecular $`S_1S_2`$ transition that is similar in energy and oscillator strength as the $`S_0S_1`$ transition responsible for PIC’s well-known J-band. We thus take $`V1000`$ cm<sup>-1</sup>. Finally, we will take $`\mathrm{\Gamma }3000`$ cm<sup>-1</sup>, corresponding to a vibrational relaxation time in the $`S_2`$ state of about 10 fs. Using these numbers, we arrive at $`w_a^{intra}3\times 10^{16}N_{}^{}{}_{}{}^{3}`$ s<sup>-1</sup>. The corresponding estimate for the inter-segment annihilation rate taken for adjacent segments ($`R=N^{}`$) reads $`w_a^{inter}10^{14}N_{}^{}{}_{}{}^{5}`$ s<sup>-1</sup>.

At low temperatures, the quantity $`N^{}`$ is found to be of the order of several tens . Letting $`N^{}=20`$, as was reported in Refs.  and , we arrive at $`w_a^{intra}4\times 10^{12}`$ s<sup>-1</sup> and $`w_a^{inter}3\times 10^7`$ s<sup>-1</sup>. Note that the magnitude of $`w_a^{inter}`$ appears to be even smaller than the spontaneous emission rate of a single molecule, which typically is of a few times $`10^8`$ s<sup>-1</sup>. Certainly $`w_a^{inter}`$ is much smaller than the spontaneous emission rate for an exciton state. It is to be noted furthermore that raising $`N^{}`$ by a factor of 2 will reduce $`w_a^{inter}`$ by almost two orders of magnitude. From the above, an important conclusion can be deduced: the inter-segment channel of exciton-exciton annihilation is in fact ineffective at low temperatures, because the radiative relaxation is much faster. On the contrary, the intra-segment annihilation rate is fairly high and should be viewed as the unique way for two weakly localized excitons to annihilate at low temperature. However, since this process occurs only for doubly excited localization segments, it will affect the entire exciton population only if the number of doubly excited segments is high, i.e., at sufficiently high laser intensities.

We note that the separation into inter-segment and intra-segment annihilation channels was in fact concluded from transient absorption experiments on PIC J-aggregates at 20 K.  The part of the observed annihilation kinetics with a decay time of 200 fs may indeed be related to intra-segment annihilation, as is clear from our estimate for $`w_a^{intra}`$. Our estimate for $`w_a^{inter}`$ shows, however, that it is unlikely that the second component of the kinetics reported in Refs.  and , with a decay time of 1.5 ps, may indeed be ascribed to inter-segment annihilation (see also the discussion in Refs.  and ).

Our findings concerning the ineffectiveness of inter-segment annihilation provide us with a way to control the exciton-exciton annihilation at low temperature. Indeed, recall that the local ground and first excited states belonging to the same localization segment are separated by the energy offset $`E_2^{}E_1^{}=3\pi ^2U/(N^{}+1)^2`$, which is of the order of the J-band width.  For typical J-aggregates, the exciton radiative rate $`\gamma `$, representing the unique relaxation constant at low temperatures, is much smaller than this energy mismatch. Hence one may get a large number of localization segments to be singly excited by applying a field with Rabi frequency smaller than $`E_2^{}E_1^{}`$, but larger than $`\gamma `$. At the same time, none of the localization segments will be doubly excited. Therefore, under such conditions, a fairly large exciton population may be created in J-aggregates, without being affected by exciton-exciton annihilation.

###### Acknowledgements.

This work has been partially supported by INTAS (project No. 97-10434). G. G. K. and V. A. M. also acknowledge support from the German Federal Ministry of Education, Science, Research, and Technology within the TRANSFORM - programme (project No. 01 BP 820/7). V. A. M. acknowledges the University of Juväskylä, where this work was started, for hospitality.


format: text

text:
# HD-THEP-01-11 hep-ph/0103049 March 2001 Nonequilibrium time evolution of the spectral function in quantum field theory

## 1 Introduction

Current and upcoming heavy-ion experiments at RHIC and LHC have been an important motivation for an extensive theoretical study of kinetic or transport equations beyond the level of the classical Boltzmann equation. In a (nonequilibrium) environment the quasiparticle structure of the theory can change and for instance the presence of nuclear matter in heavy-ion collisions may affect the mass and width of vector mesons (see e.g.  for a discussion of the broadening of the $`\rho `$ meson). In order to give a consistent treatment these effects have to be incorporated in a kinetic description. Also in other areas of physics an extension of the classical Boltzmann equation to a more general quantum Boltzmann equation is needed. The theory of baryo- and leptogenesis requires the solution of a large set of Boltzmann equations with a proper inclusion of CP violating and off-shell effects in an expanding universe . In Ref.  a description of the dynamics of Bose-Einstein condensates using a quantum Boltzmann equation can be found.

It is well known  that a consistent setup of kinetic equations derived from a microscopic quantum field theory needs to incorporate two quite different aspects: the statistical information typically described in terms of one-particle distribution functions on the one hand, and the kinetic information about the spectrum of the theory, encoded in the spectral function, on the other hand. There is a possibility to independently approximate the statistical and the spectral aspects of the time evolution equations. In the derivation of the classical Boltzmann equation one neglects all nontrivial kinetic aspects of the spectral function, which is approximated by an on-shell $`\delta `$-function. A first-order gradient expansion is then sufficient to derive the Boltzmann equation, describing the evolution of the statistical information only, from the underlying field theory . It turns out that relaxing the assumption of on-shell quasi-particles makes the problem considerably more complicated. Indeed, an analysis of a self-consistent inclusion of a nonzero width has been taken up only relatively recently . Similarly, to go beyond leading order in a gradient expansion seems to be rather cumbersome .

Recently, a successful numerical solution of the nonequilibrium time evolution without the use of a gradient expansion or on-shell approximation has been presented for a scalar field theory in $`1+1`$ dimensions . Both the early-time behavior and the large-time physics of thermalization can be described remarkably well for quantum fields , as well as for the classical field limit put to the test by Monte Carlo methods and numerical integration . The evolution equations can be obtained from the three-loop expansion of the closed time path generating functional for two-particle irreducible ($`2PI`$) Green’s functions  and comprise the classical Boltzmann equation . The full solution of these equations allows to directly investigate the influence of off-shell effects on the time evolution or the validity of a gradient expansion.

In this paper we focus on the kinetic aspects of the evolution equations and compute the nonequilibrium time dependence of the spectral function without further approximation. The paper is organized as follows. In the following section we briefly review the Schwinger-Keldysh formalism for a nonequilibrium scalar field. We derive exact time evolution equations starting from the Schwinger-Dyson equation for the two-point function and use the three-loop $`2PI`$ effective action to find a consistent approximation to the self energy. In Sec. 3 we discuss the initial conditions used in this paper. In order to make contact with heavy-ion physics we consider as initial state a so-called “tsunami”: a particle number distribution peaked in momentum space. Such an initial state is reminiscent of two colliding high energy wave packets . We hasten to say that in our $`(1+1)`$–dimensional scalar field model direct applicability to phenomology is of course limited. We study in this paper in particular the time evolution of the Wigner transformed spectral function suitable for comparison with standard frameworks for kinetic equations. Details of the Wigner transform are given in Sec. 4. We present our numerical results in Sec. 5. In Sec. 6 we make a connection with equal-time methods and moment expansions employed before. Our conclusions are summarized in Sec. 7.

## 2 Time evolution equations

Spectral function. We consider a real, scalar quantum field theory with field operators obeying the Heisenberg equation of motion $`(\mathrm{}_x+m^2)\varphi (x)+\lambda \varphi ^3(x)/6=0`$. We are interested in the nonequilibrium time evolution of the spectral function $`\rho `$ given by the expectation value of the commutator of two fields

$$\rho (x,y)=i[\varphi (x),\varphi (y)]_{}.$$
(2.1)
where $``$ denotes the average over an initial density matrix. Nonequilibrium initial value problems can be formulated using the Schwinger-Keldysh or closed-time-path formalism . For a detailed discussion we refer to Ref. . In this formalism the theory is formulated along a contour $`𝒞`$ in the complex-time plane, as shown in Fig. 1. Expectation values are defined along this contour. The two-point function is usually written as

$$G(x,y)=G_>(x,y)\mathrm{\Theta }_𝒞(x^0y^0)+G_<(x,y)\mathrm{\Theta }_𝒞(y^0x^0),$$
(2.2)
with $`\mathrm{\Theta }`$-functions defined along the closed time path $`𝒞`$ and the Wightman functions defined as $`G_>(x,y)=\varphi (x)\varphi (y)=G_<(y,x)`$.444These two-point functions were also used in Ref. . For a real scalar field $`G_{>}^{}{}_{}{}^{}(x,y)=G_<(x,y)`$. In terms of these, the spectral function is given by

$$\rho (x,y)=i[G_>(x,y)G_<(x,y)].$$
(2.3)
The spectral function can be obtained from the propagator (2.2) directly by writing

$$G(x,y)=F(x,y)\frac{i}{2}\rho (x,y)\left[\mathrm{\Theta }_𝒞(x^0y^0)\mathrm{\Theta }_𝒞(y^0x^0)\right].$$
(2.4)
Here we have separated the propagator using the commutator of two fields, $`\rho (x,y)`$, and the anti-commutator of two fields, the symmetric or statistical two-point function

$$F(x,y)=\frac{1}{2}[\varphi (x),\varphi (y)]_+=\frac{1}{2}[G_>(x,y)+G_<(x,y)].$$
(2.5)
For a real scalar field $`F(x,y)`$ and $`\rho (x,y)`$ are real functions with symmetry properties $`F(x,y)=F(y,x)`$ and $`\rho (x,y)=\rho (y,x)`$. We use $`\rho `$ and $`F`$ as the basic quantities in this paper.

Evolution equations. The full propagator (2.4) obeys the Schwinger–Dyson equation

$$G^1(x,y)=G_0^1(x,y)\mathrm{\Sigma }(x,y)$$
(2.6)
with the proper self energy $`\mathrm{\Sigma }`$ and $`G_0^1(x,y)=i[\mathrm{}_x+m^2]\delta _𝒞^{d+1}(xy)`$. A differential equation for $`G`$ which may be used to study the time evolution for given initial conditions can be obtained from (2.6) using $`_𝒞d^{d+1}zG^1(x,z)G(z,y)=\delta _𝒞^{d+1}(xy)`$ and one finds

$$\left[\mathrm{}_x+M^2(x)\right]G(x,y)+i_𝒞d^{d+1}z\mathrm{\Sigma }^{(\mathrm{nonlocal})}(x,z)G(z,y)=i\delta _𝒞^{d+1}(xy).$$
(2.7)
Here we have separated $`\mathrm{\Sigma }`$ in a local part and a nonlocal part,

$$\mathrm{\Sigma }(x,y)=i\mathrm{\Sigma }^{(\mathrm{local})}(x)\delta _𝒞^{d+1}(xy)+\mathrm{\Sigma }^{(\mathrm{nonlocal})}(x,y).$$
(2.8)
The local part $`\mathrm{\Sigma }^{(\mathrm{local})}`$ corresponds to a mass shift and we find it convenient to write $`M^2(x)=m^2+\mathrm{\Sigma }^{(\mathrm{local})}(x)`$.

Since our interest focuses on the nonequilibrium dynamics of $`\rho `$ we rewrite (2.7) such that the spectral function appears explicitly as one of the dynamical degrees of freedom. Similar to the above decomposition for $`G`$ we can decompose the self energy

$$\mathrm{\Sigma }^{(\mathrm{nonlocal})}(x,y)=\mathrm{\Sigma }_F(x,y)\frac{i}{2}\mathrm{\Sigma }_\rho (x,y)\left[\mathrm{\Theta }_𝒞(x^0y^0)\mathrm{\Theta }_𝒞(y^0x^0)\right]$$
(2.9)
into two real functions, the statistical component of the self energy obeying $`\mathrm{\Sigma }_F(x,y)=\mathrm{\Sigma }_F(y,x)`$ and the spectral component $`\mathrm{\Sigma }_\rho (x,y)=\mathrm{\Sigma }_\rho (y,x)`$. These self energies are related to the usual $`\mathrm{\Sigma }_{>,<}`$ as

$`\mathrm{\Sigma }_\rho (x,y)`$ $`=`$ $`i[\mathrm{\Sigma }_>(x,y)\mathrm{\Sigma }_<(x,y)],`$
$`\mathrm{\Sigma }_F(x,y)`$ $`=`$ $`{\displaystyle \frac{1}{2}}[\mathrm{\Sigma }_>(x,y)+\mathrm{\Sigma }_<(x,y)].`$
We also note that the retarded and advanced propagators and self energies are related to the spectral components as $`G_R(x,y)=\mathrm{\Theta }(x^0y^0)\rho (x,y)=G_A(y,x)`$ and $`\mathrm{\Sigma }_R(x,y)=\mathrm{\Theta }(x^0y^0)\mathrm{\Sigma }_\rho (x,y)=\mathrm{\Sigma }_A(y,x)`$.

It is straightforward to express (2.7) in terms of the real functions introduced above. The resulting equations are

$`\left[\mathrm{}_x+M^2(x)\right]\rho (x,y)`$ $`=`$ $`{\displaystyle _{y^0}^{x^0}}𝑑z^0{\displaystyle 𝑑𝐳\mathrm{\Sigma }_\rho (x,z)\rho (z,y)},`$ (2.10)
$`\left[\mathrm{}_x+M^2(x)\right]F(x,y)`$ $`=`$ $`{\displaystyle _0^{x^0}}𝑑z^0{\displaystyle 𝑑𝐳\mathrm{\Sigma }_\rho (x,z)F(z,y)}`$ (2.11)
$`+{\displaystyle _0^{y^0}}𝑑z^0{\displaystyle 𝑑𝐳\mathrm{\Sigma }_F(x,z)\rho (z,y)}.`$
The form of these equations is exact and they are equivalent to the standard Schwinger-Dyson identity (2.6). The equations for $`\rho `$ and $`F`$ are explicitly real and causal. A generic feature is the presence of nonlocal “memory” integrals. The time integral in Eq. (2.11) starts at some initial time taken to be $`z^0=0`$. We note that due to the canonical commutation relations the spectral function obeys the equal-time properties

$$\rho (x,y)|_{x^0=y^0}=0,_{x^0}\rho (x,y)|_{x^0=y^0}=\delta ^d(𝐱𝐲).$$
(2.12)
The time evolution of the spectral function is completely determined by (2.10) and (2.12).

Three-loop $`2PI`$ effective action. To solve the Eqs. (2.10) and (2.11) one has to find suitable approximation schemes for the self energies $`\mathrm{\Sigma }_F`$ and $`\mathrm{\Sigma }_\rho `$ which respect all symmetries and allow to describe the nonequilibrium early-time behavior as well as the large-time approach to equilibrium. Time reversal symmetry can be easily guaranteed by deriving the evolution equations from an effective action by a variational principle. A systematic way to achieve the latter is by using the loop expansion of the $`2PI`$ effective action , formulated along the closed time path . Recently, it has been shown that the equations in the three-loop approximation describe both the early-time behavior and the physics of thermalization . The self energies in this approximation , here expressed in terms of $`\rho `$ and $`F`$, are for a Gaussian initial density matrix given by

$`\mathrm{\Sigma }_\rho (x,z)`$ $`=`$ $`{\displaystyle \frac{\lambda ^2}{2}}\rho (x,z)\left[F^2(x,z){\displaystyle \frac{1}{12}}\rho ^2(x,z)\right],`$ (2.13)
$`\mathrm{\Sigma }_F(x,z)`$ $`=`$ $`{\displaystyle \frac{\lambda ^2}{6}}F(x,z)\left[F^2(x,z){\displaystyle \frac{3}{4}}\rho ^2(x,z)\right]`$ (2.14)
(cf. Ref.  for the extension to quartic initial density matrices). The effective mass term including the local part of the self energy reads

$$M^2(x)=m^2+\frac{\lambda }{2}F(x,x).$$
(2.15)
Eqs. (2.10, 2.11) with (2.132.15) form a closed set for $`\rho `$ and $`F`$ that is energy conserving and time-reversal invariant. We emphasize that the propagators in the “tadpole” diagram contribution in Eq. (2.15) and in the “setting-sun” contributions (2.13, 2.14) are full ones. The form of the self energies presented here (with free instead of full propagators) has also been found in a perturbative calculation in thermal equilibrium . We note from Eqs. (2.10, 2.11) that in local approximations, such as the Hartree approximation, derivable from the two-loop $`2PI`$ effective action, and leading-order large $`N`$ schemes, the evolution of the spectral function $`\rho `$ does not couple back to the evolution of the symmetric two-point function $`F`$.

Thermal equilibrium. It is instructive to consider the introduced quantities in thermal equilibrium. The two-point function $`G`$ in equilibrium is described by the same equation (2.7) if the closed time path is replaced by an imaginary path $`𝒞=[0,i\beta ]`$, with $`\beta `$ the inverse temperature . Since the equilibrium two-point functions depend only on the relative coordinates it is convenient to consider the Fourier transforms $`\rho ^{(\mathrm{eq})}(\omega ,𝐩)`$, etc. From the periodicity (“KMS”) condition for $`G(x,y)`$ and $`\mathrm{\Sigma }^{(\mathrm{nonlocal})}(x,y)`$ in imaginary time one finds the generic equilibrium relations relating for example $`G_>^{(\mathrm{eq})}(\omega ,𝐩)`$ and $`G_<^{(\mathrm{eq})}(\omega ,𝐩)`$ :

$$G_>^{(\mathrm{eq})}(\omega ,𝐩)=e^{\beta \omega }G_<^{(\mathrm{eq})}(\omega ,𝐩).$$
(2.16)
For the spectral and statistical components we employ here this translates into the following equilibrium relations

$`F^{(\mathrm{eq})}(\omega ,𝐩)`$ $`=`$ $`i\left(n_\mathrm{B}(\omega )+{\displaystyle \frac{1}{2}}\right)\rho ^{(\mathrm{eq})}(\omega ,𝐩),`$ (2.17)
$`\mathrm{\Sigma }_F^{(\mathrm{eq})}(\omega ,𝐩)`$ $`=`$ $`i\left(n_\mathrm{B}(\omega )+{\displaystyle \frac{1}{2}}\right)\mathrm{\Sigma }_\rho ^{(\mathrm{eq})}(\omega ,𝐩)`$ (2.18)
with $`n_\mathrm{B}(\omega )=(e^{\beta \omega }1)^1`$. While the spectral function $`\rho ^{(\mathrm{eq})}`$ encodes the information about the spectrum, one observes that the symmetric function $`F^{(\mathrm{eq})}`$ encodes the statistical aspects in terms of the particle distribution function $`n_\mathrm{B}`$. For a vanishing $`\omega `$-dependence the function $`\mathrm{\Gamma }^{(\mathrm{eq})}(\omega ,𝐩)\mathrm{\Sigma }_\rho ^{(\mathrm{eq})}(\omega ,𝐩)/2\omega `$ plays the role of a decay rate for one-particle excited states with momentum $`𝐩`$.

## 3 Initial conditions

We consider a spatially homogeneous system such that, after spatial Fourier transformation, the two-point functions depend only on $`𝐩`$ and we consider the modes $`\rho (x^0,y^0;𝐩)`$. As indicated above, the initial conditions for $`\rho (x^0,y^0;𝐩)`$ are fixed by the equal-time relations (2.12). As initial conditions for the symmetric two-point function $`F(x^0,y^0;𝐩)`$ we consider a “tsunami”–like situation  in which a spatially homogeneous collection of particles move with approximately the same momentum peaked around $`𝐩_{\mathrm{ts}}`$ and $`𝐩_{\mathrm{ts}}`$. Explicitly, we take

$`F_0(x^0,y^0;𝐩)|_{x^0=y^0=0}`$ $`=`$ $`{\displaystyle \frac{1}{\omega _𝐩}}\left[n_0(𝐩)+{\displaystyle \frac{1}{2}}\right],`$
$`_{x^0}_{y^0}F_0(x^0,y^0;𝐩)|_{x^0=y^0=0}`$ $`=`$ $`\omega _𝐩\left[n_0(𝐩)+{\displaystyle \frac{1}{2}}\right],`$ (3.1)
$`_{x^0}F_0(x^0,y^0;𝐩)|_{x^0=y^0=0}`$ $`=`$ $`0,`$
with $`\omega _𝐩=\sqrt{𝐩^2+m^2}`$ and the initial particle number

$$n_0(𝐩)=𝒩\mathrm{exp}\left(\frac{1}{2\sigma ^2}(|𝐩||𝐩_{\mathrm{ts}}|)^2\right).$$
(3.2)
Here $`\sigma `$ controls the width of the initial distribution and $`𝒩`$ is a normalization constant. The initial condition is clearly far from thermal equilibrium and reminiscent of two colliding high energy wave packets .

To study the nonequilibrium evolution of the spectral function $`\rho (x^0,y^0;𝐩)`$ we solve the coupled Eqs. (2.102.11) with (2.132.15) numerically along the lines of ref. . We use a standard discretization on a lattice in space and time for a finite spatial volume and spatially periodic boundary conditions. To remove finite size effects we increase the volume until the results become stable. Because the evaluation of the memory integrals is rather time- and memory-consuming the numerical results presented below are for a $`(1+1)`$–dimensional system. The lattice acts as a cutoff and regulates the ultraviolet divergences, which have to be renormalized. In $`1+1`$ dimensions this is straightforward since only the one-loop contribution is logarithmically divergent. We absorb this in a bare mass parameter $`\mu `$ with the replacement $`m^2\mu ^2=m^2\delta m^2`$ in Eq. (2.15). The counterterm $`\delta m^2`$ cancels the divergent vacuum contribution coming from the one-loop graph. The finite part of $`\delta m^2`$ is fixed by requiring that the renormalized one-loop mass parameter in vacuum ($`n_0(𝐩)0`$) equals $`m`$ and we express all dimensionful scales in units of $`m`$. For the plots we have used a space lattice with spacing $`ma=0.3`$ and a time lattice $`a_0/a=0.25`$. The system size is $`mL=24`$.

## 4 Wigner transformation

Derivations of transport or kinetic equations are typically performed by considering correlation functions in the so-called Wigner representation. Two-point functions are written in terms of the center-of-mass coordinate $`X=(x+y)/2`$ and the relative coordinate $`xy`$, with respect to which a Fourier transformation to momentum space is performed (Wigner transformation). Derivations based on a gradient expansion in $`_X`$ assume a slow variation of the two-point functions with $`X`$. In particular, the classical Boltzmann equation can be obtained from (2.102.15) using a first-order gradient expansion and an on-shell form of the spectral function with zero width . The full equations (2.102.15) can therefore be considered as a “quantum Boltzmann equation” including off-shell effects and resumming an infinite order of derivatives.

To analyze the evolution of the spectral function (not to solve the dynamics) we perform a Wigner transformation and write

$$i\rho (X^0;\omega ,𝐩)=_{2X^0}^{2X^0}𝑑te^{i\omega t}\rho (X^0+t/2,X^0t/2;𝐩).$$
(4.1)
The $`i`$ is introduced such that $`\rho (X^0;\omega ,𝐩)`$ is real. Because the spectral function is antisymmetric, $`\rho (X^0;\omega ,𝐩)=\rho (X^0;\omega ,𝐩)`$, we will present the positive-frequency part only. Furthermore, as a consequence of the equal-time relation (2.12) the Wigner transform obeys the sum rule $`I=𝑑\omega /(2\pi )\omega \rho (X^0;\omega ,𝐩)=1`$. In our simulations an evaluation of the sum rule typically results in $`1.000<I<1.003`$, indicating the accuracy in the numerical evaluation of the Wigner transform.

Since we consider an initial-value problem with $`x^0,y^00`$, the time integral over $`t=x^0y^0`$ is bounded by $`\pm 2X^0`$. This leads to a so-called “finite-time effect” in the Wigner transformed quantities. To understand the qualitative aspects of a finite-time interval it is instructive to consider the free spectral function $`\rho _0(x^0,y^0;𝐩)=\omega _𝐩^1\mathrm{sin}[\omega _𝐩(x^0y^0)]`$ with $`\omega _𝐩=\sqrt{𝐩^2+m^2}`$. After performing a Wigner transformation, we find

$$\rho _0(X^0;\omega ,𝐩)=\frac{\mathrm{sin}[(\omega \omega _𝐩)2X^0]}{\omega _𝐩(\omega \omega _𝐩)}\frac{\mathrm{sin}[(\omega +\omega _𝐩)2X^0]}{\omega _𝐩(\omega +\omega _𝐩)}.$$
(4.2)
For finite $`X^0`$ this spectral function shows a rapidly oscillating behavior, while its envelope is peaked at $`\omega =\pm \omega _𝐩`$. In the limit $`X^0\mathrm{}`$ it reduces to $`\rho _0(\omega ,𝐩)=2\pi \text{sgn}(\omega )\delta (\omega ^2\omega _𝐩^2)`$.

For an interacting system, however, effective damping leads to a finite correlation time $`\tau `$, such that finite-time effects are expected to vanish when $`X^0\tau `$. This can be investigated analytically by assuming that the spectral function is strictly exponentially damped and can be written as $`\rho (x^0,y^0;𝐩)=e^{\gamma _𝐩|x^0y^0|}E_𝐩^1\mathrm{sin}[E_𝐩(x^0y^0)]`$. The effective mass $`E_𝐩`$ and rate $`\gamma _𝐩`$ are allowed to depend on the time $`X^0`$. The corresponding Wigner transform reads $`\rho (X^0;\omega ,𝐩)=\rho _{\mathrm{BW}}(X^0;\omega ,𝐩)+\delta \rho (X^0;\omega ,𝐩)`$, where $`\rho _{\mathrm{BW}}`$ denotes the Breit-Wigner function

$$\rho _{\mathrm{BW}}(X^0;\omega ,𝐩)=\frac{2\omega \mathrm{\Gamma }_𝐩(X^0)}{[\omega ^2E_𝐩^2(X^0)]^2+\omega ^2\mathrm{\Gamma }_𝐩^2(X^0)}$$
(4.3)
with a width $`\mathrm{\Gamma }_𝐩(X^0)=2\gamma _𝐩(X^0)`$. The additional contribution $`\delta \rho (X^0;\omega ,𝐩)`$ vanishes exponentially as $`\mathrm{exp}(\mathrm{\Gamma }_𝐩X^0)`$, such that the finite-time effect is indeed absent for large enough (but finite) $`X^0`$. We want to emphasize that this finite-time effect is a general consequence of an initial-time problem formulated in terms of Wigner transformed quantities and should be taken into account when used to describe the early-time behaviour ($`\mathrm{\Gamma }_𝐩X^01`$) in an experimental setup.

Finally, we note that the positivity condition, $`\text{sgn}(\omega )\rho (X^0;\omega ,𝐩)0`$, can only be shown to hold in the special case that the initial density matrix commutes with the full Hamiltonian, such as in thermal equilibrium. In this case the system is of course stationary and independent of $`X^0`$. As a consequence the interpretation of the nonequilibrium spectral function as the density of states should be taken with care.

## 5 Nonequilibrium evolution

In Fig. 2 we show the time evolution of the spectral function $`\rho (X^0+t/2,X^0t/2;p)`$ as a function of the relative time $`t`$ for two different momentum modes at given $`X^0`$, as obtained from the numerical solution of Eqs. (2.102.15). From the logarithmic plot one observes an effective damping of the oscillating spectral function while the oscillations never damp out completely. The decrease of the maximum amplitude very quickly approaches an approximately exponential behavior characteristic for the decay of one-particle excited states with given momentum. One observes that the zero momentum mode is stronger damped than the higher momentum mode.555An analysis of damping times as a function of the coupling $`\lambda `$ can be found in Ref. .

The effective damping should correspond to a nonzero width of the Wigner transformed spectral function. In Fig. 3 we display the Wigner transforms $`\rho (X^0;\omega ,p)`$ for the zero momentum mode (left peak) and the higher momentum mode (right peak) for the same $`X^0`$-time as in Fig. 2. One clearly observes that the interacting theory has a continuous spectrum described by a peaked spectral function with a nonzero width. The inset shows a blow-up of the zero mode around the three-particle threshold $`3E_p/m=4.38`$. The expected bump in the spectral function is small but visible. We stress that this bump in the spectral function is kinematically forbidden for the on-shell limit and arises purely from off-shell decay. In Fig. 3 we also present fits to a Breit-Wigner spectral function. While the position of the peak can be fitted easily, the overall shape and width are only qualitatively captured. In particular, the slope of $`\rho (X^0;\omega ,p)`$ for small $`\omega `$ is quantitatively different. We also see that the Breit-Wigner fits give a narrower spectral function (smaller width) and therefore would predict a slower exponential relaxation in real time.

We now turn to the dependence of $`\rho (X^0;\omega ,p)`$ on $`X^0`$. In Fig. 4 we show the spectral function at zero momentum $`\rho (X^0;\omega ,p=0)`$ for three different values of $`X^0`$. We see that the peak slowly shifts to smaller values. In the inset to Fig. 4 this behavior is quantified by plotting the $`X^0`$-dependence of $`E_p(X^0)`$ obtained from a fit to the Breit-Wigner function (4.3). For comparison, the solid line in the same figure shows the time-dependent one-particle energy $`ϵ_p(X^0)`$ as defined via equal-time two-point functions in Eq. (5.2) below. We see a good agreement of both definitions of the one-particle energy for large enough times. We come back to this below.

A qualitative way to characterize the $`X^0`$-dependence of the shape of the spectral function is to assume that it is approximately described by a Breit-Wigner function. The time dependence of $`\mathrm{\Gamma }_p(X^0)`$ can then be followed by performing fits for different values of $`X^0`$. The resulting $`\mathrm{\Gamma }_p(X^0)/2`$ is shown in Fig. 5 for two momentum modes. After a variation for early times, the $`X^0`$-dependence becomes weaker for larger times. As indicated before (see Fig. 3), the Breit-Wigner fits result in a smaller width than what is expected by comparing with the actual spectral function. A simple way to obtain a quantitative description is by extracting a relaxation rate $`\gamma _p`$ directly in real time, i.e. from a fit of the spectral function $`\rho (X^0+t/2,X^0t/2;p)`$ to an exponential form $`\mathrm{exp}[\gamma _p(X^0)t]`$ (see Fig 2). The resulting rates are also presented in Fig. 5. We observe that the relaxation rate for the zero mode increases first, corresponding to a broadening of the Wigner transformed spectral function, before it levels off at later times. Despite the quantitative differences we see a similar qualitative $`X^0`$-dependence from the Breit-Wigner fit for large enough times.

Particle distribution function. The relatively weak $`X^0`$-dependence of the kinetic aspects of the evolution encoded in the spectral function has to be confronted with the time dependence of the statistical aspects encoded in $`F`$. The particle distribution is expected to exhibit strong qualitative changes, since we start with an initial condition (3.2) far away from thermal equilibrium and the current approximation is known to describe the physics of thermalization at large enough times .

The evolution of the statistical information can be extracted from the equal-time correlation function $`F(t,t;p)`$. If the theory exhibits a quasi-particle structure it is possible to define a quasi-particle distribution function from a combination of equal-time two-point functions . For this purpose we introduce the two-point function $`K(x,y)=\frac{1}{2}[\pi (x),\pi (y)]_+`$, where $`\pi (x)=_{x^0}\varphi (x)`$ is the canonical momentum. In the spatially homogeneous situation we consider here, $`K(x^0,y^0;p)=_{x^0}_{y^0}F(x^0,y^0;p)`$. We define the distribution function $`n(\epsilon _p(X^0),X^0)`$ at time $`X^0`$ as

$$n(\epsilon _p,X^0)+\frac{1}{2}=\left[K(X^0,X^0;p)F(X^0,X^0;p)\right]^{1/2}.$$
(5.1)
For the free theory the definition of $`n(\epsilon _p,X^0)`$ corresponds to the standard one-particle distribution. The distribution function depends on a time-dependent one-particle energy, defined as

$$\epsilon _p(X^0)=\left[K(X^0,X^0;p)/F(X^0,X^0;p)\right]^{1/2}.$$
(5.2)
The applicability of this definition to describe one-particle energies for the present nonequilibrium situation finds an impressive demonstration by comparing it with the position of the peak in the spectral function as discussed above (cf. Fig. 4).

In Fig. 6 we show the results for the effective particle number $`n(\epsilon _p,X^0)`$ as a function of $`\epsilon _p`$ at four values of $`X^0`$. The initially high occupation number in a small momentum range rearranges as time goes on. The original “tsunami” vanishes and the infrared modes become more excited. The distribution is not yet thermal. Thermalization takes place on longer time scales which is shown for this model in Ref. . The strong qualitative changes of the statistical aspects encoded in the particle distribution function have to be confronted with the comparably moderate changes of the spectral function shown in Fig. 4 and 5.

## 6 Equal-time spectral function and moments

In order to make a connection with methods employed before, we discuss some equal-time properties of the spectral function. As has been discussed earlier, the spectral function is the commutator of fields and obeys elementary equal-time properties (2.12). There is a further hierarchy of equal-time relations that follow from the exact evolution equation (2.10). The identity

$$_{x^0}^2\rho (x^0,y^0;𝐩)|_{x^0=y^0}=0$$
(6.1)
can also be traced back directly to the operator equation of motion. The following three equal-time derivatives are less trivial and read

$`_{x^0}^3\rho (x^0,y^0;𝐩)|_{x^0=y^0}`$ $`=`$ $`\left[𝐩^2+M^2(x)\right],`$
$`_{x^0}^4\rho (x^0,y^0;𝐩)|_{x^0=y^0}`$ $`=`$ $`2_{x^0}M^2(x),`$ (6.2)
$`_{x^0}^5\rho (x^0,y^0;𝐩)|_{x^0=y^0}`$ $`=`$ $`\left[𝐩^2+M^2(x)\right]^23_{x^0}^2M^2(x)`$
$`_{x^0}\mathrm{\Sigma }_\rho (x^0,y^0;𝐩)|_{x^0=y^0}.`$
These expressions are exact. One observes that the self energy $`\mathrm{\Sigma }_\rho `$, necessary to describe scattering, enters for the first time only in the fifth derivative. This suggests that in formalisms which are strictly based on equal-time correlation functions , the kinetic aspects of scattering encoded in the spectral function may only be captured in very sophisticated truncations.

The above relations can be immediately recast in terms of moments of the spectral function in the Wigner presentation. We define the $`n`$’th moment as

$$A_n(X^0;𝐩)=_{\mathrm{}}^{\mathrm{}}\frac{d\omega }{2\pi }\omega ^n\rho (X^0;\omega ,𝐩).$$
(6.3)
Due to the antisymmetry of $`\rho (X^0;\omega ,𝐩)`$ all even moments $`A_{2n}(X^0;𝐩)`$ vanish. It is straightforward to relate the odd moments directly to the equal-time derivatives given above, using $`_{x^0}=\frac{1}{2}_{X^0}+_t`$. The first moment is just the sum rule, $`A_1(X^0;𝐩)=1`$ and Eq. (6.1) provides a consistency check, $`_{x^0}^2\rho (x^0,y^0;𝐩)|_{x^0=y^0}=_{X^0}A_1(X^0;𝐩)=0`$. For the higher derivatives one finds

$`_{x^0}^3\rho (x^0,y^0;𝐩)|_{x^0=y^0}`$ $`=`$ $`A_3(X^0;𝐩),`$
$`_{x^0}^4\rho (x^0,y^0;𝐩)|_{x^0=y^0}`$ $`=`$ $`2_{X^0}A_3(X^0;𝐩),`$ (6.4)
$`_{x^0}^5\rho (x^0,y^0;𝐩)|_{x^0=y^0}`$ $`=`$ $`{\displaystyle \frac{5}{2}}_{X^0}^2A_3(X^0;𝐩)+A_5(X^0;𝐩).`$
The combination of (6.4) and (6.2) gives information on the time dependence of the lowest moments. Therefore, also in approximation schemes that are based on an expansion in moments, the lowest quantity that shows a direct dependence on the self energy $`\mathrm{\Sigma }_\rho `$ appears only at relatively high (fifth) order. We conclude that both in equal-time formalisms and in approximations based on moments a proper inclusion of scattering effects might be difficult to achieve.

## 7 Conclusion

The presence of a nonequilibrium environment affects the quasiparticle structure of the theory, encoded in the spectral function. As a relevant physical application one may think of high energy heavy-ion collisions where spectral properties of the vector mesons are affected by the presence of nuclear matter.

In this paper we have concentrated on a self-consistent determination of the spectral function out of equilibrium, using a scalar field theory in $`1+1`$ dimensions as our (toy) model. We have calculated the nonequilibrium time evolution of the spectral function by solving the equations of motion obtained from the three-loop $`2PI`$ effective action in real time. To compare the dynamics of the spectral function and the statistical aspects of the theory we also compute the evolution of an effective particle number distribution. Scattering and off-shell contributions are included in the dynamical evolution, in contrast to well-studied Hartree and leading-order large $`N`$ approximations. Even for moderate couplings and initial conditions far away from equilibrium, we observe a rather weak dependence of the Wigner transformed spectral function on $`X^0`$. This fact provides a necessary condition for a successful gradient expansion typically employed in the derivation of kinetic equations. After variations at early times we observe that the spectral function approaches a rather stable shape characterized by a nonvanishing width. A consistent “quantum Boltzmann equation” certainly has to relax the zero-width assumption employed in the derivation of the standard Boltzmann equation. We have argued that it might be notoriously difficult to derive the kinetic aspects of scattering encoded in the spectral function from equal-time formalisms or approximation schemes based on an expansion in moments. We emphasize that the employed framework based on the $`2PI`$ effective action allows for a consistent description of both the early and the late time behavior. Since a change in the position of the peak and width of the spectral function can be described consistently in this framework, it would be interesting to extend the study presented here to more realistic models for heavy-ion collisions.

Finally the analysis performed in this paper was largely motivated by the physics of relativistic heavy-ion collisions. We would like to stress, however, that the approach using the $`2PI`$ effective action is general and that possible applications of the formalism range from the early universe to condensed matter systems.

Acknowledgements

J.B. thanks Jürgen Cox for collaboration on related subjects. We would like to thank Tomislav Prokopec for many discussions. This work was supported by the TMR network Finite Temperature Phase Transitions in Particle Physics, EU contract no. FMRX-CT97-0122.