\section{Decomposition of tensors}


Here, we will discuss the decomposition of $\tilde{M}$ and $\tilde{E}$ into products of lower dimensional tensors and the underlying ideas.

We start with $\tilde{M}$ and its decomposition:
\begin{align}
\tilde{M}_{p s k\alpha_{p}\alpha_{s}l}^{k' \alpha'_{p}\alpha'_{s}l'}(n)=&\sum_{\substack{\tilde{s}_1 \tilde{p}_1 \\ \tilde{p}_2 m_2 m_t}} \tilde{M}^{(1),\alpha'_{p}}_{\tilde{p}_2,\alpha_{p}m_t m_2} \tilde{M}^{(2) k' \tilde{p}_2 \alpha'_{s}l'm_t}_{\tilde{s}_1  \tilde{p}_1  k\alpha_{s} m_2}(n)
\tilde{M}^{(3)l'}_{m_2 m_t l} \tilde{M}^{(4)}_{\tilde{s}_1 s m_t}\tilde{M}^{(5)}_{\tilde{p}_1 p m_t}\tilde{M}^{(6)}_{s p m_t}. \label{m_decomp}
\end{align}
As seen in Eq.~\eqref{mtildetensor} the indices $k$ and $k'$ discriminate the different interactions handled in the tensor $\tilde{M}$. If $k$ is not changed ($\delta_{k,k'}$) $\tilde{M}$ does not change anything and just passes the information through present in the indices $\alpha_p$ $\alpha_s$, $k$ and $l$, in this case the interaction, if any,  happens in a subsequent tensor. If the index $k$ is changed from $1$ to $-1$ in $k'$, the double time integrated correlation function in $G$ describes a first system-bath interaction in the interval of the prior $\tilde{E}$ tensor, which is then finalized in the time interval connected to the $\tilde{M}$ tensor. Then the index $l$ passes the information, if the system-bath interaction happened on the left or right side of the Liouville space in $\tilde{E}$ to the tensor $\tilde{M}$, whereas $\alpha_s$ and $\alpha_p$ pass information about the changed system Hilbert space states described in $\tilde{E}$ to the tensor $\tilde{M}$. Therefore $\alpha'_s$ and $\alpha'_p$, which are passed to subsequent tensors $\tilde{M}$, are set to $1$ (note we use indices for s and p starting with $1$), which means the interaction already occurred and no subsequent interaction can occur on one of the remaining $\tilde{M}$.

The tensors used for the decomposition of $\tilde{M}$ reduce the number of indices per tensor and dimension of each of the indices.
We start with $\tilde{M}^{(1)}$:
\begin{align}
\tilde{M}^{(1),\alpha'_{p}}_{\tilde{p}_2,\alpha_{p}m_t m_2}=&(\delta_{\alpha_{p},\alpha'_{p}}\delta_{\tilde{p}_2,1}\delta_{m_t,3}\delta_{m_2,1}
+\sum_{m_1=1,2}\delta_{m_2,1}\delta_{m_t,m_1}\delta_{\tilde{p}_2+1,\alpha_{p}}\delta_{\alpha'_p,1}).
\end{align}
 The value of the index $m_t$ determines 
whether $\tilde{M}$ does not act on the system ($m_t=3$) and thus just passes  $\alpha_{p}$ through to $\alpha'_{p}$, or whether  $\alpha_p$ is passed via the index $\tilde{p}_2$ to the integrated correlation function $G$ ($m_t=1,2$). $m_1$ indicates if the second interaction is on the left or right side and is passed to $m_t$.
In the next tensor $\tilde{M}^{(2)}$, we include the double integrated correlation function $G$:
\begin{align}
\tilde{M}^{(2),k',\tilde{p}_2,\alpha'_{s}l'm_t}_{\tilde{s}_1, \tilde{p}_1, k,\alpha_{s} m_2}(n)=&  \left(\delta_{\alpha_{s},\alpha'_{s}}\delta_{\tilde{s}_1,1}\delta_{\tilde{p}_1,1}\delta_{\tilde{p}_2,1}\delta_{m_2,1} \delta_{m_t,3}\delta_{k,k'}\right.
+\sum_{m_1=1,2}\sum_{\tilde{s}_2}\delta_{\tilde{s}_2+1,\alpha_{s}}\delta_{\alpha'_s,1} \delta_{m_t,m_1}\delta_{k',-1}\delta_{k,0}\nonumber\\ &\left.\qquad \cdot G_{\tilde{s}_1,\tilde{p}_1,\tilde{s}_2,\tilde{p}_2 m_1m_2}(n)\right). \label{G_tensor1}  
\end{align}
Again, for $m_t=3$ the incoming $\alpha_{s}$ is just passed through to $\alpha'_{s}$.
Whereas for $m_t=m_1=1/2$  a second system bath interaction on the left/right side in Liouville space in $G$ occurs. Therefore $\alpha_s$ is connected to $G$ via the index $\tilde{s}_2$.
$G$ is included in $\tilde{M}^{(2)}$ and its other indices are passed to the other tensors in the product: $\tilde{p}_2$ to $\tilde{M}^{(1)}$, $\tilde{p}_1$ to $\tilde{M}^{(5)}$, $\tilde{s}_1$ to $\tilde{M}^{(4)}$, $m_2$ to $\tilde{M}^{(3)}$.
The next tensor $\tilde{M}^{(3)}$ handles the distribution of the $l$ index:
\begin{align}
\tilde{M}^{(3)l'}_{m_2 m_t l}=&\delta_{m2,1}\delta_{l,l'}\delta_{m_t,3} +\sum_{m_1=1,2} \delta_{m_1,m_t}\delta_{m_2+1,l}\delta_{l',1}
\end{align}

As usual, for $m_t=3$, $l$ is passed through to $l'$, and for $m_t=m_1=1,2$  $l$ is passed to  $G$.
The next tensors starting with $\tilde{M}^{(4)}$ connect indices coming from $G$ either with the left or right part of the indices $s$, $p$: 
\begin{align}
\tilde{M}^{(4)}_{\tilde{s}_1 s m_t}=& \delta_{m_t,1} \delta_{s_L,\tilde{s}_1}-\delta_{m_t,2} \delta_{s_R,\tilde{s}_1}+\delta_{m_t,3}\delta_{\tilde{s}_1,1}, 
\end{align}
so for $m_t=1$ $\tilde{s}_1$ is connected to $s_L$ and for $m_t=2$ $\tilde{s}_1$ is connected to $s_R$. $m_t=3$ is the pass-through mode and effectively does nothing.
$\tilde{M}^{(5)}$ is constructed in an analogous way:
\begin{align}
\tilde{M}^{(5)}_{\tilde{p}_1 p m_t}=& \delta_{m_t,1} \delta_{p_L,\tilde{p}_1}-\delta_{m_t,2} \delta_{p_R,\tilde{p}_1}+\delta_{m_t,3}\delta_{\tilde{p}_1,1},
\end{align}
since for $m_t=1$ $\tilde{p}_1$ is connected to $p_L$ and  for $m_t=2$ $\tilde{p}_1$   to $p_R$.
The next tensor $\tilde{M}^{(6)}$  connects  $s$ and $p$ either to the left or right part (on opposite site to the previous tensors!), and is controlled again by $m_t$:
\begin{align}
\tilde{M}^{(6)}_{s p m_t}=&\delta_{m_t,1}\delta_{s_R,p_R}+\delta_{m_t,2}\delta_{s_L,p_L}+\delta_{m_t,3}.
\end{align} 
Thus for $m_t=3$ nothing happens, and for $m_t=1$ the right sides of $s$ and $p$ are connected, whereas for $m_t=2$ the left sides are connected.

Before we discuss the decomposition of $\tilde{E}$, we first discuss its contributions from Eq.(\ref{etildetensor}):
For $k=0$ no interaction is happening in the interval $\Delta T$, so we do not have to pass information to the $\tilde{M}$ tensors.
On the other hand $k=-1$ means, so far nothing is applied to the time interval, the contribution will occur in one of the $\tilde{M}$, which are later added to interval $\Delta T$, therefore $k=-1$ will not be included in the summation at the edge of network.
Also $k=2$ means no system bath interactions in subsequent $\tilde{M}$ inside $\Delta T$, which are connected to other $\Delta T$ intervals.
Finally for $k=1$  one system-bath interactions occur in the $\Delta T$ interval associated with $\tilde{E}$, whereas the second occurs in another $\Delta T$ interval involving its $\tilde{M}$ tensor. Here $l$ encodes the side of the Liouville space (left or right), on which the system bath interaction occurs.

Again the tensors used for the decomposition of $\tilde{E}$ into low rank and dimensional tensors read:
\begin{align}
&\tilde{E}_{s,p,k}^{\alpha_{p},\alpha_{s}l }= \sum_{m_1 m_2 q_1 q_2 s_1 p_2} E^{(1)}_{p_Ls_Lkm_1m_2}E^{(2)}_{p_Rs_Rkm_1m_2}E^{(3)}_{p_Lp_R q_1 q_2}
E^{(4)}_{s_Ls_R m_1 m_2  q_2} E^{(5)}_{m_1 m_2 s_1 q_1 q_2 p_2 l} E^{(6)}_{s_Ls_R\alpha_s s_1 m_1 m_2} E^{(7)}_{p_Lp_R\alpha_p p_2 m_1 m_2}, \label{etildecomp}
\end{align}
 which reduces the number and dimension of the used indices, and thus made the numerical calculations more feasible.

The first tensor, $E^{(1)}$, in the product
\begin{align}
E^{(1)}_{p_Ls_Lkm_1m_2}=&\delta_{k,0}\delta_{p_L,s_L}\delta_{m_1,1}\delta_{m_2,1}+\delta_{k,-1}
+\delta_{k,2}(1-\delta_{m_1,2}\delta_{m_2,2})  
+\delta_{k,2}\delta_{m_1,2}  \delta_{m_2,2} \delta_{p_L,s_L}
+\delta_{k,1} (\delta_{m_2,1}-\delta_{m_2,2}\delta_{p_L,s_L}),
\end{align}
connects  the left part $p_L$ and $s_L$,  if no system-bath interaction in $\Delta T$ occurs ($k=0$), if two interactions at the left side occur ($k=2$ with $m_1=2$ and $m_2=2$) or if one interaction on the right side takes place ($k=1$ with $m_2=1$).
The second tensor $E^{(2)}$:
\begin{align}
E^{(2)}_{p_Rs_Rkm_1m_2}=&\delta_{k,0}\delta_{p_R,s_R}\delta_{m_1,1}\delta_{m_2,1}+\delta_{k,-1} 
+\delta_{k,2}(1-\delta_{m_1,1}\delta_{m_2,1}) 
&+\delta_{k,2}\delta_{m_1,1}  \delta_{m_2,1} \delta_{p_R,s_R} 
+\delta_{k,1} (\delta_{m_2,1}\delta_{p_R,s_R}+\delta_{m_2,2})
\end{align}
works in the same way, only that it acts on the right side instead of the left side, thus right and left (also encoded in $1$ and $2$ for $m_1$ and $m_2$) is exchanged and a sign difference caused by the commutator.
The term,
$E^{(3)}$, is only important for two system bath interaction per $\Delta T$
\begin{align}
E^{(3)}_{p_Lp_R q_1 q_2}=& (\delta_{k,0}+\delta_{k,-1}+\delta_{k,1}) \delta_{q_1,1}\delta_{q_2,1} 
+\delta_{k,2}(\delta_{m_1m_2}\delta_{q_1,q_2} 
+\delta_{m_1,1}\delta_{m_2,2}\delta_{q_1,p_L} 
+\delta_{m_1,2}\delta_{m_2,1}\delta_{q_1,p_R}),
\end{align
and only for $k=2$ is something  done.
It either connects the indices $q_1$ $q_2$ of the correlation function $G$, if the first and second system-bath interaction occurs both on the left or on the right side, or $q_1$ is connected either to the left or right side of $p$, these contributions are a result of the double commutator of the system bath coupling.
The next factor $E^{(4)}$ completes the work for the correlation function $G$:
\begin{align}
E^{(4)}_{s_Ls_R m_1 m_2  q_2}=& \big(\delta_{k,0}+\delta_{k,-1}+\delta_{k,1})\delta_{q_2,1}\delta_{m_1,1}\delta_{m_2,1}
-\delta_{m_1,2}\delta_{m_2,1} \delta_{q_2,s_L} +\delta_{m_1,m_2} ) 
\end{align}
and connects for $k=2$ with interaction on left and right the index $q_2$ either to the left or right side of $s$.
The next factor $E^{(5)}$ serves two purposes:
\begin{align}
E^{(5)}_{m_1 m_2 s_1 q_1 q_2 p_2 l}=&  ((\delta_{k,0}+\delta_{k,-1})\delta_{m_2,1} \delta_{l,1} +\delta_{k,1}\delta_{m_2,l+1}) 
-\delta_{k,2}  \delta_{l,1} G_{s_1,q_1,q_2,p_2,m_1,m_2}, \label{G_tensor2}
\end{align}
which contains the double integrated correlation function $G$ for the two interaction in $\Delta T$ (case $k=2$) and for $k=1$ it connects $l$ shifted by one to $m_2$ to pass it on to the $M$ tensors.
The indices $q_1$ and $q_2$ are passed to $G$ from the previous tensors, where as $s_1$ and $p_2$ will come from the proceeding tensors, starting with $E^{(6)}$:
\begin{align}
E^{(6)}_{s_Ls_R\alpha_s s_1 m_1 m_2}=& (\delta_{k,0}+ \delta_{k,-1} ) \delta_{\alpha_s,1} \delta_{m_1,1}\delta_{m_2,1}\delta_{s_1,1}
+\delta_{k,2} \delta_{\alpha_s,1} (\delta_{m_1,1}\delta_{s_1,s_L}+\delta_{m_1,2}\delta_{s_1,s_R} ) 
\nonumber\\ &
+\delta_{k,1}\delta_{m_1,1}(\delta_{m_2,1}\delta_{s_L,\alpha_s+1}\delta_{s_1,1} 
+\delta_{m_2,2}\delta_{s_R,\alpha_s+1}\delta_{s_1,1}),
\end{align}
where for $k=2$, $s_1$ is depending on $m_1$ either to the left or right side of $s$.
For $k=1$ the left or right side of $s$ (depending on $m_2$) is connected to $\alpha_s$ shifted by $1$, to connect it to a $\tilde{M}$ tensor for the second system-bath interaction.
The last factor $E^{(7)}$ works in a similar way:
\begin{align}
E^{(7)}_{p_Lp_R\alpha_p p_2 m_1 m_2}=&(\delta_{k,0}+ \delta_{k,-1} ) \delta_{\alpha_s,1} \delta_{m_1,1}\delta_{m_2,1}\delta_{s_1,1}
+\delta_{k,2} \delta_{\alpha_p,1} (\delta_{m_2,1}\delta_{p_2,p_L}+\delta_{m_2,2}\delta_{p_2,p_R} ) 
\nonumber\\ &
+\delta_{k,1}\delta_{m_1,1}(\delta_{m_2,1}\delta_{p_L,\alpha_p+1}\delta_{p_2,1}
+\delta_{m_2,2}\delta_{p_R,\alpha_p+1}\delta_{p_2,1}),
\end{align}
for $k=2$ it connects either the left or right side (depending on $m_2$) of $p$ to $p_2$ for passing it to the correlation function $g$. For $k=1$ the left or right side of $p$ is passed  to $\alpha_p$ shifted by $1$. 

\section{Green functions and model for the photonic crystal cavity system}

The electromagnetic response functions can be defined in terms of the photonic Green functions, defined from
\begin{equation}
\boldsymbol{\nabla}\times\boldsymbol{\nabla}\times\mathbf{G}(\mathbf{r},\mathbf{r}',\omega) -\frac{\omega^2}{c^2}\epsilon(\mathbf{r},\omega)\mathbf{G}(\mathbf{r},\mathbf{r}',\omega)=\mathbb{1}\delta(\mathbf{r}-\mathbf{r}')\label{eq: GFHelmholtz2},
\end{equation}
with suitable boundary conditions.
Following the approach of
Refs.~\onlinecite{Yao:09,PhysRevLett.98.083603}, these can be obtained
in an analytical form using mode expansion techniques and exploiting the completeness relation of the Green function, and we consider the two cavity modes and the waveguides modes as the only modes of interest in the problem. We assume identical cavities, that are side-coupled to a waveguide mode, and separated by a distance $L$;
we also consider a TLS embedded in each cavity at a field antinode position, aligned with the polarization of the
cavity mode. Subsequently, the 
 correlation functions,
 defined from
 $G_{nm}(\omega)=\mathbf{d}^{(n)}_{12}{\cdot}\mathrm{Im}({\bf G}(\mathbf{r}_1,\mathbf{r}_1,\omega))   \cdot \mathbf{d}^{*(m)}_{12}$, 
 take the form \cite{Yao:09}
\begin{align}
    &G_{11}(\omega)=\frac{\omega^2 E_1 E_1}{\omega_{c}^2-\omega^2-i \omega\gamma_{c}(1+e^{i2\phi}R_2(\omega))-i\omega\gamma_0}, \\
    &G_{22}(\omega)=\frac{\omega^2 E_2 E_2}{\omega_{c}^2-\omega^2-i \omega\gamma_{c}(1+e^{i2\phi}R_1(\omega))-i\omega\gamma_0}, \\
    &G_{12}(\omega) = \frac{\omega^2R_1(\omega) e^{i\phi}E_1 E_2}{\omega_c^2-\omega^2-i \omega \gamma_c(1+e^{i2\phi}R_1(\omega))-i \omega \gamma_0}, \\
    &G_{21}(\omega) = \frac{\omega^2R_2(\omega) e^{i\phi}E_1 E_2}{\omega_c^2-\omega^2-i \omega \gamma_c(1+e^{i2\phi}R_2(\omega))-i \omega \gamma_0}
\end{align}
with $E_{1/2}=\frac1{\epsilon_rV_{\rm eff}}$ (with $V_{\rm eff}$  the effective mode volume and $\epsilon_r$  the dielectric constant of the photonic crystal slab), $R_{1}(\omega)= R_2(\omega) = \frac{i\omega\gamma_c}{\omega_c^2-\omega^2-i\omega(\gamma_c+\gamma_0)}$, and $\gamma_c=Q_c/\omega_c$. The cavity decay rates include
on-chip decay from the cavity to the waveguide, $\gamma_c$, as well as off-chip
decay, $\gamma_0$, where the latter is typically much smaller than the former (to maintain a waveguide mode beta factor). The round trip phase
is determined from
$2\phi = 2k_BL + \varphi$, $k_B=\omega n_g/c$
and we set $\varphi=0.8 \pi$.



Parameters used in  the main text (if not others specified): 
$\epsilon_r=12$, $L=2 \lambda$ (waveguide length in wavelengths)
$n_g=10$ (group velocity index), 
$\omega_c=200\mathrm{THz}$ frequency of the cavities, $Q_c=2000$ quality factor of the cavity, $V_{\rm eff}=5\cdot 10^{7} {\rm nm}^3$ mode volume, $\varphi=0.8 \pi$ 
out of plane loss $1\%$ of cavity loss ($\gamma_0=0.01\gamma_c$).


\section{Dipole correlation function and localized detection}
\begin{figure}[bt]
	\includegraphics[width=0.6\textwidth]{figqnmspec_sub.png}
	\caption{Dipole correlation function $\langle P(t)  P^{\dagger}(0)\rangle $ for a system initially on an excited density at emitter 1 for far field detection (detection in $P(t)$ contains both emitters) and for localized detection at the quantum emitter (detection in $P(t)$ contains only one emitter).  }

  
	\label{fig_qnmspec_sub}
\end{figure}
The single quantum 
$\langle P(t)  P^{\dagger}(0)\rangle $ and double quantum  dipole correlation function $\langle P(t) P(t)P^{\dagger}(0) P^{\dagger}(0)\rangle$ allow a direct inspection of the system eigenenergies, and spectral resonances, for one or two excitations in the system, if the system is initially in the ground state. For example, $\langle P(t)  P^{\dagger}(0)\rangle $ describes a coherence between the ground state and single excited states created by polarisation flip operator $P^{\dagger}(0)$ at time $0$ and inspected at time $t$ by $P(t)$. Thus by varying $t$, the coherent oscillation can be scanned, and its Fourier transform corresponds to the linear optical response and we can deduce from the individual resonances the bright single excited states of system.


On the other hand, for $\langle P(t) P(t)P^{\dagger}(0) P^{\dagger}(0)\rangle$ at time $0$,  $P^{\dagger}(0) P^{\dagger}(0)$ described the coherence between the ground state and double excited states; similar to the linear, single excitation case, the double excited states can be scanned and a Fourier transform retrieves the energies of the  double excited states available for two photon excitation.
As in the linear case, experimental techniques exist to access this correlation function \cite{kim2009two}.

However, for the example system
in the main Letter,
the situation is a bit different to the usual experimental situation, since we include only polarization operators of the quantum emitters, for the operators $P^\dagger$. But in the system, the coupling to the photons can be  strong, so that eigenstates are in fact polariton states, that also include a photon mode part.
Assuming photon mode operators $b^\dagger_j$, then the single excitation eigenstates $|e\rangle$ are a superposition state,
\begin{align}
    |e\rangle=c^1_e \sigma_1^\dagger |0\rangle +c^2_e \sigma_2^\dagger |0\rangle+ \sum_j c^{pj}_e b^\dagger_j , 
\end{align}
with $\sigma_l^\dagger= |2\rangle_l \langle 1|_l$ and an overall ground state $|0\rangle$. So it is clear that we miss the photon part in the spectra.
In an analogues way, the double excited eigenstates $|f\rangle$ can be expanded in the local states.

In an experimental pump probe experiment, the pulses are also only acting on the quantum emitters and are much faster compared to the transfer process to the photon modes, thus that local densities $|2\rangle_l \langle 2|_l$ are prepared in the excited system.
We can access this contribution also by inspecting $\langle P(t)  P^{\dagger}(0)\rangle $ for a system where the first quantum emitter is initially in the local excited state. Of course this is not an eigenstate of the system but a superposition of $|e\rangle$, furthermore after the polarisation operator at $t$ the system is in a coherence between states  $|f\rangle$ and $|e\rangle$, creating the second excitation in the other emitter.
We see in Fig.~\ref{fig_qnmspec_sub}  precisely the same features as in the nonlinear signal in the main text.  The negative features are of course puzzling, but can be a result, that the system is in the local basis and thus expansion coefficients such as $c^1_e$ or $c^{12}_f$ etc, and their phase enter the signal. While our method does not allow to access the bath---the photon part---we can show, that if we detect only one quantum emitter \cite{localized_excitation,PhysRevB.86.085308,Schlosser_2013,specht2016reconstruction}, formally dark state can give negative contributions in the spectrum (cf.~Fig.~\ref{fig_qnmspec_sub}) and likely this is also the case for the negative contributions with respect to the photon part.

\section{Convergence analysis and consequences of approximate unitary}

\begin{figure}[hbt]
	\includegraphics[width=0.6\textwidth]{figqnmsuppconv.png}
	\caption{Reduced density matrix element $\rho_{2222}$ for the example from Fig. \ref{fig_qnmspec} d) for different discretizations $\Delta t$ of the quantum path integral. }

  
	\label{fig_qnmspec_sub_conv}
\end{figure}
Main parameters for the numerically convergence are the truncation bond dimensions and accuracy for the tensor network and the quantum path integral time discretization $\Delta T$.

For truncating bonds the MPS after MPO application, we use a maximum bond dimension as well as a accuracy for determining the truncated dimension. This is necessary since the different layers of the network within a time step have quite different bond dimensions (one order magnitude possible).
So for some of them bond dimension is the best parameter to do the truncation and for the others accuracy is the best for keeping the important information in the network. This is at least our subjective experience while designing the algorithm. Some intermediate layers require relative low bond dimension.  There relative useless information for the next will be included without an accuracy truncation parameter. This will prevent a fast calculation in subsequent layers, which require higher bond dimensions cause by the additional useless information.
Even worse, it could prevent convergence at all.

The time propagation  $U_i$ over an interval $\Delta T$ is expanded in Eq. \eqref{uexpansion} in orders of the system-bath interaction, thus it is only $O(\Delta t^2)$ unitary. In principle this can have consequences on semi-positivity and trace preservation of the reduced density matrix thus the density matrix elements over time as plotted in Fig. \ref{fig_qnmspec} d).

In fact, if we choose a $\Delta T$, which is too large,
we see large negative values e.g. for the reduced density of the double excited state matrix $\rho_{2222}$ as shown in Fig. \ref{fig_qnmspec_sub_conv}, but these problems disappear for sufficiently small $\Delta T$.
Problems with violation of density matrix positivity also occur if too strong truncation removes two many quantum path ways.
The violated positivity in $\rho_{2222}$ can also be seen in the numerical calculation used for Fig. \ref{fig_qnmspec}(a), but it occurs after the polarizations are created, so it can not be the origin of the negative signatures, since only secular terms are included in the system bath coupling, which do not convert polarizations into densities and vice versa. Therefore the densities are completely decoupled from the plotted observables.

\end{widetext}



 



\end{document}