Evolution of a Non-Hermitian Quantum Single-Molecule Junction at Constant Temperature

This work concerns the theoretical description of the quantum dynamics of molecular junctions with thermal fluctuations and probability losses. To this end, we propose a theory for describing non-Hermitian quantum systems embedded in constant-temperature environments. Along the lines discussed in [A. Sergi et al., Symmetry 10 518 (2018)], we adopt the operator-valued Wigner formulation of quantum mechanics (wherein the density matrix depends on the points of the Wigner phase space associated to the system) and derive a non-linear equation of motion. Moreover, we introduce a model for a non-Hermitian quantum single-molecule junction (nHQSMJ). In this model the leads are mapped to a tunneling two-level system, which is in turn coupled to a harmonic mode (i.e., the molecule). A decay operator acting on the two-level system describes phenomenologically probability losses. Finally, the temperature of the molecule is controlled by means of a Nosé-Hoover chain thermostat. A numerical study of the quantum dynamics of this toy model at different temperatures is reported. We find that the combined action of probability losses and thermal fluctuations assists quantum transport through the molecular junction. The possibility that the formalism here presented can be extended to treat both more quantum states (∼10) and many more classical modes or atomic particles (∼103−105) is highlighted.


Introduction
Molecular junctions are nano-devices composed of metal or semiconductor electrodes known as leads [1,2]. A single molecule simulates a conducting bridge between the leads. Such a nanostructure is almost perfectly suited to study non-equilibrium quantum transport [1,2].
Various approaches satisfactorily befit the investigation of the phenomenology of molecular junctions [2]. In this work, we introduce a simple toy model of a molecular junction [2] in order to perform a qualitative study aimed at singling out universal features of such systems. To this end, the operator-valued Wigner formulation of quantum mechanics [3][4][5] is particularly useful to embed quantum toy-models in quantum phase space baths since it can simultaneously describe dissipative phenomena and thermal fluctuations. A quantum-classical molecular junction model has been studied in [6].
In general, theories of quantum transport offer an acceptable representation of excitations and state population transfer. They usually assume that what is transferred does not disintegrate or disappear from the system. Such a physical property is expressed through the conservation of probability. On the contrary, by definition, the probability in quantum systems with gain and loss is not conserved. Such a feature can be incorporated in the theory by introducing non-Hermitian Hamiltonians [7][8][9][10][11][12]. Recently, the development of non-Hermitian quantum mechanics (NHQM) has been propelled by various experiments on many systems (see [13,14] and references therein). While in this work a non-equilibrium quantum statistical mechanics of molecular junctions is adopted, linear equations of motion are used in [10][11][12]. Conceptually, this can be understood in terms of the different theoretical targets: In [10][11][12] the target is to study non-equilibrium transport. Instead, our work is aimed at showing the possibility of a statistical mechanical description for non-Hermitian quantum systems in classical bath [15].
Various investigations on thermal effects in molecular junctions are found in the literature [26][27][28]. However, to the best of our knowledge, thermal effects have been investigated separately from probability losses of molecular junctions. In light of the previous considerations, the purpose of this work is to combine the non-linear equation for the density matrix of a system with a non-Hermitian Hamiltonian and temperature control in quantum phase space dynamics [29]. Temperature control is technically implemented by means of the Nosé-Hoover chain [30] thermostat, formulated in terms of quasi-Lie brackets [3]. The unified formalism here reported is developed for studying a greater variety of phenomena than those accessible to theories treating thermal fluctuations and probability losses separately. In detail, the theory is obtained through the density matrix approach to non-Hermitian quantum mechanics [15][16][17][18][19][20][21][22] and it is expressed in terms of the operator-valued formulation of quantum mechanics [3][4][5]15]. For clarity, we highlight the three mathematical ingredients that must be combined for building our formulation. First, we need an approach for describing quantum operators in terms of operator-valued Wigner phase space functions. The conceptual framework of this specific theory assumes an Eulerian picture where a Hilbert space, or a space of quantum discrete states, is defined for each point X = (Q, P) of phase space. Accordingly, the operator-valued Wigner function [3][4][5] may also be called phase-space-dependent density matrix. The consistent coupling of quantum and phase space degrees of freedom requires that quantum dynamics in the Hilbert space associated to X also depends non-locally on the quantum dynamics associated to different points X . Likewise, the dynamics of X does not only depend on its quantum space but it also depends on the quantum spaces of other points X and vice versa. In practice, such a dependence can be implemented through a variety of algorithms [31][32][33][34]. When the Hamiltonian is quadratic, one can expand the operatorvalued Moyal bracket [35,36] up to linear order inh and obtain quantum equations of motion [37].
The second ingredient entering our theory is the requirement that the degrees of freedom in phase space evolve in time satisfying the constraint determined by the constant thermodynamic temperature. The temperature control of the classical degrees of freedom is achieved in silico by means of the Nosé-Hoover chain algorithm [30]. We find such an approach advantageous because the formulation of the Nosé-Hoover chain algorithm can be realized within the theoretical framework given by the quasi-Lie brackets [3]. Such brackets allow one to generalize the Nosé-Hoover chain algorithm, originally formulated for classical systems only, to the more general case of the operator-valued formulation of quantum mechanics [29,35,36].
The third ingredient, of course, is that the modeling of the leads provides a non-Hermitian quantum system. When all these three theoretical features are simultaneously taken into account, we obtain an equation that generalizes the one arising from the quantum-classical quasi-Lie bracket [3] because it can describe probability gain/loss [15] too. Moreover, since such an equation keeps the thermodynamic temperature of the phase space degrees of freedom constant by means of the Nosé-Hoover chain algorithm [3,30,35,36], it generalizes the equations used in [15][16][17][18][19][20][21][22], which did not implement any thermodynamic constraint.
In this paper we investigate the time evolution of a single-molecule junction toy-model represented by a non-Hermitian Hamiltonian with a parametric phase space dependence. The isolated leads are represented by a two-level system. In this situation, the transport of the occupation probability of one or the other state can only occur through quantum tunneling. The molecule connecting the two leads is modeled by a harmonic oscillator. When we consider only the coupling between the harmonic mode and the two-level system, the total system is considered closed. In this case, transport can take place both via tunneling and via the channel given by the oscillator. In order to build a complete non-Hermitian quantum single-molecule junction model at constant temperature, the closed spin-boson system is transformed into an open quantum system by means of two theoretical steps. The first step is to consider a decay operator acting on the two-level system. This operator introduces another transport channel that does not conserve probability and then breaks time-reversibility.
The second step is to put the harmonic mode into contact with a heat bath. Of course, this condition forces the molecule to have the same temperature of the bath and to experience fluctuations according to the canonical ensemble probability distribution. In practice, the action of the constant-temperature bath on the molecule is implemented by means of the Nosé-Hoover chain thermostat. For this reason, in the following we will use wording such as "heat bath", "constant-temperature environment", or "Nosé-Hoover chain thermostat" interchangeably.
The final complete non-Hermitian Hamiltonian is designed to represent two lossy leads interacting with a thermalized molecule.
The organization of the paper is the following: In Section 2 we give a short summary of the density matrix approach to the dynamics of a system governed by a non-Hermitian Hamiltonian. In Section 3 we describe briefly how to treat non-Hermitian quantum systems interacting with classical degrees of freedom. We present the main theoretical results of this work in Section 4, where we introduce the mathematical formalism describing the quantum dynamics of a non-Hermitian quantum system interacting with classical degrees of freedom at constant temperature. In Section 5 we introduce the non-Hermitian quantum single-molecule junction model. The results of the numerical calculations are reported in Section 6. Finally, our conclusions are given in Section 7.

Quantum Dynamics of Non-Hermitian Systems
Consider a quantum system described by a non-Hermitian Hamiltonian operator: whereĤ andΓ are Hermitian operators (Γ is the decay operator). In the presence of a probability sink or source, the dynamics of the system is defined in terms of the following equations (we seth = 1): where |Ψ(t) and Ψ(t)| are state vectors of the system. Equations (2) and (3) reduce to the standard Schrödinger equation whenΓ = 0. The density matrix approach to non-Hermitian quantum systems [15][16][17][18][19][20][21][22] is obtained upon introducing at any t > 0 a Hermitian, semipositive defined, but non-normalized density matrix: The equation of motion forΩ [16,20,22] is obtained by taking the time derivative of Equation (4) and using Equation (1): where [..., ...] is the commutator and [..., ...] + is the anticommutator. Equation (5) is linear but, nevertheless, it describes an open quantum system. This is possible since it is defined in terms of a decay operator that physically describes the presence of probability sources/sinks, arising from the hidden system with which the open quantum system interacts. Equation (5) is non-Hermitian, breaks time reversibility and does not conserve the trace ofΩ(t). In fact, from Equation (5) one can easily derive: Tr Ω (t) = −2Tr ΓΩ (t) .
In order to obtain a well-founded statistical mechanics of non-Hermitian quantum systems, one can introduce the normalized density matrix: Tr Ω (t) .
This operatorρ(t) obeys the non-linear equation: As a consequence of its very definition and of Equation (8), Tr(ρ(t)) = 1 at all times. Equation (8) breaks time reversal symmetry and is non-linear. In practice, the formalism trades off the non-conservation of probability in Equation (5) with the non-linearity of Equation (8). As the linear Equation (5), the non-linear Equation (8) describes the dynamics of an open quantum system. Consistently, averages of arbitrary operators, e.g.,χ, can be calculated by means of the normalized density matrixρ as: Equation (9) establishes the foundations of the statistical interpretation of non-Hermitian quantum mechanics.

Non-Hermitian Quantum Systems Set in Phase Space
Non-Hermitian quantum systems can be embedded in phase space [15]. In the following, we sketch the derivation given in [15]. Let us consider a system described by a hybrid set of coordinates (x, X), wherex are quantum coordinates (operators) and X = (Q, P) are phase space coordinates (an obvious multidimensional notation is used in order not to clutter formulas with too many indices). In this specific case, the Xs describe the phase space coordinates of the system. As we have already discussed in the Introduction, the description of the system, using the operator-valued-formulation of quantum mechanics, can be established via an Eulerian point of view according to which a space of state vectors is defined for any point X of phase space [3][4][5] so thatx (and operators defined in terms ofx) can act on it. In such a hybrid space there are operators that depend only onx, which we denote with the symbolˆon top of the operator, e.g.,Γ. There are pure phase space functions which we denote by their arguments, e.g., H B (X), or G(X, t) when there is also an explicit time dependence. Finally, there are operators that depend both onx and (parametrically) on X. We denote the latters with the symbol˜on top, e.g.,Ω(t) ≡Ω(X, t).
In the non-Hermitian case, we consider a Hamiltonian operator with parametric dependence on phase space:H whereĤ is the Hermitian Hamiltonian of the quantum subsystem, H B (X) is the Hamiltonian of the degrees of freedom in phase space,H I describes the interaction between the quantum system and the bath, andH S =Ĥ + H B (X) +H I . Of course,H is the total non-Hermitian Hamiltonian of the system withΓ as the decay operator. When the operatorΓ acts on the quantum dynamical variables of the non-Hermitian quantum system, it has been shown [15] that, as a generalization of the equations given in [3][4][5], the equation of motion becomes: whereΩ(t) is the phase-space-dependent non-normalized density matrix of the system and is the symplectic matrix written in block form. Using it, one can write the Poisson bracket in the form adopted in Equation (18): The dimension of phase space is 2N. The trace of the phase-space-dependent density matrixTr is defined as: Tr where Tr denotes a partial trace over the quantum dynamical variablesx and dX is the integral in phase space. The traceTr Ω (t) obeys the equation of motion: The result given in Equation (15) is proven in Appendix A.
Let us now introduce a phase-space-dependent normalized density matrixρ(X, t): Quantum statistical averages of an arbitrary phase-space-dependent operatorχ can now be calculated as: The phase-space-dependent normalized density matrix obeys the non-linear equation of the motion below:

Non-Hermitian Quantum Systems Set in a Constant-Temperature Phase Space
Equations (11) and (18) describe the dynamics of non-Hermitian quantum systems embedded in phase space [15]. When considering open quantum systems in more general settings, it is common to consider the effects of thermal fluctuations. Phase space coordinates undergo thermal fluctuations when their microscopic equilibrium state is described by the canonical distribution function. In this ensemble the thermodynamic temperature is constant. Within the framework of the operator-valued Wigner formulation of quantum mechanics [3][4][5], temperature constraints are efficiently formulated by means of quasi-Lie brackets [3,35,36]. In the classical case, deterministic thermostats, such as the Nosé-Hoover chain thermostat [30], are also implemented by means of quasi-Lie brackets [38][39][40]. Below, we briefly explain how this is achieved. In silico (i.e., on the computer), temperature control is realized by augmenting the dimensions of phase space using just a few additional degrees of freedom. The new higher dimensional phase space is known as extended phase space. We indicate points in the extended phase space by means of the symbol X e . In the case of the Nosé-Hoover chain thermostat, upon introducing two fictitious coordinates Λ 1 , Λ 2 with their associated momenta Π 1 , Π 2 , we have the following point in the extended phase space: To lighten the notation, the phase space coordinates of the Nosé-Hoover chain thermostat are denoted as Y = (Λ 1 , Λ 2 , Π 1 , Π 2 ). For the reader, the symbol Y denotes a calligraphic Y. A generic operator acting upon the extended phase space is identifiable by the apex e. Considering a classical system with a potential term V(Q) and Hamiltonian H B (X) = P 2 /2 + V(Q), the extended phase space Hamiltonian (see [38][39][40] and references therein) includes the Nosé-Hoover chain energy [30,38] and is defined as: where µ K are fictitious inertial parameters, T is the termodynamic temperature, k B is the Boltzmann constant, and N is the number of physical coordinates Q that are thermalized. The equation of motion in extended phase space can be written as: where the first and second equality define a quasi-Lie bracket. The last equality enlighten the matrix structure of the equations of motion. This last form also allows one to compactly define the compressibility of phase space: The antisymmetric matrix B e , entering Equations (21) and (22), is defined as: is an antisymmetric phase-space-dependent matrix generalizing the symplectic matrix in Equation (12). Equation (21) is called a quasi-Hamiltonian equation because, while conserving H e (X e ), it cannot be derived from the Hamiltonian formalism alone: In order to write the equation of motion (21), together with H e (X e ), one also needs the antisymmetric matrix B e in Equation (23). The equation of motion for the distribution function f (X e , t) is: Given the discussion above, our task becomes that of generalizing Equation (24) to a non-Hermitian quantum system. The phase-space-dependent Hamiltonian of the quantum system must be defined on the extended phase space of the Nosé-Hoover chain thermostat: A phase-space-dependent non-normalized density matrixΩ e (t) can be introduced as well. We postulate thatΩ e (t) obeys the equation of motion: (26) In the absence of the quantum degrees of freedom, Equation (26) reduces correctly to Equation (24). Moreover, when (Λ I , Π I ) → 0, with I = 1, 2, (which means that the Nosé-Hoover thermostat is not applied), Equation (26) reduces unerringly to Equation (11).
The traceTr e of the density matrixΩ e (t) involves an integral over the extended phase space:T r e Ω e (t) = Tr dX eΩe (X e , t) .
The equation of motion of the trace defined in Equation (27) has the same structure of Equation (15): Appendix B shows how Equation (28) is obtained.
The phase-space-dependent normalized density matrix in the extended phase space, ρ e (X e , t), is defined in analogy with Equation (16) as: Tr e Ωe (X e , t) .
In analogy with Equation (17), quantum statistical averages are calculated as: The phase-space-dependent normalized density matrix in Equation (29) obeys the nonlinear equation of motion: Tr e Γρ e (t) We remark that Equations (26) and (31) are two important results of this work since these equations define the dynamics of a non-Hermitian quantum system embedded in a classical bath at constant temperature.

Model of a Non-Hermitian Quantum Single-Molecule Junction at Constant Temperature
In this section we formulate a molecular junction toy-model described by a non-Hermitian Hamiltonian denoted byH m . The various energy contributions entering the Hermitian terms in part of the Hamiltonian model are: The two-level system HamiltonianĤ m describes quantum transport between the leads in terms of the variation of the population of a ground and an excited state. The extended phase space point X m is defined as in Equation (21) but the model considers a single Q coordinate and its conjugate momentum Π. Transport in the isolated two-level system takes place through tunneling. The inclusion in the model description of the decay operator Γ m opens a channel through which the population of the states can disappear forever from the two-level system. The explicit form ofΓ m will be given later on when we transform to the adiabatic basis. This subsystem is also coupled to a harmonic mode, with free Hamiltonian H m B (Q), by means of an interaction HamiltonianH m I . The harmonic mode opens another channel for the transport of the population between the leads. In addition to this, the harmonic mode is embedded in a thermal bath via the Nosé-Hoover chain thermostat with energy H m nhc . The complete non-Hermitian Hamiltonian modelH m reads: where the last equality also defines the Hermitian part of the Hamiltonian model, i.e., To the authors' knowledge, the Hamiltonian model in Equation (36) provides a novel non-linear approach to the modeling of lossy molecular junctions in thermal baths. Figure 1 displays a pictorial representation of the non-Hermitian quantum single-molecule junction (nHQSMJ) model at constant temperature. As the TLS moves toward the right lead, the arrow rotates until it reaches the upward position, which implies a complete transfer to the right lead. The sink below the right lead absorbs the TLS probability in an irreversible way. In the formalism, the sink is represented by a Hermitian decay operator. The decay operator acts on quantum states that in turn determine probabilities. Hence, the action of the sink unfolds upon an ensemble of trajectories of the TLS system. The picture also shows a box separating the molecular junction from the environment. The drawing of the box resembles on purpose that of an oven with a thermostat on top. The thermostat temperature is set at the same value of the T R temperature of the environment. The thermostat controls the temperature inside the oven so that it is equal to that of the environment. The yellow color inside the oven and around it, in the environment, conveys the idea of thermal equilibrium.
The abstract dynamics of the nHQSMJ model is obtained upon using the non-Hermitian Hamiltonian (36) in Equation (31). The abstract equations of motion can be represented using different basis sets. The approach described in [31] is based on the representation in the adiabatic basis. We also adopt such an approach here. The adiabatic Hamiltonian is defined as:h and the adiabatic basis is introduced through the eigenvalue equation: In practice, the adiabatic states are a kind of dressed states of the quantum system when it interacts with approximately "frozen" classical degrees of freedom. The operatorΓ m is chosen in a phenomenological way. In the adiabatic basis, it is taken as: In this way we describe a model where the excited state is in contact with a sink while the ground state is left unperturbed. The abstract equation of motion for the non-normalized density matrix of the model is: The antisymmetric matrix B m has the same structure of that in Equation (23), but it considers only the physical degrees of freedom of the harmonic modes. Accordingly, the phase space compressibility of the model becomes κ m S = −(Π 1 /µ 1 ) − (Π 2 /µ 2 ). Equation (41) introduces the two Liouville operators: The equation of motion (41) is represented into the adiabatic basis as: In Equation (44) we introduced the two Liouville super-operators −iL m αα ,ββ and −iL Γ m αα ,ββ . We first considerL m αα ,ββ . We have: where ω αα = E α (Q) − E α (Q) is the Bohr frequency. The operator: is a classical-like Liouville operator generating the dynamics of the physical coordinates X under the feedback of the fictitious coordinates Y. In Equation (46) we have also introduced the Hellmann-Feynman force: This is the force acting on the harmonic oscillator when state α is occupied.

Calculations and Results
We consider a situation where the molecule is weakly coupled to the leads. This implies that neglecting the action of the transition operatorT αα ,ββ in Equation (51), the phase-space-dependent density matrix is propagated by means of the sequential short-time propagation (SSTP) algorithm [31,32] as: where τ = t/n step and the non-Hermitian adiabatic Liouville super-operator is −iL Quantum statistical averages can be calculated as: where the bracket ... m stands for the average in the extended phase space of the model. We implement Equation (53) in the following way [31,32]. Once the quantum initial state is assigned for every point of the extended phase space, the calculation of quantum averages can be performed by sampling the initial X e with a Monte Carlo algorithm. The point X m is then propagated over the assigned quantum state for a time length t. The time step τ = t/n step in Equation (53) must be chosen and small enough to minimize the numerical error. In the calculation we take the following uncorrelated form of the initial phase-space-dependent density matrix: whereΩ Calculations were performed with fixed values of γ = 0.1, ω = 1/3, ∆ = 1, and c = 0.007, time step τ = 0.005, number of time steps n step = 10 4 , and number of Monte Carlo steps n mcs = 25 · 10 3 . All the values of the parameters are given in a dimensional units. A total of 20 different values of the inverse temperature were considered upon chosing β l+1 = 1/T l+1 = β 0 (1 + l), with l = 0, ..., 19 and β 0 = 0.0005. For each β l we have investigated four different types of dynamics: (i) Unitary dynamics, (ii) constanttemperature dynamics, (iii) non-unitary dynamics, and (iv) non-unitary dynamics at constant temperature. Initial conditions are chosen with the two-level system in its ground state and the bath in a thermal state. It is useful to introduceΞ m andX m , the nonnormalized and normalized reduced density matrices, respectively. We write their matrix elements in the adiabatic basis as: Both reduced density matrices can be easily calculated by means of our numerical approach. In order to check the numerical scheme, we calculated the evolution in time of Tr(Ω m (t)). The results are shown in Figure 2. When the decay operator is set to zero, Panels (a) and (b), the trace of the density matrix is conserved with extremely good numerical precision both when the dynamics is unitary and when the temperature is controlled. This provides a convincing indication that our algorithm conserves important dynam-ical invariants. Panels (c) and (d) in Figure 2 show the decrease in time ofTr(Ω m (t)) non-unitary dynamics. Moreover, there is no appreciable difference between unitary and constant-temperature dynamics. In Figure 3 we show the time-evolution of Ξ m αα (t) and X m αα (t) m , with α = 1, 2, defined in Equations (59) and (60), respectively. Panels (a) and (b) display Ξ m 11 (t) and Ξ m 22 (t), respectively. Initially, adiabatic states are equally occupied. However, according to Equation (40), only the occupation of the excited state experiences depletion, see in Panel (a). Panels (c) and (d) of Figure 3 show X m 11 (t) and X m 22 (t), respectively. The behavior of X m 11 (t) is hardly distinguishable from that of Ξ m 11 (t). However, the situation is different for X m 22 (t) in Panel (d). As a matter of fact, while Ξ m 22 (t) is constant, X m 22 (t) increases so that Tr m (ρ m (t)) = X m 11 (t) + X m 22 (t) = 1.   Figure 2. There is no difference between Panel (a,c) since the trace decreases because of the decrease of Ω 11 . The comparison between Panel (b,d) and further clarifies that the cause of the population loss of the excited state is the action of the decay operator and not the energy transfer.
In Figure 4 we show the real part of the reduced off-diagonal element Re(X m 12 (t)).    One at a small frequency and one at a higher frequency. The Rabi-like oscillations of the population difference at a distant time are signs of the absence of stable quantum transport in the nHQSMJ model. Nosé-Hoover chain dynamics, shown in Panel (b), suppresses the peak at small frequency. The Rabi-like oscillations of the population difference at distant time are less pronounced: They take place at a relatively high frequency, as the Fourier transform shows, around zero. This indicates that in the presence of thermal noise the population difference oscillates around zero so that transport is stable. Figure 6 shows the time evolution of the population difference β = 0.0050, which corresponds to a temperature higher than that in Figure 5. All the effects shown by Figure 5 appear more marked in Figure 6. This can be seen as a particular instance of environmentassisted quantum transport: Remarkably, the noise provided by an environment can defeat Anderson localization [41] and assist quantum transport. The idea that the noise of the environment could enhance quantum transport, instead of suppressing it, has been postulated to explain the high efficiency of energy transport in photosynthetic systems [42].
Obviously, the noise-enhancement of quantum transport can take place only below a certain threshold. Above it, transport is suppressed by the quantum Zeno effect [43]. In Panel (c), Hamiltonian dynamics of the bath is combined with a non-zero decay operator acting on the two-level system. As seen, the decay operator suppresses the high frequencies in the Rabi-like oscillations of the population difference at distant time. Such oscillations show that, even in this case, quantum transport is not very stable. Comparison of the dynamics in Panel (b) and in Panel (c) shows that there is a difference between thermal and probability sink dissipation: The first suppresses low frequency oscillation while the second damps high frequency oscillations. As one can expect, the combination of the two dynamics, shown in Panel (d), leads to the suppression of both frequencies, increasing the efficiency and the stability of the population transport. Calculations for β = 0.0025 are reported in Figure 7. The results are indistinguishable to the human eye from those obtained at β = 0.0075. This reinforces the conclusion that quantum transport in the nHQSMJ model is enhanced by coupling the two-level system with a probability sink and controlling the temperature of the molecule.    it is shown that a higher temperature damps high frequency oscillations in SMJ+NHC dynamics. Statistical error appears to be somewhat greater. Panel (c) shows that the effect of the decay operator is to damp high frequency oscillations in SMJ+nH dynamics. Panel (d) shows that the combined action of the NHC thermostat and the decay operator, i.e., SMJ+nH+NMJ dynamics determines a stable transfer process. Both low and high frequency oscillations are strongly damped, except for a very weak oscillation at a very low frequency.  . Panel (b) shows that the thermostat (SMJ+NHC dynamics) facilitates the population transfer and stabilizes it more than in the case of lower temperatures. The low frequency terms are highly damped while the high frequency ones are unaltered. Panel (c) shows that the effect of the decay operator (case of the SMJ+nH dynamics) is to cancel the high frequency term. Finally, Panel (d) shows that the SMJ+nH+NHC dynamics displays a stable transfer: Both the low and the high frequency terms are strongly damped.

Conclusions
In this work, we constructed a theory for studying non-Hermitian phase-spacedependent quantum systems at constant temperature. This theory is based on an operatorvalued Wigner formulation of quantum mechanics or, in other words, on a phase-spacedependent density matrix. The condition of constant temperature in phase space was achieved by means of the Nosé-Hoover chain thermostat. A mathematical result of the formalism is the derivation of the non-linear equation of motion for the normalized phasespace-dependent density matrix of the system. We remark that this result, considering a constant temperature bath, improved the theory developed in [15] where the temperature fluctuations were not constrained. These latter situation is somewhat unrealistic in common experiments.
This theory was applied to a model of a non-Hermitian quantum single-molecule junction. We emphasize that this model treated probability loss and thermal fluctuations on the same level. In detail, the model comprised of a two-level system with a probability sink coupled to a thermalized harmonic mode. The structure of our model was conjectured upon drawing an analogy with the process of noise-assisted transport [44][45][46]. In our case, we expect that the assisted quantum transport arises from the combined action of the probability sink and the constant-temperature fluctuating molecule. Indeed, we observed a transport enhancement in silico upon simulating numerically the non-Hermitian quantum single-molecule junction model for a range of temperatures, while keeping the values of the other parameters constant. The Fourier transformed signal in frequency space showed that the Nosé-Hoover chain thermostat suppresses the slow frequencies while the probability sink damps the high frequency oscillations of the population difference of the non-Hermitian quantum single-molecule junction model at a constant temperature.
The non-Hermitian quantum single-molecule junction model at constant temperature introduced in this paper could more accurately be classified as a particular instance of a class of models. As a matter of fact, it could be easily generalized by considering a greater number of quantum states (∼10) and an even greater number of classical modes (∼10 3 − 10 5 ). Of course, more sophisticated algorithms would have to be used for calculating the evolution of the density matrix. We also note that spin models that were similar but simpler than the non-Hermitian quantum single-molecule junction model introduced in this paper could be solved analytically when obeying certain symmetry properties [47][48][49][50][51][52][53][54][55]. It would be interesting to investigate whether the extension of these models to non-Hermitian quantum mechanics would still be analytically treatable. We defer such generalizations to future work.