# Dynamics of exciton formation and relaxation in photoexcited semiconductors

###### Abstract

We investigate the dynamics of the exciton formation and relaxation on a picosecond time scale following a pulsed photoexcitation of a semiconductor. The study is conducted in the framework of the density matrix theory complemented with the dynamics controlled truncation scheme. We truncate the phonon branch of the resulting hierarchy of equations and propose the form of coupling among single–phonon–assisted and higher–order phonon–assisted density matrices so as to ensure the energy and particle–number conservation in a closed system. Time scales relevant for the exciton formation and relaxation processes are determined from numerical investigations performed on a one–dimensional model for the values of model parameters representative of a typical organic and inorganic semiconductor. The exciton dynamics is examined for different values of central frequency of the exciting field, temperature, and microscopic model parameters, such as the strengths of carrier–carrier and carrier–phonon couplings. We find that for typical organic semiconductor parameters, formation of bound excitons occurs on a several–hundred–femtosecond time scale, while their subsequent relaxation and equilibration takes at least several picoseconds. These time scales are consistent with recent experimental studies of the exciton formation and relaxation in conjugated polymer–based materials.

###### pacs:

71.35.-y,71.10.-w## I Introduction

The continual and ever–increasing demand for economic and efficient ways of utilizing solar energy drives a huge part of current research activities. In particular, organic solar cells have developed rapidly in the past decade and have become promising candidates for economically viable large–scale power generation due to their flexibility, cost effectiveness, relatively simple fabrication techniques, and mass production. Tress (2014); Deibel and Dyakonov (2010) Processes upon which the operation of solar cells is based are the light absorption in a semiconducting material and the subsequent conversion of photons into mobile charge carriers that produce an electric current. Brédas et al. (2009); Clarke and Durrant (2010) An optical excitation of a semiconductor creates an exciton, i.e., an electron–hole pair in which Coulomb attraction between oppositely charged electron and hole prevents their separation. In a conventional inorganic semiconductor, relatively weak Coulomb interaction (primarily due to large dielectric constant) results in the exciton binding energy of the order of 10 meV. Knox (1963); Singh and Ruda (2006); Axt and Kuhn (2004) Thus, thermal excitations are likely to split the exciton in an electron and a hole. On the other hand, in a typical organic semiconductor, the attraction between an electron and a hole is much stronger (mainly due to low dielectric constant), the exciton binding energy being of the order of or larger than 500 meV. Bässler and Köhler (2012); Brédas et al. (2009) Therefore, while optical absorption in an inorganic semiconductor results in almost immediate generation of free charges, in an organic semiconductor it leads to formation of tightly bound electron–hole pairs, which should be separated in order to generate current. Brédas et al. (2009); Tress (2014); Clarke and Durrant (2010) This last conclusion has an enormous impact on the design and geometry of organic photovoltaic devices.

Photoexcitation of a semiconductor creates electron–hole pairs in a highly nonequilibrium state. Apart from the Coulomb interaction, which primarily induces correlations, the carrier–phonon interaction is also vital for a thorough understanding of nonequilibrium processes taking place in photoexcited semiconductors. Theoretical approaches for treating these processes are most often based on the density matrix theory Kuhn (1998); Rossi and Kuhn (2002) or the nonequilibrium Green’s functions formalism. Haug and Jauho (1996) Density matrix theory has become the preferred technique in the treatment of experiments with ultrashort pulses, since it deals with quantities that depend on one time argument and are directly related to observables.

Previous theoretical studies of the exciton formation process after an ultrafast optical excitation of a semiconductor were typically focused on inorganic semiconductors. Early studies were conducted in the framework of the semiclassical Boltzmann approach. Selbmann et al. (1996); Gulia et al. (1997) The fully microscopic and quantum theory for the interacting system of electrons, holes, photons, and phonons, capable of treating a wide variety of optical and excitonic effects after an ultrafast optical excitation of a semiconductor, was elaborated in Refs. Kira and Koch, 2006, 2012; Kira et al., 2001; Hoyer et al., 2002, 2003. On the other hand, the exciton formation from an initial state of two opposite charges in organic semiconductors was typically modeled by simulating the time evolution of empirical Hamiltonians applied to small systems, where the effects of the lattice are not included or are treated classically. Kobrak and Bittner (2000); Tandon et al. (2003)

The main aim of this work was to investigate the dynamics of exciton formation on short (up to several ps) time scale. This time scale is of particular relevance for the operation of organic solar cells, since it has been well established that the exciton separation at the interface of donor and acceptor materials occurs on a subpicosecond time scale. Sariciftci et al. (1992); Hwang et al. (2008) However, the details of the exciton formation and separation process and the factors that determine its efficiency are still not well understood. In recent years, significant insights have been obtained from subpicosecond time–resolved experiments performed both on neat materials Banerji et al. (2010, 2011) and blends. Cowan et al. (2012); Grancini et al. (2013); Bakulin et al. (2012); Jailaubekov et al. (2013); Gélinas et al. (2014) The results of all these experiments highlight the importance of nonequilibrium nature of excitons formed after photoexcitation.

In our study, we employ the Hamiltonian which includes all relevant physical effects in the system: electronic coupling which leads to band formation, electron–hole interaction which causes exciton formation, electron–phonon interaction that leads to relaxation, and the interaction with external electromagnetic field. We do not, however, include the effects of stimulated emission which lead to radiative recombination of excitons, since we are interested in the exciton dynamics on a short time scale, where these effects are negligible. From the time evolution of relevant quantities, we identify the time scale of the processes of formation of free charges and bound excitons and their subsequent relaxation. Rather than focusing on the details of one particular material system, we have chosen a Hamiltonian whose parameters can be easily varied so that we can identify the influence of different physical effects on relevant time scales. The study is conducted in the framework of the density matrix formalism combined with the so–called dynamics controlled truncation (DCT) scheme, firstly developed in 1994 by Axt and Stahl. Axt and Stahl (1994a, b) This method is particularly suited for a system described by a pair–conserving Hamiltonian which is initially unexcited and was successfully applied to study the dynamics of exciton formation for near–band–gap excitations and low–excitation densities. Siantidis et al. (2001, 2002a, 2002b) Here, we truncate the phonon branch of the hierarchy so as to ensure that the resulting equations are compatible with the energy and particle–number conservation in a closed system. Furthermore, we propose the form of coupling between single–phonon–assisted and higher–order phonon–assisted electronic density matrices which is compatible with the energy conservation in a closed system.

The paper is organized as follows. In Sec. II the general form of the Hamiltonian, along with the equations which describe the exciton formation process, is presented. Section III is devoted to the results of our numerical investigations of the exciton formation process which are carried out on a one–dimensional model system. The discussion of our results in light of recent experimental investigations of ultrafast exciton dynamics is presented in Sec. IV, whereas concluding remarks are given in Sec. V.

## Ii Theoretical framework

We use the standard two–band semiconductor model which takes into account the interaction of carriers with the external electromagnetic field applied to the semiconductor, as well as carrier–carrier and carrier–phonon interactions. We will work in the electron–hole picture which is particularly suited for describing the effects which arise after the optical excitation of an initially unexcited semiconductor. Notation from Ref. Axt and Mukamel, 1998 will be used. The Hamiltonian has the form

(1) |

where describes interacting carriers

(2) |

(3) |

is the free–phonon Hamiltonian, describes the carrier–phonon interaction

(4) |

whereas the coupling to the optical field is given by

(5) |

Fermi operators () create (annihilate) an electron of energy in the single–particle state in the conduction band, while Fermi operators () create (annihilate) a hole of energy in the single–particle state in the valence band. Matrix elements of the Coulomb interaction potential are defined as

(6) |

where are single–particle eigenfunctions for an electron in the state and in the band . Bose operators () create (annihilate) a phonon in mode , while are carrier–phonon matrix elements. We neglect intraband contributions to the carrier–field interaction and retain only interband dipole matrix elements.

We note that the Hamiltonian of interacting carriers [Eq. (2)] includes the limiting cases of Wannier and Frenkel excitons. Namely, when single–particle eigenfunctions are of the Bloch form labeled by a wave vector , then under suitable approximations, described e.g., in Ref. Hanamura and Haug, 1977, we obtain the Hamiltonian describing the limiting case of Wannier excitons. On the other hand, if single–particle eigenfunctions are taken to be atomic states labeled by a position vector , then using approximations that exploit localization properties of this basis set the Hamiltonian appropriate for the limiting case of Frenkel excitons is obtained. Combescot and Pogosov (2008)

We study the dynamics of exciton formation in photoexcited semiconductors in the framework of the density matrix theory. Differential equations for dynamic variables are formed and, due to the many–body nature of the problem, an infinite hierarchy of differential equations is obtained. The main approximation is then the truncation of the hierarchy, which can be based upon different physical pictures. The Hamiltonian defined by Eqs. (1)(5) has the property that only the interaction with the optical field can change the number of pair excitations. The DCT scheme relies upon the aforementioned property of the Hamiltonian and classifies higher–order density matrices according to their leading order in the optical field. Axt and Stahl (1994a); Axt et al. (1996); Axt and Mukamel (1998) Namely, when the system is initially in the ground state represented by the vacuum of electron–hole pairs, the expectation value of the normal–ordered product of electron operators and , hole operators and and an arbitrary number of phonon operators and is at least of the order in the applied field. Therefore, higher–order density matrices are also of higher order in the optical field and only a finite number of electronic density matrices contributes to the optical response at any given order in the optical field. The DCT scheme truncates only the electronic branch of the hierarchy and can be used along with any strategy to deal with the phonon–assisted branch of the hierarchy. Axt and Kuhn (2004) We limit ourselves to the case of weak optical field and low carrier densities, in which it is justified to neglect biexcitonic effects and keep only contributions up to the second order in the optical field. In Refs. Siantidis et al., 2001; Axt and Mukamel, 1998 a reduced treatment of the phonon branch of the hierarchy, which can be combined with the DCT scheme for the electronic branch of the hierarchy, was presented. This treatment includes correlation expansion for phonon–assisted variables combined with the Markov approximation. As a result, phonon–assisted variables are eliminated from the formalism and only two types of electronic density matrices remain. These are the interband transition amplitude (excitonic amplitude)

(7) |

and the electron–hole pair density (excitonic population)

(8) |

In this study, we adopt a different strategy for dealing with the phonon–assisted density matrices. In order to facilitate the truncation of the phonon–assisted branch of the hierarchy, the following generating functions for the phonon–assisted electronic density matrices are defined: Axt and Mukamel (1998)

(9) |

(10) |

(11) |

where and are arbitrary sets of real parameters. As a consequence of the generating–function property, all phonon–assisted electronic density matrices can be obtained as derivatives of these functions taken at . The electron and hole populations and correlations and , as well as their phonon–assisted counterparts, do not have to be considered as independent variables in the formalism since they can be eliminated in favor of by identities (contraction identities) that are exact within the second–order treatment. Axt and Mukamel (1998); Axt et al. (1996) The differential equations for variables and are given in Appendix A.

The most general form of an electron–hole pair state is Hanamura and Haug (1977)

(12) |

where represents the state in which the conduction band is completely empty and the valence band is completely filled. The excitonic basis is defined by the eigenvalue problem which can be transformed into equations for amplitudes :

(13) |

The excitonic basis is orthonormal

(14) |

We perform all calculations in the excitonic basis and expand all density matrices in the excitonic basis, for example

(15) |

(16) |

and similarly for the corresponding phonon–assisted electronic density matrices; in the case of single–phonon assistance, the explicit definitions are

(17a) | |||||

(17b) |

The creation operator for the exciton in the state can be defined as

(18) |

The number of excitons in the state , after performing the decoupling (which is exact up to the second order in the optical field) , where stands for the correlated part of the electron–hole pair density, can be expressed as the sum

(19) |

where . The first term in Eq. (19) describes the so–called coherent excitons, whereas the second term describes the incoherent excitons. Namely, an optical excitation of a semiconductor firstly induces single–particle excitations in form of optical polarizations and carrier densities. Optical polarizations decay very fast due to various scattering mechanisms present. Kira and Koch (2012) Therefore, their squared moduli, which are usually referred to as coherent excitonic populations, Siantidis et al. (2001) do not provide information about the true excitonic populations, which are the consequence of Coulomb–induced correlations between electrons and holes and which typically exist in the system for a long time after the decay of optical polarizations. Axt and Kuhn (2004) In order to describe true excitons, which are atomlike complexes of electrons and holes bound by the Coulomb attraction, we have to consider two–particle correlations between them, and not single–particle quantities. Kira and Koch (2012) The last conclusion justifies identification of the term with the incoherent excitonic populations.

The dynamic equations for the relevant variables should be compatible with the energy conservation in a system without external fields. Our system, however, interacts with external optical field, but, since we consider a pulsed excitation, the energy of the system should be conserved after the field has vanished. The total energy of the system, i.e., the expectation value of the Hamiltonian defined in Eqs. (1)–(5), is expressed as

(20) |

where the carrier energy is

(21) |

the phonon energy is

(22) |

the carrier–phonon interaction energy is

(23) |

and the carrier–field interaction energy is

(24) |

In Eqs. (20)–(24) we have kept only contributions up to the second order in the external field and transferred to the excitonic basis. We also introduce excitonic dipole matrix elements

(25) |

as well as matrix elements of the carrier–phonon interaction in the excitonic basis which describe the coupling to the phonon mode :

(26) |

Within previous approaches to solving the hierarchy of equations obtained after performing the DCT scheme, single–phonon–assisted density matrices , which appear in Eq. (23), were not explicitly taken into account, but the respective differential equations were solved in the Markov and adiabatic approximations. However, it can be shown that the total energy under these approximations is not exactly conserved after the external field has vanished. In order to satisfy the energy conservation, we retain density matrices as independent dynamic variables in the formalism.

The dynamics should also conserve the particle number after the external field has vanished, since all the other terms in the Hamiltonian given by Eqs. (1)–(5) commute with the total particle–number operator. The number of electrons (and also the number of holes, since carriers are generated in pairs in this model), with accuracy up to the second order in the external field, is given as

(27) |

The equations for the purely electronic relevant variables and phonon distribution function are

(28) |

(29) |

(30) |

Even at this level, without specifying the form of equations for one–phonon–assisted electronic density matrices, using Eq. (29) with vanishing electric field it is easily shown that, in the absence of external fields, our dynamics conserves the total number of particles.

We will neglect hot–phonon effects and assume that in all the equations for , , and their phonon–assisted counterparts the phonon numbers assume their equilibrium values . We will, however, retain Eq. (30) in the formalism because it is necessary to prove the energy conservation.

In equations for phonon–assisted electronic density matrices we neglect the coupling to the light field, i.e., we neglect contributions arising from the combined action of the phonon coupling and the interaction with the light field (so–called cross terms). Kuhn and Rossi (1992); Axt and Mukamel (1998) The equations for the electronic density matrices with one–phonon assistance contain electronic density matrices with two–phonon assistance, from which we explicitly separate the factorized part and the correlated part, for example

(31) |

(32) |

We should bear in mind that the two–phonon–assisted electronic density matrices with two creation (annihilation) phonon operators, whose factorized part vanishes, should be considered on this level of truncation of the phonon branch. Schilp et al. (1994) Further comments on the factorization performed in Eq. (31) are given in Appendix B. The following equations for the electronic density matrices with single–phonon assistance are obtained:

(33) |

(34) |

(35) |

The correlated parts of two–phonon–assisted density matrices appearing in Eqs. (33) (), (34), and (35) can be obtained solving their respective differential equations, in which all three–phonon–assisted density matrices have been appropriately factorized and their correlated parts have been neglected, in the Markov and adiabatic approximations. This procedure closes the phonon branch of the hierarchy. However, the full solution to these equations, when combined with Eq. (33), is cumbersome to evaluate, so further approximations are usually employed. The most common one is the so–called random phase approximation, which neglects sums over correlated parts of one–phonon–assisted electronic density matrices (which are complex quantities) due to random phases at different arguments of these density matrices. Kuhn (1998) After performing all the discussed approximations, the last two summands in Eq. (33), which represent the rate at which changes due to the coupling to electronic density matrices with higher phonon assistance, reduce to

(36) |

where is given as

(37) |

(38) |

Details of the procedure employed to close the phonon branch of the hierarchy are given in Appendix B.

It was recognized that this form of the coupling to higher–order phonon–assisted electronic density matrices is at variance with the energy conservation. Schilp et al. (1995); Kuhn (1998); Rossi and Kuhn (2002) In this work, we will use the following form of the coupling to higher–order phonon–assisted density matrices:

(39) |

where is, as before, defined by Eqs. (37) and (38). This form of is compatible with the energy conservation, as long as excitonic matrix elements of the carrier–phonon interaction are purely real, which is the case relevant for our numerical investigation in Sec. III. Namely, as is shown in Appendix C, the rate at which the total energy changes after the pulse is equal to the rate at which the carrier–phonon interaction energy changes due to the coupling of the single–phonon–assisted electronic density matrices to density matrices with higher–order phonon assistance,

(40) |

It is then clear that, if all are real, the form of given in Eq. (39) does not violate the energy conservation. Furthermore, as describes the elementary process in which an exciton initially in the state is scattered to the state emitting the phonon from the mode , the reverse microscopic process, described by , is also possible, so in the differential equation for the quantity may appear. In Appendix C, we comment on the energy conservation in greater detail.

Similar strategy can be adopted to simplify the coupling to electronic density matrices with higher phonon assistance in (34) and (35), with the final result

(41) |

where

(42) |

and is defined in Eq. (38).

An alternative route to derive equations for the relevant variables is to rewrite the Hamiltonian given in Eq. (1) in terms of operators [see Eq. (18)], keeping only contributions whose expectation values are at most of the second order in the optical field. The result is

(43) |

The excitonic operators (up to the second order in the optical field) satisfy Bose commutation relations In this representation, Zhang et al. (1999) the excitons are treated as noninteracting bosons and the form of the exciton–phonon interaction is transparent, with exciton–phonon coupling constants defined in Eq. (26).

## Iii One–dimensional semiconductor model and numerical results

Numerical computations will be carried out on a two–band one–dimensional semiconductor model. We use a tight–binding model on a one–dimensional lattice with sites and lattice spacing to describe the semiconductor. Periodic boundary conditions are used. The Hamiltonian describing interacting carriers is given as

(44) |

It is assumed that the carrier transfer integrals are non–zero only among nearest–neighbor pairs of sites. The Coulomb interaction is taken in the lowest monopole–monopole approximation, Huhn and Stahl (1984) and the interaction potential is taken to be the Ohno potential

(45) |

is the on–site carrier–carrier interaction, while is the characteristic length given as , where is the static relative dielectric constant. This form of carrier–carrier interaction is an interpolation between the on–site Coulomb interaction and the ordinary Coulomb potential (in which the static relative dielectric constant is taken) when (see, e.g., the discussion on the effective electron–hole interaction in Ref. Knox, 1963). The interaction with phonons is taken to be of the Holstein form, where a charge carrier is locally and linearly coupled to a dispersionless optical mode

(46) |

where the free–phonon Hamiltonian is

(47) |

Despite the fact that the carrier–phonon interaction in real materials has a more complicated form, we choose for our numerical investigations its simplest possible form (Eq. (46)) capable of providing the energy relaxation of the electronic subsystem. The interaction with the electric field is

(48) |

As the system described is translationally symmetric, we can transfer to the momentum space and obtain the same Hamiltonian as described in Eqs. (1)–(5) with the following values of parameters:

(49a) | |||

(49b) | |||

(49c) | |||

(49d) |

The signs of the transfer integrals are . The constant energy , while is chosen so that the maximum of the valence band is the zero of the energy scale. is the Fourier transformation of the Ohno potential and it is computed numerically as

(50) |

The translational symmetry of our model enables us to solve efficiently the eigenvalue problem (13) which defines the excitonic basis. Instead of solving eigenvalue problem of a matrix, we can solve independent eigenvalue problems of matrices of dimension , thus obtaining excitonic eigenstates and their eigenenergies, which are counted by the center–of–mass wave vector and the band index . Thus, the general index of an excitonic state should be, in all practical calculations, replaced by combination . This has the following consequences on the matrix elements in the excitonic basis: dipole matrix elements reduce to

(51) |

whereas carrier–phonon interaction matrix elements reduce to

(52) |

Due to the translational symmetry of our model, only the dynamic variables for which the total created wave vector is equal to the total annihilated wave vector will have nontrivial values in the course of the system’s evolution. For example, from all density matrices only those with can have non–zero values.

Our objective is to analyze, in the framework of this relatively simple model, the characteristic time scales of exciton formation and relaxation in a photoexcited semiconductor, along with the impact that various model parameters have on these processes. Basic parameters in our model are transfer integrals (which determine bandwidths of the conduction and valence bands), electron–phonon coupling constants , the phonon energy , the dielectric constant , and the on–site Coulomb interaction . We will, throughout the computations, assume for simplicity that and .

As is well known, the main differences between a typical organic and inorganic semiconductor can be expressed in terms of bandwidths, dielectric constant and the carrier–phonon interaction strength. Namely, inorganic semiconductors are characterized by wide bands and high value of dielectric constant, whereas organic semiconductors have narrow bands and small value of dielectric constant. The carrier–phonon interaction is stronger in organic than in inorganic semiconductors. Having all these facts in mind, we propose two sets of model parameters which assume values typical of an organic and inorganic semiconductor. Values of our model parameters are adjusted to material parameters of GaAs for the inorganic case and pentacene for the organic case. Values of carrier–phonon coupling constants are chosen to correspond to typical values for mobility and/or typical values for the polaron binding energy.

Typical bandwidths in organic semiconductors are , Bässler and Köhler (2012) which corresponds to the transfer integral , whereas inorganic semiconductors usually exhibit bandwidths of several electronvolts Bässler and Köhler (2012) and we take in our calculations the value of the transfer integral . In both cases, the lattice constant was fixed to . The dielectric constant in a typical inorganic semiconductor is of the order of 10 and in the calculations we take the value of static dielectric constant of GaAs . For a representative value of the dielectric constant in organic semiconductors we take . Bässler and Köhler (2012); Clarke and Durrant (2010) The value of the on–site Coulomb interaction is chosen to give the correct order of magnitude for the exciton binding energy, which is calculated numerically. For the organic parameter set, we set meV, which gives the exciton binding energy around 320 meV, while for the inorganic parameter set meV and the corresponding exciton binding energy is roughly 10 meV.

The carrier–phonon coupling constants for the inorganic case are estimated from the mobility values. The mobility of carriers is estimated using the relation , where is the scattering time and is the effective mass of a carrier. For cosine bands considered in this work, in the vicinity of the band extremum. The scattering time is estimated from the expression for the carrier–phonon inelastic scattering rate based on the Fermi’s golden rule, which around the band extremum assumes the following form

(53) |

where . Therefore, the carrier–phonon coupling constant in terms of the carrier mobility reads as

(54) |

Using the value for the electron mobility in GaAs at , Levinshtein et al. (1996) we obtain .

We can also estimate the carrier–phonon coupling constants from the polaron binding energy. As an estimate of this quantity, we use the result of the second–order weak–coupling perturbation theory at in the vicinity of the point : Cheng and Silbey (2008)

(55) |

It is known that polaron binding energies in typical inorganic semiconductors are and we used this fact along with Eq. (55) to check our estimate for from the value of mobility; for , we obtain . The polaron binding energies in polyacenes lie in the range between and . Pope and Swenberg (1999) The value of in the set of model parameters representative of organic semiconductors was estimated from the polaron binding energy in pentacene, which is around . We obtain that . The values used for the organic/inorganic set of parameters are listed in Table 1.

Parameter | Inorganic | Organic |
---|---|---|

(meV) | 1519 Siantidis et al. (2001) | 2000 Cudazzo et al. (2015) |

(meV) | 500 | 125 |

12.9 Siantidis et al. (2001) | 3.0 Bässler and Köhler (2012) | |

(meV) | 25 | 40 |

(meV) | 36.4 Siantidis et al. (2001) | 10.0 Troisi and Orlandi (2006); Fujii et al. (2007) |

(meV) | 15 | 480 |

The form of the electric field is assumed to be a rectangular cosine pulse

(56) |

where is the central frequency of the field and is the Heaviside step function. Time is chosen large enough so that the pulse is so spectrally narrow that the notion of the central frequency makes sense. On the other hand, the pulse should be as short as possible, so that after its end we observe the intrinsic dynamics of our system, the one which is not accompanied by the carrier generation process, but merely shows how initially generated populations are redistributed among various states. Trying to reconcile the aforementioned requirements, we choose fs. The amplitude of the electric field and the interband dipole matrix element are chosen so that we stay in the low–density regime; particularly, we choose them so that the corresponding Rabi frequency assumes the value of meV, which is smaller than any energy scale in our problem and ensures that the excitation is weak.

In order to quantitatively study the process of exciton formation after a pulsed excitation of a semiconductor, we solved the system of quantum kinetic equations for electronic density matrices and their single–phonon–assisted counterparts [Eqs. (28), (29), (33), (34), and (35) supplemented with Eqs. (36) and (41)] using the fourth–order Runge–Kutta algorithm. The computations are performed for the temperature K and the central frequency of the pulse equal to the single–particle gap (). The exciton is considered bound (unbound) if its energy is smaller (larger) than the smallest single–particle energy difference . Cudazzo et al. (2015) The equation of the boundary line which separates bound from unbound pair states reads as

(57) |

An unbound exciton may be considered as (quasi)free electron and hole, so this way it is possible to distinguish between bound excitons and free carriers.

The pulsed excitation of a semiconductor leads, in the first step, to the generation of coherent electron–hole pairs that are described in our formalism by the coherent pair amplitudes . The decay of the coherent pair occupation

(58) |

is due to the scattering processes which initiate already during the generation of the pairs and gives a direct measure of the loss of coherence. Siantidis et al. (2001) At the same time, incoherent pair occupations start to grow, driven by the loss rate of coherent pair occupations. Axt and Mukamel (1998); Siantidis et al. (2001) In order to quantify the process of exciton formation, we will follow the time dependence of the total number of incoherent bound excitons

(59) |

This quantity represents the number of truly bound electron–hole pairs which exist even after the optical field has vanished and as such is the direct measure of the efficiency of the exciton formation process. We will, when useful, also consider the number of incoherent excitons in a particular band , . The quantities and will be normalized to the total number of excitons defined in Eq. (27).

### iii.1 Numerical results: Organic set of parameters

We start this section by an overview of properties of the excitonic spectrum, shown in in Fig. 1(a), which will be relevant for further discussions of the exciton formation process.

The lowest excitonic band is energetically well separated from the rest of the spectrum, the energy separation between the minima of the bands and being around , which is much larger than both the value of at room temperature and the phonon energy in our model (see Table 1). As a consequence, downward transitions that end at the lowest excitonic band start almost exclusively from the states on band and an exciton, which is at some instant in a state on the band, cannot be scattered to an unbound excitonic state.

We briefly comment on the size of the exciton for these values of model parameters. From the exciton wave function in space, we can obtain the exciton wave function in real space performing the Fourier transformation

(60) |

The exciton wave function in real space is a product of the plane wave which describes the motion of the center of mass with the wave vector and the wave function of the relative motion of an electron and a hole:

(61) |

The latter part is directly related to the exciton size. We calculated squared modulus of the wave function of the relative motion of a pair for states in various bands. The result is shown in Fig. 1(b). It is clearly seen that an electron and a hole are tightly bound in these states and their relative separations are of the order of lattice constant, which is the typical value for the exciton radius in organic semiconductors. We point out that this does not mean that an exciton is localized; due to the translational symmetry of our system, it is delocalized over the whole lattice, as described by the plane–wave factor in the total wave function of a pair. Moreover, we note that the system size is large enough for the results to be numerically accurate, as it is much larger than the typical size of the exciton in a bound state.

The impact that different parameters have on the exciton formation process is studied by changing one parameter, at the same time fixing the values of all the other parameters to the previously mentioned ones. We performed all the computations for a limited number of lowest excitonic bands, which crucially depends on the central frequency of the excitation. For the given excitation, we took into account all the bands whose minima lie below , where is a numerical constant.

We will firstly discuss the exciton formation process for different central frequencies of the exciting pulse. We have considered central frequencies in resonance with state, state, single–particle gap and the central frequency which is above the band gap. As can be noted from Fig. 2, raising the central frequency of the laser field leads to lower relative number of incoherent bound excitons. Namely, the higher is the central frequency, the higher (in energy) are the bands in which the initial coherent excitonic populations are created and the slower is the conversion of these coherent populations to incoherent populations in lower excitonic bands.

However, in the long–time limit, the relative number of incoherent bound excitons should not depend on the central frequency of the laser, but tend to the value predicted by the Maxwell–Boltzmann distribution, which is above 99%. Such a high value is due to the large energy separation between the lowest excitonic band and the rest of the spectrum. We can thus infer, based on Fig.