Figure 4. An example snippet to run the HEOM calculations within a Jupyter notebook.
The user needs to define the parameters of the bath, simulation, output, and choose the approximation control options. The user can also setup a desired type of electron-phonon coupling (system-bath correlation). In this example, we define the coupling matrices (representing the operators Eq. 2) to reflect one of the simplest scenarios – the coupling of a single bath mode to a single electronic state, such that\(Q_{m}=\left.\ |m\right\rangle\left\langle m|\right.\ \),\(m=1,2\). The zero’s mode is not coupled to any of the electronic states, \(Q_{0}=0\). The \(Q_{m}\) matrices define the type of correlation of bath and quantum degrees of freedom and may be a subject of particular investigation. For instance, one may want to couple the bath’s phonon \(0\) to the lowest two states at the same time, in which case one could define\(Q_{0}=\left.\ |0\right\rangle\left\langle 0|\right.\ +\left.\ |1\right\rangle\left\langle 1|\right.\ \). One can also anti-correlate the response of states 0 and 1 to the same phonon, say 1, in which case one would define\(Q_{1}=\left.\ |0\right\rangle\left\langle 0|\right.\ -\left.\ |1\right\rangle\left\langle 1|\right.\ \).
The parameters dictionary “params” in the above snippet defines the critical parameters affecting the whole HEOM set – the number “KK” that defines the number of Matsubara terms (KK + 1), the highest tier level of the hierarchy (“LL”), as well as the properties of the bath – reorganization energy (“eta”), system-environment interaction rate, which can be associated with the decoherence rate13(“gamma”), and bath temperature (“temperature”).
The snippet also illustrates the concept of the on-demand print-out/storage of the results. With the variable “mem_output_level” set to 3, all possible properties (so far “timestep”, “time”, and “denmat”) could be saved in the output files. However, one needs to explicitly include these keywords in the list of “properties_to_save”. The output level variable determines the outputs of which dimensionality/complexity can be saved. The value of 0 would normally save the scalar data, the value of 1 would save everything up to 1 D arrays (e.g. total energy vs. time), the value of 3 would save a 3D data such as a series of density matrices (each is a 2D entry) as a function of time (adding another dimension). However, one may have many irrelevant properties (as far as the goals of a particular simulation are concerned) one doesn’t want to save. The “properties_to_save” variable allows one to exclude saving such data by not listing the corresponding data names in it. Note also that listing the properties to save is required for them to be saved. Thus, just setting the “mem_output_level” parameter to 3 will not save anything if the “properties_to_save” list is empty.
A few other variables of note are the “do_scale” which turns on/off the scaled HEOM, Eq. 13, the “truncation_scheme”, which chooses one of four truncation schemes, Eqs. 5d, 13d, the “adm_tolerance”, which sets the threshold for discarding the ADMs,\(\left|{\tilde{\rho}}_{\mathbf{n}}\right|<tol_{1}\), and the “adm_deriv_tolerance”, which sets the threshold for skipping the integration of some HEOM equations,\(\left|\frac{d}{\text{dt}}{\tilde{\rho}}_{\mathbf{n}}\right|<tol_{2}\). With the variable “filter_after_steps” one can control the frequency of the updates of the active HEOM equations to be integrated as well as the updates of the “no-zero” ADMs to be used in computing time-derivatives. Ideally, this update can be done at every step, but one can gain some computational savings when increasing this parameter. It should be noted that the convergence with respect to these parameters should always be tested to ensure the meaningfulness of the computed results.