Quasiparticle representation of coherent nonlinear optical signals of multi-excitons

Multi-exciton Green's functions and scattering matrices in many-fermion systems are calculated using a quasiparticle approach based on a generalized Bethe–Salpeter equation. The simulation protocol only requires numerical diagonalization of the single-exciton manifold. Using coboson algebra all many-body effects are recast in terms of two tetradic exciton–exciton interactions: direct Coulomb scattering and Pauli exchange. The tedious equations-of-motion derivations and calculations of multi-exciton manifolds are avoided. The approach is applied to calculate the third- and fifth-order signals generated by sequences of femtosecond optical pulses. Several coherent fifth order optical signals that directly probe three-exciton states and their projections on double and single-exciton states are predicted.


Introduction
Elementary excitations of many-fermion systems can be described as quasiparticles (QP) and are viewed as independent harmonic oscillators. The quasiparticle approach is widely used in the calculation of transport and optical properties of metals, semiconductors and strongly correlated quantum materials [1]. Examples include the Landau theory of metals, the Drude model, plasmons, magnons, etc [2][3][4]. Much effort has focused on the dispersion of single QP by utilizing a mean-field formalism such as the semiconductor Bloch equations [5]. Inherently, this level of theory only takes into account two-point, Coulomb-induced correlations between interacting particles. In smaller systems, such as quantum wells and dots, the interacting particles become 'strongly correlated' and interactions between them may no longer be treated as local scattering events [6]. Various nonlinear optical techniques can however create and manipulate a controlled number of QP [11,12] whose interactions are mapped into effective anharmonicities [13] and provide an interesting window into the many-body physics. These techniques are sensitive to multi-point correlations. Sequences of ultrafast optical pulses can be designed such that the signals vanish for non-interacting QP, thus providing a backgroundfree probe for anharmonic deviations from the independent particle picture. Four and six wave mixing techniques probe four-and six-point correlations respectively [6][7][8][9][10].
We propose a Green's function approach based on composite boson (coboson) algebra [22] for calculating such correlations and apply it onwards the study of two-and three-exciton states. The boson picture is attractive since it provides a convenient computational tool and a simple, physical, almost classical, picture of weak excitations that is exact in the linear response regime. The boson description gradually breaks down as the excitation level is increased and the response becomes nonlinear. Considerable effort has been devoted to developing various bosonization schemes which go beyond the linear response regime by transforming the original fermion degrees of freedom into weakly interacting bosons. These schemes attempt to treat higher levels of excitations by lumping all interactions into an effective 'direct' Coulomb scattering [14]. However, they cannot fully describe Pauli-driven exciton exchange by an effective Hermitian Hamiltonian. A proper many-body treatment requires various adjustments to the bosonization approaches. Several approaches treated the problem to first order in quasiparticle density by an effective scattering between two boson-excitons. For semiconductors, Hanamura and co-workers [15,16] employed the Usui transformation, which uses a non-Hermitian Hamiltonian, to calculate the nonlinear response. The boson commutation relations were q-deformed [17]. This theory was adopted in quantum optics to describe Dicke super-radiance. The Jordan-Wigner approach to Pauli exchange transforms the Pauli excitations in one-dimensional systems to non-or weakly interacting fermions [24]. This method has been successfully applied to J-aggregates [25]. Another example of such deformation was adopted by Mukamel and co-workers [18,19] for the Frenkel excitons and by Knoester et al for hardcore bosons [20]. Combescot and Betbeder-Matibet [21,22] developed a coboson formalism that works for higher quasiparticle densities. This requires an infinite expansion of both the Hamiltonian and commutation relations in coboson operators. The coboson algebra can be formulated in compact form by introducing two quantities: an exciton scattering operator V and an operator D representing the deviation from bosonicity. It is not necessary to derive the exact Hamiltonian in the coboson form in order to calculate various physical observables [23]. This paper focuses on developing a Green's function approach for the nonlinear optical response and identifying signatures of multiple excitons in multidimensional signals. Two protocols based either on sum-over-states (SOS) or QP are primarily used to calculate the nonlinear optical responses. Both start with interacting many-fermion systems whose ground state |g is given by a single Slater determinant which can be obtained at the Hartree-Fock or DFT level. The block-diagonal Hamiltonian is assumed to conserve the number of excitons (i.e. the off-diagonal blocks vanish). This Heitler-London approximation for Frenkel excitons is the only approximation on the Coulomb interaction that we will make. It is expected to hold for typical molecules and semiconductors. We consider the following Hamiltonian blocks: a single, non-degenerate ground state, N single electron-hole pair states, N (N −1) 2 two electron-hole pairs, etc. The n-pair blocks has C n N states, with N being the system size. In the SOS formalism the relevant blocks are diagonalized to yield multi-exciton states [18]. The nonlinear optical response is then recast in terms of multi-exciton Green's functions and transition dipole elements between the multi-exciton states. This protocol has several drawbacks since the multiexciton wavefunctions contain much more information than needed for typical observables. The numerical effort grows rapidly with the number of excitations. Moreover, massive cancellation of ∼N 2 scaling terms in various pathways to yield final ∼N signals complicates the calculation and analysis [26].
The QP approach mitigates these difficulties. A QP description has been developed earlier based on nonlinear exciton equations (NEE) [19,27], also known as dynamics controlled truncation approach [6]. In this picture, Pauli exclusion is exactly accounted for by the deformation of boson commutators. Products of uncorrelated electron/hole pairs are used as a basis set for multiple excitations. The nonlinear response originates from the scattering matrix ( (n) ), defined via the Bethe-Salpeter equation, which uses n quasiparticle boson Green's functions G (n,0) as a reference. Calculating the scattering matrix requires only simple algebraic manipulations (products and inversions of matrices); expensive matrix diagonalizations are not necessary. The order-by-order derivation of the NEEs is the bottleneck which so far had prevented its application to more than two interacting excitons. For example, in order to find (3) for the third-order response, one needs to solve C 1 N single-exciton equations coupled to C 2 N , C 3 N equations for double-and triple-excitons respectively. By treating the excitons as coboson QP we can calculate high-order (χ (5) and above) nonlinear optical responses. Our QP approach relies solely on the coboson algebra without resorting to equations of motion (NEE) or numerical multi-exciton diagonalization as in the SOS [18]. The required parameters for this algebra require knowledge of the two-exciton Hamiltonian block only. Explicit forms of the three-exciton and higher blocks are not needed. Several nonlinear, high-order optical techniques in the coboson representation will be derived by using loop diagrams. All necessary ingredients, such as boson Green's functions, coboson scattering matrices and transition dipole moments, will be directly obtained from the coboson algebra.
We shall further address some difficulties in recasting the SOS formalism in the coboson basis. In particular we shall focus on those arising from over-completeness of the two coboson manifold (N 2 size and with some unphysical states, compared with N (N − 1)/2 true twoexciton space). The multi-exciton Hamiltonian blocks become singular in this representation.
Enlarging the system size results in a singular part of the Hamiltonian, known as the self energy. Usually, in semiconductors, such effects are attributed to the system size reduction. The formalism may be extended to solve more complex models, e.g. by adding coupling to phonons. This will be done in appendix G by proper averaging of products of Green's functions.

The coboson representation of multi-excitons
The coboson algebra is the key element for our approach. It is instructive to derive its elements starting with the many-body Hamiltonian of a system of interacting fermions. These fermions are electrons at discrete position m1 created from the ground state c † m1 |g and holes at position m2 created via d † m2 |g . The respective creation/annihilation operators obey the fermionic anti-commutation rules: {c m1 , c † n2 } = {d m1 , d † n2 } = δ m1,n2 . We use the notation of [18] whereby subscripts in the form of letter+number define the electron and hole indices. For a collection of N two level chromophores those indices run over their positions and neglect spins. We adopt the tight-binding Hamiltonian in Heitler-London approximation to describe elementary excitations [19]. It is given by the sum of the single particle contribution H 0 , the Coulomb interaction H C and the dipole coupling with the radiation field H I : Here the matrices t describe the hopping between the sites. The tetradic matrices V (1) , V (2) describe electron/electron and hole/hole Coulomb repulsion, while W (1,2) represent electron-hole attraction. Equation (4) represents the interaction with right/left components of the optical field E ± (t) via transition dipole moments matrix µ. Hereafter matrices are denoted by bold letters.
To develop the QP picture for electronic excitations we shall switch from the above fermion picture to uncorrelated electron/hole pairs and eventually to the coboson QP picture. In order to transform the fermion Hamiltonian to the uncorrelated electron-hole representation, we first note that the matter Hamiltonian (equations (2) and (3)) conserves the number of excitations (electron-hole pairs). It is therefore possible to focus on one subspace of n-exciton states and construct a electron/hole pair Hamiltonian that coincides with the original fermion Hamiltonian in that space. Below we only consider the single |ψ x = B † n |g and two-exciton |ψ x x = B † m B † n |g spaces, but the method can be extended to higher spaces. The respective electron/hole pair creation operator is given by , with the subscript m standing for the pair m = m1, m2 in the site basis. When m1 = m2 it creates the Frenkel exciton and the charge transfer state otherwise. It is worth noting that physical excitons form a subspace of electron/hole pairs obtained by imposing Pauli constraints on the over-complete electron/hole pair basis. For the double excitations, the constraints assume m > n form, thus reducing the size of |ψ x x from N 2 to N (N − 1)/2. After some algebra, the resulting electron-hole pair Hamiltonian assumes the from [18,25] The above Hamiltonian is block-diagonal in the basis |ψ x , |ψ x x and conserves the number of excitations. To recast the Hamiltonian (equation (6)) in the coboson representation, we first diagonalize the single-particle block H 1 and obtain the single-exciton states |1 and energies E 1 : where we label the single-exciton states in increasing order of energy (i.e. |1 , |2 , etc). This allows one to formally introduce the coboson creation operators by projecting the uncorrelated electron/hole pair |m onto the single coboson states |1 . Back and forth transformations between the cobosons (excitons) and the electron-hole pairs are given by B † 1 |g = |1 . We shall adopt the following notation for subscripts: letter subscripts B † ν (ν = k, l, m, n, etc) are reserved for the original electron/hole pair basis while numerical subscripts B † j ( j = 1, 2, 3, etc) are used for the cobosn basis. The coboson basis for the electron/hole pairs is over-complete, yielding the multi-coboson states: | . . . , 2, 1 = . . . B † 2 B † 1 |g . Observable multi-exciton states are obtained via projection technique used in appendix F in order to recover the SOS blockdiagonal Hamiltonian (figure 1).
The basic elements of the coboson algebra were developed in [36] and summarized in appendix A. These are: Pauli scattering (A.1), direct Coulomb scattering (A.4) and the interaction expansion (A.11). The over-completeness of the coboson basis results in a specific form of the identity operator (A.9). The algebra is fully described by two types of coboson-coboson interactions given by tetradic matrices: Pauli exchange λ and the Coulomb scattering ξ . These carry all many-body information necessary for computing the many-exciton Green's functions and can be graphically visualized by Shiva diagrams which allow for rapid calculations and interpretation of their products and combinations [22]. These interactions can be obtained by direct comparison of the Hamiltonian and commutator relations for the cobosons with those for the uncorrelated particle bosonized Hamiltonian as described below.
Corresponding Shiva diagrams are shown in figure 2(b). The dipole coupling is obtained by transforming equation (4): In summary, the required coboson algebra parameters were calculated by projecting the system Hamiltonian in the fermion representation onto the uncorrelated electron-hole pair and coboson representations and then comparing the two. Our starting point was the fermionic Hamiltonian which describes any number of excitations in the system. We then derived the Hamiltonian in the electron-hole pair representation by truncating the fermion Hamiltonian to include up to two excitations in the system. The coboson form of the Hamiltonian was obtained by using the corresponding algebraic rules. By comparing the electron-hole pair and coboson representations, we expressed the ξ matrix of the coboson algebra in terms of U matrix of the electron-hole pair representation, which in turn was written in terms of the Coulomb matrix elements of the original electron/hole Hamiltonian. The energies/eigenfunctions of the single-exciton manifold were obtained by comparing of the action of the electron-hole pair Hamiltonian acting on this manifold. Finally, we expressed the λ matrix of the coboson algebra in terms of the electron-hole exchange P matrix. This was accomplished by comparing the action of the commutator in both pictures on the single-exciton manifold using the coboson and the fermion operators.

Green's functions and scattering matrices of cobosons
We next turn to calculating the Green's functions and exciton scattering matrices. The Green's function matrix elements in the frequency domain are defined by with η → 0 + . The main advantage of the coboson over the SOS formalism is that the expensive diagonalizations of multi-exciton blocks of H are not necessary for calculation of the Green's functions G (n) . In order to describe the coboson scattering, we start with the Bethe-Salpeter equation which describes the deviation of the cobosons from free bosons: Here G (n,0) is free-boson zero-order Green's function matrix. It is obtained by treating the excitons as structureless (λ → 0) and non-interacting (ξ → 0) bosons: where elements of the identity matrix I are equal to one for the diagonal elements and zero otherwise. The 2n-adic matrix ∆ (n) is given by (2) 12,34 = δ 4,2 δ 3,1 + δ 3,2 δ 4,1 , 123,456 ≡ δ 4,1 δ 5,2 δ 6,3 + δ 4,1 δ 5,3 δ 6,2 + δ 4,2 δ 5,1 δ 6,3 + δ 4,2 δ 5,3 δ 6,1 + δ 4,3 δ 5,1 δ 6,2 + δ 4,3 δ 5,2 δ 6,1 (25) permutes the energies in equation (23) and does not change the bosonic nature of G (2,0) . It is introduced for consistency with the modified Dyson equation derived below. Γ (n) is the coboson scattering matrix yet to be determined. Utilizing the coboson algebra rules, we can move 1 ω−H +iη in the multi-coboson Green's function (21) to the right until it eventually acts on the ground state. This will generate the Coulomb terms from the exciton scattering operator V and the Pauli blocking via the operator D. This procedure leads to the generalized Dyson equation where Ξ (n) and Λ (n) are known functions of the Coulomb direct interaction ξ and exciton exchange λ. Λ (n) accounts for Fermi statistics, which are not included in the ordinary Dyson equation. The matrix ∆ (n) makes sure that bosons of the same energy are indistinguishable. In appendix B, we worked out the details for the two coboson case (n = 2): (2) 12,34 = λ 21,43 + λ 21,34, 12,34 = ξ 21,43 .
Similarly, the case of three-coboson interactions (n = 3) was obtained in appendix C: =ξ (3,1) +ξ (3,2) +ξ (3,3) , where we have introduced several auxiliary quantities denoted by tildes:  Their physical significance is best described by the Shiva diagrams in figure 3 and is further discussed in appendix C. Note that all A (n) quantities are tensors of rank 2n. They may be mapped onto the regular matrices (rank 2) of A N n ,N n size, as described in appendix A. Comparing the modified Dyson (equation (26)) with the Bethe-Salpeter equation (22) in their matrix form and carrying out some straightforward algebra (see appendix B) yields the exciton scattering matrix Γ (n) can be mapped back onto the rank 2n tensor and is one of the key quantities of this paper.
In the coming sections, we shall apply these Green's functions and scattering matrices to calculate nonlinear optical response functions. Our protocol consists of the following steps: 1. Starting with the many-fermion Hamiltonian, derive the effective Hamiltonian and the commutation rules in the coboson representation only for the two exciton manifold. This allows one to calculate the two elementary tetradic matrices: a two exciton Coulomb scattering matrix ξ and a Pauli exclusion matrix λ. These are the only ingredients required for calculating the response functions to all orders.

Write down loop diagrams for the response function of interest and use them to construct
Green's function expressions for the response. 3. Use the Bethe-Salpeter equation to eliminate the terms representing the harmonic boson nonlinear response, whose sum must vanish identically.
4. Calculate the necessary Green's functions and scattering matrices using algebraic properties of the coboson operators. 5. Insert the Green's functions in the expressions from step 3 and finally obtain the response.
These steps are simple and, even though step 4 requires the coboson algebra, it can be carried out with minimal effort. The coboson algebra may be used for recasting these signals in SOS form (see appendix F).

Probing the two-exciton manifold by third-order signals
We first consider a continuous-wave (CW), frequency-domain experiment in which the system is subjected to three stationary beams and calculate the signal with frequency ω s = ω 1 + ω 2 − ω 3 . This is the only third-order signal allowed by our dipole selection rules. The susceptibility, χ (3) , is represented by many diagrams (for diagram rules see appendix D.2). However, we know that the free-boson contributions must cancel out. Since the free-boson system is linear, we need only consider terms which contain two-exciton states. We are thus left with only the two diagrams given in figure 4. The susceptibility, which can be read off the diagrams, is whereĜ,ˆ are the operators corresponding to the Green's functions (free-boson form) and the scattering matrix respectively;μ ± are the raising/lowering part of the dipole moment operator. Since, in CW experiments, we do not have control over the time ordering of various interactions, we will have to permute the order of the three field frequencies ω 1 , ω 2 and ω 3 . p{ω i } represents the summation over those 3! permutations. Expanding equation (34) in the one-exciton state basis yields our final compact expression for the third-order susceptibility Here the dipole moments µ are given by equation (20). It is worth noting that the dipole moments above are of the free-boson model. Justification of this non-trivial fact is given in appendix E. The scattering matrix (2) is given in equations (33), (27). Note that the indices of (2) appear in the reverse order of their respective states along the loop; this follows from of our convention for identifying matrix elements of the Green's function (equation (B.1). We have further introduced the following abbreviated notation for the single-exciton free-boson Green's function: Loop diagrams for the third-order susceptibility. These diagrams are frequency domain counterparts to diagrams III a and III b that appear later in figure 6. The other two anharmonic diagrams, I a and II a , are both identical to III b in the frequency domain because the order of interactions is not controlled. Equation (34) can be read directly off these diagrams.
where γ 1 are the pure dephasing rates for the coboson states. In the frequency-domain, these become We further introduce similar abbreviated notation for the frequency-domain free-boson Green's function for two and three excitons (the latter will be used in the next section): In these expressions, all Green's functions represent free-boson evolution in the relevant exciton-manifold as given by the number of indices; G 1 , G 21 and G 321 represent evolution in the single-, double-and triple-exciton manifolds respectively. We next turn to time-domain techniques. We first consider an experiment in which the system is subjected to a sequence of three, temporally well-separated pulses with wave-vectors k 1 , k 2 , k 3 (in chronological order). The third-order polarization is then created along eight possible directions: where k s = ±k 1 ± k 2 ± k 3 is the signal direction; ω s = ±ω 1 ± ω 2 ± ω 3 its frequency; r the average distance from the system to the detector; τ 4 the detection time; and t 3 , t 2 , t 1 are the time intervals separating the interactions with optical fields. The nonlinear polarization is related to the response function S (3) s via (42) We further recast the signal in the frequency domain The k 1 + k 2 + k 3 signal vanishes due to the dipole selection rules in our model. We further have the symmetry . This leaves us with three independent techniques k I = −k 1 + k 2 + k 3 , k II = k 1 − k 2 + k 3 and k III = k 1 + k 2 − k 3 .
After drawing the relevant diagrams, we can write the expressions for the corresponding responses directly following the diagram rules given in the appendix D.1. These rules suffice for all third-order diagrams considered in this paper. Higher order (e.g. fifth and above) signals can be derived similarly by writing the time-domain expressions from the relevant diagram and performing the time-convolutions analytically. Note that the dipole operators (and the states they create or annihilate) are labeled by the order in real-time in which they operate.
In appendix E, we apply these rules in order to derive the third-order signals shown in loop diagrams of figures 5 and 6. For the k I technique (figure 5) the contribution from diagrams I b and I c diagrams vanishes (compensated by I a 0 1 , I a 0 2 ) as from harmonic oscillator. The remaining diagram (I a ) reads Similarly, the only diagram contributing to k II is shown in figure 6 Two diagrams contribute to the k III technique (diagrams III a , III b in figure 6 (lower panel)) and the resulting polarization is given by Equations (44)-(46) are fully specified by the structure of Γ (2) in equations (33) and (28). They give succinct expressions for the various coherent signals obtained from a time-domain experiment of temporally well-separated pulses. Similar equations were obtained for the Frenkel exciton model in [19]. For completeness, these quasiparticle expressions are compared to the SOS expressions given in appendix F. The effects of coupling to a bath can be incorporated in the above signals by properly averaging the relevant products of Green's functions for each time interval. This is described in appendix G (figure 7).

Probing the three-exciton manifold by fifth-order signals
The general expressions for signals resulting from a series of temporally ordered pulses become cumbersome as one proceeds to higher-order optical processes with multiple time intervals. Nevertheless, the coboson expressions can readily be derived. This will be demonstrated below for four distinct five-pulse techniques.
We start with the k s = 3k 1 − 2k 2 signal that results from three simultaneous interactions 3k 1 followed by evolution in the three-exciton manifold, then two simultaneous interactions 2k 2 followed by further evolution and then finally signal detection. This technique provides a clean look at the three-exciton manifold. In the rotating wave approximation, it is given by the three We may then easily read the response off the three diagrams in terms of correlation functions following our rules. This gives . Bottom: loop diagrams for the technique k s = +2k 1 + k 2 − 2k 3 which correlates the double-and triple-exciton manifolds. Note that these three diagrams also differ only by their sign which is determined by the number of interactions from the right. Thus, diagram S (5) β may be interpreted to cancel diagram S (5) γ leaving equation (54).
Note that diagram S (5) b is split into two in the same way as diagram III b of figure 6 (lower panel); however, this time both expressions are identical.
We next turn to the technique k s = +k 1 + 2k 2 − 2k 3 in which the signal pulse k s coincides temporally with −2k 3 . This is given by three diagrams shown in figure 9 (top). The diagrams differ only in their overall sign which is determined by (−1) N where N is the number of interactions from the right. Hence, the contribution from diagram β cancels that from γ (or α). We may thus write the response in the time and frequency domains as This technique probes correlations between the single ( 1 ) and triple ( 2 ) exciton manifolds.
We now turn to a third technique, k s = +2k 1 + k 2 − 2k 3 ( figure 9, bottom). The diagrams also differ only in their overall sign in exactly the same way as the diagrams from the previous technique. Thus, once again there is a cancellation of the contributions from diagrams β and γ (or α). The time and frequency-domain responses are This experiment probes correlations between the double (in 1 ) and triple (in 2 ) exciton manifolds. Finally, we consider the fifth-order signal generated at k s = +2k 1 − 2k 2 + k 3 . This technique is more complicated but can also look at two-exciton populations | f f |. The signal is given by the four diagrams shown in figure 10  (III b0 , III b1 and III b2 ). III b0 contains the non-interacting part of the two-exciton Green's function. III b1 is the part of the Dyson equation with the scattering matrix, when τ < t 3 , while III b1 corresponds to the case when τ > t 3 . Exciton annihilation of | f f | can be added to these expressions during t 2 . This would be an interesting future extension (figures 11 and 12).

Numerical examples and discussion
We now demonstrate how the quasi-particle coboson formalism works in practice by applying it to a model system. The implementation is based on an approximate, tight-binding electronic Hamiltonian suited to describe single-and double-excitations in molecular aggregates [28]. We assume a linear chain of n coupled two-level systems separated by distance r . Electron spin is neglected. Double quantum coherence (DQC) signals are our technique of choice since they can be calculated by the SOS expressions without invoking the coboson formalism. We shall compare our results with those obtained in uncorrelated electron-hole basis by NEGF or SOS methods. We start with the Hamiltonian of single-excitations and double-excitations in the site basis given by equation (1). Simulations are presented for N = 2, 3, 4 atomic sites. The single-exciton manifold is constructed as a linear combination of on-site excitations (Frenkel) and charge-transfer excitations. The two-exciton manifold is represented by a basis set of non-interacting bosons. Pauli-allowed physical states are only a subspace of the two-exciton manifold. For the molecular dimer (N = 2) we represent the manifolds graphically in figure 14.
As indicated earlier, the coboson basis for the two-exciton manifold is over-complete and nonorthogonal; this poses some difficulties in the SOS formalism by making the doubly excited Hamiltonian blocks H (2) singular (see appendix F for details). More precisely, the block is ill-conditioned. For bright, on-site excitations we assume a transition moment µ = 1 au and charge-transfer excitations are assumed to be dark (µ = 0 au). In other words, the transition dipole moment is proportional to the overlap between electron and hole wave functions. The sizes of the tetradic matrices λ and ξ are 4 4 , 9 4 = 6561 and 81 4 = 4.3 × 10 7 for n = 2, 3 and 4, respectively. These are also the sizes of the extended doubly excited manifolds. The electron and hole hopping amplitudes in our model t (1) m1,n1 and t (2) m2,n2 are assumed to decay exponentially with the separation r between chromophores t 1 m1,n1 = t 2 m2,n2 = A * exp(−|r m − r n |/ρ) for m = n. The diagonal contributions are given by 5.0 * k e /η for m = n. Loosely speaking, this is an extension of the tight-binding model. The Hamiltonian parameters are summarized in table 1. In all calculations of the two-dimensional (2D)-DQC, spectra are based on the results of section 4 and we assume a common dephasing rate of single excitons γ i = 400 cm −1 . The Coulomb attraction as well as electron-electron and hole-hole repulsion are assumed to have the truncated Coulomb form V (i) k1,n1 = W (1,2) n1,m2 = (k e /η)/(|r k − r n |/η + 1) where η is a cutoff distance. The modification of the denominator accounts for possible residence of the interacting particles on the same site which is necessary for the extended basis and is automatically removed by the Pauli blocking of those states. The tetradic matrices V (1) l1k1,n1m1 , V (2) l1k1,n1m1 , W (1,2) l1k1,n1m1 describe electron-electron and hole-hole repulsion and electron-hole attraction, respectively. Assuming that the atomic orbitals are strongly localized [29], we define these tetradic matrices as m2 . The parameters of the electron-hole Hamiltonian (equation (1)) defined above allow us to calculate the elements of the electron-hole pair Hamiltonian (equation (6)). The choice of parameters results in a single electron-hole pair block of the Hamiltonian (equation (7)) with dominant diagonal coupling (t) and weak off-diagonal W contributions. On the other hand, the two exciton block given by equation (8) is sparse. The single electron-hole pair block can be diagonalized numerically (equation (11)) yielding the eigenfunctions 1|m . The Pauli exchange matrix can be projected on those states as in equation (19) yielding λ. The Coulomb part of the two exciton block projected on those states via equation (16) gives the Coulomb direct scattering ξ . Together, these tetradic matrices constitute the coboson algebra.
Before discussing any measurable signals, we shall briefly survey the structure and properties of the main ingredients required for their simulations: the Green's function and the scattering matrix. Specifically we shall discuss the deviation of the cobosons from ideal bosons on the level of the spectral function defined via the imaginary part of Tr(G (2,0) (ω)) for bosons and Tr(G (2) (ω)) for cobosons. The repeated actions of Pauli blocking and Coulomb scattering is demonstrated for a model dimer of coupled two-level systems in figure 13. The top panel depicts the two-boson spectral function. Here, resonances appear at the sum of the energy of the harmonic oscillators (arising from U (1) ). Those are a useful measure of the number of possible permutations leading to a given bi-exciton like formation composed of two harmonic oscillators. Due to the bosonic properties, all possible permutations of double-excitations contribute. The lowest energy peak corresponds to bi-excitons composed of two subsequent charge-transfer excitations ( figure 14). The high-energy peak correspond to two bright, on-site excitations. In between, we see all combinations of two subsequent excitations leading to a double excitation with two particles residing on the same site. The action of the Pauli-exclusion matrix Λ (2) is only demonstrated in the middle panel of figure 13. It removes all double excitations with two particles residing on the same site leaving only two resonances (two-exciton states) in the spectral function. Coulomb interaction between cobosons is taken into account in the bottom panel of figure 13. Due to the anharmonic contributions contained in Γ (2) (via Ξ (2) as the projection of U (2) ), the bright, two-exciton resonance degeneracy is lifted.
We can now turn to the DQC spectra depicted in figure 15. The signal vanishes for pure bosons. Thus, it provides a good measure of the anharmonicities which enter via two-and single-exciton coherences as cross-peaks. Formal analysis of the spectra is done by using the SOS expression given in appendix F. For a dimer, the effect of Pauli scattering is only shown in the top row of figure 15, as real (left) and imaginary (right) values. A single, two-exciton cross-resonance along the 2 -axis appears exactly at twice the energy of the bright singleexciton state (given along the 3 -axis). Coulomb scattering effects, as contained in Ξ (2) , are demonstrated in the second row of figure 15. The energy of the two-exciton resonance along 2 is reduced as compared to that of double single-exciton transitions showing formation of a weakly bound bi-exciton. The cross-peak is split into two along 3 showing the splitting of the fundamental g → e ground-state to single-exciton transition from the e → f single-to bi-exciton transition. The second possible two-exciton resonance (the low energy pole in the coboson spectral function) is dark in the DQC signal. The signals for N = 3 and 4 are depicted in the third and fourth row of figure 15 respectively. Due to increased exciton de-localization in can be identified for N = 4. Now the g → e fundamental appears red-detuned from the e → f single-to bi-exciton transition.

Concluding remarks
We have developed a quasi-particle coboson, diagrammatic technique for computing high-order, nonlinear response functions. The approach is based on the Bethe-Salpeter and modified Dyson equations which define the coboson scattering matrices and relies solely on the coboson algebra. It does not suffer from the massive pathway cancellations that complicate the conventional SOS expressions. The algebra is completely defined by two-exciton interactions and may be readily derived from the original fermionic Hamiltonian. For any given number of excitons, this makes it possible to calculate high-order nonlinear response in terms of a sparse coboson scattering matrix Γ (n) . Applications are made to the third-and fifth-order responses which probe the twoand three-exciton manifolds, respectively. The single-exciton manifold, which can be extracted from standard calculations such as TDHF, TDFT [30][31][32][33][34], or the Bethe-Salpeter equation for fermions [35], fully defines the algebra. Many-body interactions are fully accounted for using products and combinations of two-body interactions: direct Coulomb scattering ξ , and Pauli exchange λ. Cobosons are more closely connected to optical measurements than the original fermion operators since each transition can be described by a single coboson rather than two fermions. In addition, taking into account the finite bandwidth of laser pulses, it should be possible to make drastic truncations of the relevant cobosons in a given measurement. This is not possible in the fermion representation and makes the calculations much more affordable even at the SOS level. • Coulomb (direct) scattering: • Interaction expansion: • Overcompleteness: • P the electron-hole exchange tensor [18]: lk,nm = 1 2 δ m2,l2 δ n2,k2 δ m1,n1 δ k1,l1 . Tetradic/hexadic matrix notation. In Combescot's work and elsewhere, matrix elements are sometimes specified by Those are rank 2n tensors and can be mapped on the regular (rank 2) matrices by adopting the notation for an arbitrary matrix A row,col where row and col are given as elements of n ⊗ · · · ⊗ 1. Practically, we construct the tensors in Mathematica as the following this yields the BlockMatrix whose elements are 2n-2-rank tensors and so on till the regular matrices. Function ArrayFlatten provides the mapping to the regular matrix manifold. The back-mapping to the tensor is provided by

Appendix B. The two-exciton Green's function and the scattering matrix
Here, we derive the Green's functions, the exciton scattering matrices and transition dipole moments in the coboson representation. For clarity, we present the main results first and then elaborate on the significance of the various terms and parameters. This section is essential for the subsequent calculations of all nonlinear signals. The key matrix equation to obtain the n 2 coboson Green's functions G (n) is the modified Dyson equation (26). Compared with the Bethe-Salpeter equation (22), this yields the exciton scattering matrix (equation (33)) whose components are the main subject of this appendix.
We start with the two-exciton manifold (n = 2). The matrix elements of the two-exciton Green's function are defined via The zero-order coboson Green's function G (2,0) is obtained by treating the excitons as structureless and non-interacting bosons The full form of the two-coboson Green's function requires the identity which is obtained by (i) moving B 3 to the right to finally act on the ground state and (ii) using the Pauli scattering commutation relations (A.1). By means of the interaction expansion in appendix A, we can shift that it finally acts on the ground state. After some straightforward but lengthy algebra utilizing equation (A.7), this procedure leads to By utilizing equation (B.3) the above expression (B.4) can be recast in the more compact matrix notation of modified Dyson equation (26) with parameters given by equation (27) From equation (19), we can see that λ is symmetric to simultaneous interchange of both pairs of indices (i.e. λ 21,43 = λ 12,34 ). The Green's function (equation (B.4)) assumes a simpler form when ξ vanishes (as in small enough quantum dots) and all nonlinearities are caused by Pauli exclusion alone: G (2) = G (2,0) (∆ (2) − Λ (2) ). From the form of the Λ (2) matrix, it is clear that the Green's functions vanish when the two excitons try to occupy the same state (i.e. G (2) 11,11 = 0). Therefore the Pauli exchange matrix is responsible for removing those unphysical states caused by over-completeness of the coboson basis.

Appendix C. The three-excitons Green's function
The three-exciton Green's function is defined by At this point we introduce several axillary quantities, expressed in the form of sixtic matrices. The first is the non-interacting exciton Green's function The second is the Pauli blocking indicator (PBI) matrix, similar to equation (B.3) with the parameters given by equations (25) and (30). Equation (C.3) has been obtained by moving B 4 in equation (C.1) all the way to the right so that it finally acts on the ground state. Pauli scattering (equation (A.1)) is used and the corresponding Shiva diagrams are given in figure 2.
When the Coulomb coupling is neglected, the Green's function becomes The matrix elements ofλ can be viewed as flags switching on and off the Green's function matrix elements in order to maintain Pauli blocking and conservation of energy. The first term in the PBI (equation (C.3)) assures energy conservation in the absence of external perturbations. That is, only diagonal elements are flagged on. There are 3! ways to exchange the excitons and still conserve energy. The second PBI term can switch off some of those flags in accordance to Pauli blocking modeled by single-exchange (two-coboson) elements λ. This occurs when two cobosons try to share the same site. The latter case is described by the part of the λ matrix with n1 = n2 and l1 = l2 so that λ 12,34 = δ 1,3 δ 2,4 . However in the case of three cobosons we have 12 (seeλ (0) of equation (30)) uncompensated negative terms. These terms are compensated by the third part of the PBI of equation (31). Contribution fromλ (1) restores those negative flags back to zero making the exchange via virtual excitons in the three coboson manifold. The net Green's function thus vanishes when the underlying electrons and holes have the same quantum numbers and sites.
In the general case (retaining the Coulomb interactions), after lengthy but straightforward application of the algorithm for shifting the Green's function operator 1/(ω −Ĥ + iη) all the way to the right and using the interaction expansion, we obtain coupled system of Dyson equations for G (3) and an auxiliary matrixG (3) : is given by equations (32). We next use the hexadic matrices multiplication convention in appendix A in order to recast equations (C.6), (C.5) in a compact matrix form: with the parameters˜ (3) and (3) given by equations (28)- (30). The Shiva diagrams in figure 3 illustrate the single and two coboson interactions involved in the above system of equations.
We can excludeG (3) from equation (C.8) by recasting equation (C.7) in the form This yields a single Dyson equation for the three-exciton Green's function in equation (26) with parameters given in equation (28).

Appendix D. Algorithm for computing the fifth-order techniques
1. Write the dipole operators, Green's functions and any relevant scattering matrices, in the order in which they appear on the loop as we move clockwise per the general instructions given in section D.1. 2. For all periods of evolution during which multi-exciton manifolds are not accessed, simply replace the time arguments of the Green's functions and any scattering matrices for that evolution with their Fourier-conjugate frequencies when they appear on the left side of the loop and with the negatives of their conjugate frequencies when on the right side. As an example, compare the t 1 and t 2 evolutions in diagram I a of figure 4 with the corresponding Green's functions in equation (44). 3. Periods of evolution during which a scattering matrix is active on one branch of the loop and an interaction occurs on the other must be treated differently. In these cases, the Green's functions and scattering matrices (if any) corresponding to the state resulting from the interaction have their arguments shifted by the Fourier frequency that is conjugate to the time period during which the original scattering matrix was active but before the interrupting interaction. For a concrete example of how this works, observe the shifting of the arguments of the Green's functions and equation (49) that correspond to the bra side of the t 3 evolution and t 2 evolution in diagrams III b and S (5) b respectively in figures 5 and 6. Finally, the arguments of the Green's function and scattering matrix corresponding to the 'interrupted' multi-excitonic evolution are simply the Fourier-conjugate frequency to the period of time before the interrupting interaction (compare the t 2 evolution and t 1 evolution in diagrams III b and S (5) b respectively in figures 5 and 6 with equation (49)). 4. For all periods of multi-exciton evolution along the loop during which the density matrix is in a coherence with the single-exciton manifold (e.g. the coherence with the single-exciton state 1| in diagram I a of figure 4), add the energy and dephasing of the relevant singleexciton state to the Fourier-conjugate frequency of the Green's functions and scattering matrix corresponding to the multi-exciton time period. This incorporates the evolution of the single-exciton state; therefore, we then omit the Green's function corresponding to its evolution. Note that this step applies also to those cases described in step 3 for which the interrupting interaction takes us to the first excitonic manifold as in diagram III b of figure 5.

D.1. Loop diagram rules in the time domain
TD1. The loop represents the density operator. Its left branch stands for the ket, the right corresponds to the bra. TD2. Each interaction with a field mode is represented by an arrow pointing to the left or right and acting on either the ket or the bra. TD3. Arrows pointing to the right represent field annihilation operators and arrows pointing to the left represent field creation operators. This is made explicit by dressing the arrow with appropriate wave vectors ±k j . TD4. Within the rotating wave approximation, each interaction with a field annihilation operator is accompanied by the dipole operator µ † , which leads to the excitation of the state represented by the ket and de-excitation of the state represented by the bra. Arrows pointing 'inwards' (i.e. pointing to the right on the ket and to the left on the bra) consequently cause absorption of a photon by exciting the system, whereas arrows pointing 'outwards' represent de-exciting the system by photon emission. TD5. The interaction at the observation time is always the last. As a convention, it is chosen to occur from the left. This choice is arbitrary and does not affect the result. TD6. The overall sign of the correlation function is given by (−1) N R where N R is the number of interactions from the right. TD7. Since the loop time goes clockwise along the loop, periods of evolution on the left branch correspond to retarded Green's functions while periods of evolution on the right branch correspond to advanced Green's functions. FD3. Arrows pointing to the right represent field annihilation operators and arrows pointing to the left represent field creation operators. This is made explicit by dressing the arrow with appropriate wave vectors ±k j . FD4. Within the rotating wave approximation, each interaction with a field annihilation operator is accompanied by the dipole operator µ † , which leads to excitation of the material system. Arrows pointing to the right cause absorption of a photon by exciting the molecule, whereas arrows pointing to the left represent de-exciting the system by photon emission. FD5. The interaction at the observation time is fixed to be with the detected mode and is always the last. It is chosen to occur on the left branch of the loop. This choice is arbitrary and does not affect the result. FD6. Since the loop time goes clockwise along the loop, periods of evolution on the left branch correspond to retarded Green's functions while periods of evolution on the right branch correspond to advanced Green's functions. FD7. The frequency arguments of the various Green's functions are cumulative. That is, they are given by the sum of all 'earlier' interactions along the loop. FD8. A diagram carries a factor of (−1) N R where N R is the number of interactions from the right.

D.2. Loop diagram rules in the frequency-domain
Examining first the third-order response in the k I technique, it is clear from the diagrams that the excited-state emission (ESE) and ground-state bleaching (GSB) processes (diagrams I b and I c respectively of figure 5) never have more than one exciton. Their evolution is thus harmonic. However, the excited-state absorption (ESA) process (I a in figure 5) includes a period where two excitons are present. We will use the Bethe-Salpeter equation for the two-exciton Green's function. The first term in this equation represents free-boson evolution in the twoexciton manifold (represented by diagram I a 0 in figure 5) while the second term represents the anharmonicity. Since the nonlinear response from a harmonic oscillator vanishes, we know that the harmonic evolution of the ESA process must cancel with the sum of the ESE and GSB processes. We may thus express the nonlinear response in the k I direction by a single term (diagram I a of figure 5) Similarly, in the k II technique, a harmonic cancellation occurs and the signal is merely given by the anharmonic portion of the only diagram to include multiple excitons (diagram II a in figure 6). Hence For the k III signal, both diagrams access the two-exciton manifold. Still, since the harmonic portion must cancel, the signal may be written as the sum of the anharmonic contributions from each diagram (III a and III b in figure 6): To evaluate these correlation functions, we simply insert the resolution of the identity matrix given in the appendix (equation (A.9)) between all Green's functions and dipole operators. First however, we will use the fact (shown in equation (F.9)) that, for the two-exciton dipole moments, µ = µ (0) (I − λ), where, according to our convention, µ (0) is the dipole for a harmonic oscillator.
Using the over-completeness relation defined in the appendix and the definition of the Green's function defined in equation (B.1), it is easily verified that G (2) = −λG (2) . Since the two-exciton Green's function enters sandwiched between two two-exciton dipole operators, which may be simplified via ) This factor of 4 will be canceled by the factor of 1 4 introduced by one of the two-exciton blocks of the identity operator when it is inserted into the above correlation functions (the two-exciton block of the identity is inserted on both sides of G (2,0) Γ (2) G (2,0) ). Thus, we may simply replace the dipole moments with their harmonic counterparts in the above expressions for the signals. Note that there are two dipoles accessing the two-exciton manifold, each contributing a factor of √ 2, the result is a factor of 2 for each of the above signals. Since the remaining dipoles and Green's functions are in the single-exciton space, we may evaluate simply by inserting the identity matrix between the operators in the correlation functions above. This requires the evaluation of matrix elements: (G (2,0) (2) G (2,0) ) 12,34 = G (2,0) Since the result is to be plotted with respect to two of the frequencies conjugate to the time arguments, we will Fourier transform the signals according to P I,II ( 3 , t 2 , 1 ) = P I,II (t 3 , t 2 , t 1 ) e i 3 t 3 e i 1 t 1 dt 1 dt 3 , (E.6) This is shown in detail in sections E.1 and E.2 and results in equations (44)-(46). Turning now to the fifth-order diagrams under consideration in section 5 and shown in figure 6, the time-domain response functions are easily read off from the rules in appendices B and C: ×ˆ (2), † (τ 2 − τ 2 )Ĝ (2,0), † (τ 2 )μ −Ĝ(3,0) (t 1 + t 2 − τ 1 ) Following the same method as for the third-order signals and evaluating these correlation functions, one inserts the identity given by equation (A.9) between all operators in the above correlation functions. The result is then Fourier transformed with respect to the two time intervals. Once again, the convolutions can be carried out analytically by making the substitution in equation (E.12). These manipulations result in equations (48)-(50).

E.1. Derivation of the k I
We begin with equation (E.1) (see figure 5) and perform the simplification of the dipole moments to those of single-excitation described above. This crucial simplification for the scattering matrix approach also holds for NEE formalism in [27] thus not being specific for the coboson representation. Expanding equation (E.1) in the coboson basis yields In order to carry out the convolution integrals, we Fourier transform the scattering matrix according to (2) and explicitly write out the Green's functions as exponentials: Since each forward (backward) propagation through time is accompanied by a factor −i (i), an overall factor of (i)(−i)(−i) is introduced. In the above equation the time depend integrals can be carried out analytically as . (E.14) Formally, our next step will be to Fourier transform with respect to t 3 and t 1 via In anticipation of our subsequent integration over dω by closing the contour in the upper half plane (which will introduce a factor of 2π i by the residue theorem), we will drop the first, second and fourth terms in equation (E.14) on the grounds that their poles are all in the lower half plane. The only term which will contribute therefore is .
Note that the residue integral and the Green's functions together introduce a factor of (i)(i)(−i)(−i) to the numerator of the above expression while the integrals over τ , τ each introduce a factor of i to the denominator and the integrals over t 1 , t 3 each introduce a factor of −i to the denominator. Thus, all these imaginary factors cancel and we are left with equation (E.16). There is a lone pole in the upper half plane at ω = E 1 + 2iγ 1 + 3 (see figure 11(B)). The origin of this pole is the factor of e 2iγ 1 t 3 introduced by the product . We may thus carry out the dω integration, arriving at Because the signal is related to the polarization by equation (42), we may obtain This immediately gives equations (44)-(46) are obtained similarly.

E.2. Derivation of k III
Due to cancellation of the diagrams in figure 12 (III a0 cancels out III b0 and III a1 cancels out III b1 ) we obtain the signal in the time domain This form of the expression is basis independent. At this point, we introduce the coboson algebra by inserting the identity operator in equation (E.9) and, using equations ((E.4), B10, (36), (37)), we obtain the signal in coboson representation: In the above expression we have used the fact that G 43 (t 3 + t 2 − τ ) = iG 43 (t 2 + t 3 )G 43 (−τ ).

III
( 3 , 2 , t 1 ) = −i e 2 ,e 1 , f µ g,e 1 µ e 1 , f µ * e 2 , f µ * g,e 2 G e 2 (t 1 )G f ( 2 + ω 1 + ω 2 )G e 1 ( 3 + ω 1 + ω 2 − ω 3 ) +i e 2 ,e 1 , f µ g,e 1 µ e 1 , f µ * e 2 , f µ * g,e 2 G e 2 (t 1 )G f ( 2 + ω 1 + ω 2 )G f e 1 ( 3 + ω 1 + ω 2 − ω 3 ), where G e = G 1 and G f,e = 1/(ω − E f + E e + iγ ). These expressions are easily read off the corresponding diagrams in figures 3 and 4. Note that the states e i in the SOS formula are equivalent to the states i from the quasiparticle expressions since they were obtained in either case from diagonalization of the one-exciton manifold. Here, we use the notation e 1 for these states to reinforce that the expressions are in the SOS picture. Unlike in the quasiparticle picture, implementing the SOS expressions requires diagonalizing the second exciton manifold in order to obtain the two-exciton energies {E f }. The major effort is now shifted to computing the energies {E f }. Once these are known, the expressions for the nonlinear response are very simple.
In the coboson representation, those are given by diagonalization of the n-exciton block. The single-exciton block (diagonalized in exciton basis |1 ) as well as two-and three-exciton blocks are = (n) − (n) and λ (n) is the norm of the matrix defined in a regular way as the maximum of its singular values. The above blocks can be written in a compact matrix notation as λ (n) H (n) = G (n,0) (ω = 0) −1 + (n) λ (n) .
In order to obtain the nonlinear third-order response functions, we also need to know the dipole moments between the ground state and the single-exciton manifold as well as from the single-to double-exciton manifold. The single-exciton dipole moment can be readily obtained from equation (20) µ e = g|B 1 2 µ 2 B † 2 |g = µ 1 . (F.8)