ence phenomena in excitation energy transfer has been developed long ago in the framework of the reduced den- sity matrix (RDM) [21, 22]. Equations of motion result- ing from this scheme are characterized by the presence of time retarded terms responsible for energy relaxation and decoherence processes. Equations of this type will be denoted as time non-local in this paper. Later, an alter- native approach to the derivation of the RDM equations of motion has emerged which yields time local equation of motion [23, 24]. Both theories express the relaxation term in form of an infinite series in terms of the system- bath interaction Hamiltonian, but differ in time ordering prescriptions for the cumulant expansion of the evolu- tion operator. The time local equations correspond to so-called partial time ordering prescription of the cum- mulant expansion, while the time non-local equations re- sult from so-called chronological time ordering [25, 26]. Although the two schemes yield formally different equa- tions of motion for the RDM, they are in fact equivalent as long as the complete summation of the corresponding infinite series is performed. When the infinite series are truncated at a finite order, the two theories yield equa- tions that predict different RDM dynamics. This is a result of different statistical assumptions about the bath that are implicitly made in the two cases [25, 26]. In all orders of expansion, so-called Markov approximation can be used to transform the time non-local equation of motion into a certain time-local form. This has to be re- garded as an additional approximation which simplifies the numerical treatment of the time non-local equations. Interestingly, in the second order the time-local equations and the time non-local equations with Markov approxi- mation have exactly the same form. Until recently, most experiments were not sensitive to coherence between electronic levels. This allowed a host of further approximations to simplify equations of mo- tion. Most notably, the secular approximation , which amounts to decoupling RDM elements oscillating on dif- ferent frequencies from each other, has limited the energy transfer phenomena to separate dynamics of population transfer and coherence dephasing [27]. Even on very short times scale, experiments aimed at studying pop- ulation dynamics (pump probe) are not sensitive enough to coherences between electronic levels to require non- secular theory, although it was suggested that measured relaxation time can be distorted by non-secular effects [28]. Consequently, most of the theory developed for eval- uation of experiments has been aimed at improving calcu- lation of the relaxation rates [29–31]. With experiments now uncovering new details about the role of electronic coherence, theoretical methods beyond rate equations for probabilities which are both accurate and numerically tractable are required. Although schemes for construct- ing equations of motion for the RDM beyond second or- der, based on co-propagation of the RDM with auxiliary operators, seem feasible and promising [32, 33], second order theories might still be the only option for treatment of extended molecular systems. It was suggested previ- where | g$_{m}$ 〉 , | e$_{m}$ 〉 denote electronic ground and excited states, and ε ( m ) g , ε ( m ) e represent electronic energies of The investigated molecular system is an aggregate composed of N monomers embedded in protein environ- ment. Let us first consider a monomeric molecule (a chromophore) embedded in the environment, but insu- lated from interaction with its neighboring monomers. The monomer Hamiltonian has a form The paper is organized as follows. The next section introduces Hamiltonian description of an aggregate of small molecules embedded in a protein or solid state environment. In Section III we describe the details of second order theory of system–bath interaction and we derive four different equations of motion for the RDM de- scribing electronic states of a molecular aggregate. Two- dimensional coherent spectroscopy, and non-linear spec- troscopy in general are introduced in Section IV. In Sec- tion V we present and discuss numerical results compar- ing different theories of relaxation on calculations of co- herence life time and 2D spectra. In this paper we study the following four different sec- ond order theories: (a) full time non-local (full TNL) equation of motion resulting from the Nakajima-Zwanzig identity or equivalently from the chronological ordering prescription in the cummulant expansion, (b) the full time local (full TL) equation of motion resulting from the partial ordering prescription in the cummulant ex- pansion, or equivalently from Markov approximation ap- plied to TNL equation, (c) time non-local equation with secular approximation (secular TNL), and (d) time local equation with secular approximation (secular TL). We discuss the applicability of these equations to the descrip- tion of the energy relaxation and decoherence dynamics in small systems of molecular excitons with the empha- sis of on recent 2D spectroscopic experiments and the dynamics of coherence between electronic excited states. Note that, in this paper, full refers to the equations where no secular approximation has been applied. These equa- tions are still of second order of perturbation theory. ously that second order perturbation theory with respect to system-bath coupling provides a suitable framework for development of such methods [34]. This notion is also supported by the fact that in the special case of so-called spin-boson model, second order time-local equation of motion already represents an exact equation of motion for the RDM [35]. 2 II. MODEL HAMILTONIAN H m = ( ε ( m ) g + T ( P$_{m}$ ) + V ( m ) g ( Q$_{m}$ ) | g$_{m}$ 〉〈 g$_{m}$ | + ( ε ( m ) e + T ( P$_{m}$ ) + V ( m ) e ( Q$_{m}$ ) | e$_{m}$ 〉〈 e$_{m}$ | , (1)