Figure 1: Illustration of the level structure of an excitonic system. The excited states | e$_{n}$ 〉 of N monomers with tran- sition frequency Ω (left part of the figure) split due to the resonance interaction into N one exciton states | ¯ u$_{n}$ 〉 (right). Absorption spectra of ensembles of non-interacting (left) and interaction (right) monomers. The system also exhibits higher excitons states (two-exciton states | ¯ U$_{n}$ 〉 are depicted here), with Ω being the mean transition frequency from the one- to two-exciton bands. A pictorial 2D spectrum with peaks re- sulting from transitions between the ground- and one-exciton states (red arrows) and one- and two-exciton states (blue ar- rows) is presented in the upper left corner of the figure. The transitions between the ground- and one-excitons states lead to positive contributions to the 2D spectrum (absorption and ground state bleach), while the transitions between the one- and two-exciton states result in a negative contribution (ex- cited state absorption). A. Projection operator technique and Nakajima-Zwanzig identity B. Convolutionless equation of motion We start with the density operator W , which describes the state of the system composed of the aggregate and its surroundings. This operator fulfills Liouville-von Neu- mann equation Here, [ A, B ]$_{−}$ ≡ AB − BA , and we introduced Liouville superoperator L by its action on an arbitrary operator A as L A = 1 $_{ℏ}$[ H, A ]$_{−}$ . Liouville superoperators L$_{S}$ , L$_{B}$ and L$_{S}$$_{−}$$_{B}$ that correspond to Hamiltonian operators of This superoperator is of course unknown, but it allows us to turn a time non-local differential equation (14) for- mally into a time-local one. We can write Before we introduce approximations into Eq. (14), we can make one more formal step, which enables us to over- come its time non-local character. Let us assume an evo- lution superoperator U $^{′}$( t ) , which satisfies where the subscript I denotes interaction picture with re- spect to the bath Hamiltonian H$_{B}$ . Eq. (14) is known as Nakajima-Zwanzig identity [27, 36]. Although one can- not practically solve Eq. (14), because it is as difficult as the original complete problem, one can use it as a starting point for approximations that will follow. where w is a pure bath operator satisfying the condition tr$_{B}$ { w } = 1 . For convenience we also define a comple- mentary projection operator Q ≡ 1 − P . For the pro- jected density matrix P W ( t ) one can write a formally exact equation of motion as where tr$_{B}$ represents trace over all bath DOF, is the quantity carrying just the relevant information [27]. To derive equation of motion for ρ ( t ) we can use projection operator approach as following. We define a projection operator P by its action on an arbitrary operator A as We are often not interested in the complete information carried by W ( t ) , but only in the information related to the electronic states of the aggregate. The reduced den- sity matrix (RDM) defined as the system, the bath and their interaction will be used later in this paper. A formal solution of Eq. (10) can be written using evolution superoperators such that 4 W ( t ) = U ( t − t$_{0}$ ) W ( t$_{0}$ ) . (11) ρ ( t ) = tr$_{B}$ { W ( t ) } , (12) P A = tr$_{B}$ { A } w, (13) ∂ ∂t P W ( t ) = − i PL$_{I}$ ( t ) P W ( t ) −PL$_{I}$ ( t ) exp { − i ˆ t t 0 dτ QL$_{I}$ ( τ ) Q } Q W ( t$_{0}$ ) − t − t 0 ˆ 0 dτ $^{′}$PL$_{I}$ ( t ) exp { − i ˆ τ ′ t 0 dτ $^{′′}$QL$_{I}$ ( τ $^{′′}$) Q } × QL$_{I}$ ( t − τ $^{′}$) P W ( t − τ $^{′}$) , (14) P W ( t ) = U $^{′}$( t − t$_{0}$ ) P W ( t$_{0}$ ) . (15) ∂ ∂t P W ( t ) = − i PL$_{I}$ ( t ) P W ( t ) ∂ ∂t W ( t ) = − i $_{ℏ}$[ H, W ( t )]$_{−}$ = − i L W ( t ) . (10)