where |k⟩ is the vacuum's Fock state (eigenvector of Eq. 36), and Îex = |G⟩⟨G| + m |m⟩⟨m| + m̸ =n |m, n⟩⟨m, n| is the identity operator for the electronic exciton subspace up to double excitation. The decay rate associated with L is denoted as Γ. Note that the jump operator in Eq. 43, in principle, includes all possible transitions between Fock states |n + 1⟩ and |n⟩. In a practical simulation, one can choose to include the physically relevant decay channels. In the model calculation presented in this paper, we will only consider one decay channel (see Eq. 44). The system Lindblad jump operator in this subspace of Fock states can be written as