## Appendix D Translating Master Equations into Matrix Differential Linear Equations  
We will develop the translation from the HEOM master Equation (28) into the superoperator version through the application of the rule (29). For other quantum master equations, we simply define the supervector \(\vec{\rho}=(\rho_{11},\rho_{12},...,\rho_{1N},\rho_{21},...,\rho_{NN})\) (note other versions can be advised varying the order of the elements of \(\rho\)) and we go directly on the application of (29). In the HEOM method, due to the coupling among \(\rho\) and the auxiliary matrices \(\rho_{\mathbf{n}}\) through \(\rho_{\mathbf{n}_{k\pm}}\), we proceed similarly by moving each \(\rho_{\mathbf{n}}\) into its corresponding \(\vec{\rho}_{\mathbf{n}}\). Then, the HEOM is a system of linear differential equations for the supervector set \(\Omega\equiv\{\vec{\rho}_{\mathbf{n}}\}\) labeled with the vectors \(\mathbf{n}\) of deep \(D\) as was defined in Section 3.3.3, each one of size \(N^{2}\).  
Then, the application of the rule (29) transforms each term in (28) into the matrix product between a matrix of size \(N^{2}\times N^{2}\) (superoperators) and each one of the supervectors in \(\Omega\):\[[H_{\rm S}+\sum_{k=1}^{N}\lambda_{k}|k\rangle\langle k|_{\rho_{\bf n }}] \rightarrow \hat{n}\sum_{kJ}\kappa_{{\bf n}}^{(H,1)}{}_{ij,kl}\quad\rho_{{\bf n }_{kl}} \tag{111}\] \[\mbox{with}: \hat{n}\kappa_{{\bf n}}^{(H,1)}=(H_{\rm S}+\sum_{k=1}^{N}\lambda_ {k}|k\rangle\langle k|)\otimes{\bf 1}-{\bf 1}\otimes(H_{\rm S}+\sum_{k=1}^{N} \lambda_{k}|k\rangle\langle k|)^{T}\] \[-\sum_{k=1}^{N}n_{\rm F}\gamma_{k}\rho_{{\bf n}} \rightarrow \sum_{kJ}\kappa_{{\bf n}}^{(H,2)}{}_{ij,kl}\quad\rho_{{\bf n}_{kl}}\] (112) \[\mbox{with}: \kappa_{{\bf n}}^{(H,2)}=-\sum_{k=1}^{N}n_{\rm F}\gamma_{k}{\bf 1 }\otimes{\bf 1}\] \[(\{|3\rangle\langle 3|_{\rho_{\bf n}}\}+\{|4\rangle\langle 4|_{ \rho_{\bf n}}\}) \rightarrow -\frac{1}{r_{\rm trap}}\sum_{kJ}\kappa_{{\bf n}}^{(H,3)}{}_{ij,kl} \quad\rho_{{\bf n}_{kl}}\] (113) \[\mbox{with}: -\frac{1}{r_{\rm trap}}\kappa_{{\bf n}}^{(H,3)}=(|3\rangle\langle 3| +|4\rangle\langle 4|)\otimes{\bf 1}+{\bf 1}\otimes(|3\rangle\langle 3|+|4 \rangle\langle 4|)\] \[i\sum_{s=1}^{N}\lambda_{s}[V_{s,r}\rho_{{\bf n}_{s+}}] \rightarrow \sum_{kJ}\sum_{s}^{N}\kappa_{{\bf n}}^{(H,+)}{}_{ij,kl}\quad\rho_ {{\bf n}_{s+},kl}\] (114) \[\mbox{with}: \kappa_{{\bf n}}^{(H,s)}=i\Delta_{\rm S}(V_{s}\otimes{\bf 1}-{ \bf 1}\otimes V_{s})\] \[i\sum_{s=1}^{N}\alpha_{s}\alpha_{s}\rho_{{\bf n}_{s-}} \rightarrow \sum_{kJ}\sum_{s}^{N}\kappa_{{\bf n}}^{(H,s-)}{}_{ij,kl}\quad\rho_ {{\bf n}_{s-},kl}\] (115) \[\mbox{with}: \kappa_{{\bf n}}^{(H,s-)}=i\sum_{s=1}^{N}\Delta_{\rm S}n_{\rm S} ((V_{s}\otimes{\bf 1}-{\bf 1}\otimes V_{s})-i\theta_{s}(V_{s}\otimes{\bf 1}+{\bf 1} \otimes V_{s}))\] \[\Theta_{s}=[V_{s},\cdot]-i\theta_{s}\{V_{s},\cdot\},\quad\theta_ {s}=\frac{\beta\hbar\gamma_{s}}{2}\]  
where \((A\otimes B)_{ij,kl}=A_{ik}B_{jl}\). Finally, for all \({\bf n}\in\Omega\), it becomes:  
\[\dot{\vec{\rho}}_{\bf n} = \sum_{i=1}^{3}\kappa_{{\bf n}}^{(H,i)}\vec{\rho}_{\bf n}+\sum_{s=1 }^{N}\kappa_{{\bf n}}^{(H,s+)}\vec{\rho}_{{\bf n}_{s+}}+\sum_{s=1}^{N}\kappa_ {{\bf n}}^{(H,s-)}\vec{\rho}_{{\bf n}_{s-}} \tag{116}\]  
This task could appear complex but it is computationally convenient in a first instance due \(\kappa_{{\bf n}}^{(H,s\pm)},\kappa_{{\bf n}}^{(H,i)}\), for \(s=1,...,N,i=1,2,3\) are constant matrices. The numerical easiest process to integrate the equation is to consider simply that \(\rho_{\bf n}(t+\delta t)\approx\rho_{\bf n}(t)+\rho_{\bf n}(t)\delta t\) (better integration methods could be implemented to avoid tiny values of \(\delta t\)). The main problem with this approach to HEOM is the polynomial growing number of equations, for \(D=3\) and arbitrary \(N\), there are \(1+N+\frac{N(N+1)}{2}+\frac{N(N+1)(N+2)}{6}\) equations. They are 120 for \(N=7\) and 165 for \(N=8\) involving matrices of size \(49\times 49\) and \(64\times 64\), respectively, on a fine grained time period of tens of picoseconds.