Quantum Exact Response Theory Based on the Dissipation Function

1. Introduction

Statistical physics provides a successful description of systems in thermodynamic equilibrium through probability distributions known as ensembles. However, most systems of interest are not in equilibrium. Response theory studies their behavior when subject to external actions, usually adopting a perturbative approach that received a substantial impulse in the 1950s with the works of Callen, Green, Kubo and others [1,2,3,4,5,6]. The perturbation is assumed to be small compared to the Hamiltonian, and can be time dependent. In classical mechanics, paralleled by quantum mechanics, the initial unperturbed state of the system is represented by a probability distribution in phase space, that is postulated to be affected by the external action like a fluid in real space is affected by external forces, but preserves its mass. Therefore, a continuity equation for the probability, known as Liouville equation, is taken to hold in phase space [7]. If the external action is small, one assumes that the effect on the probability density can be approximated to the first order in the magnitude of the action, using the time correlation function of the perturbation and the observable of interest, computed with respect to the equilibrium probability. Kubo provided a single formalism to treat both classical and quantum dynamics, concerning probability densities in the first case and density matrices in the second, that, in the linear approximation, yields susceptibilities for systems not far from the thermodynamic equilibrium. The need to develop suitable techniques in the quantum branch of this field is ever more important [8,9,10]. For instance, the phenomenon of decoherence, typical of open quantum systems, represents a major challenge in quantum computing [11,12]. A deeper understanding of quantum dissipative dynamics could lead to significant advances in several fields. In this regard, Kubo’s theory has recently been generalized to non-equilibrium situations described by general time-local master equations [13,14,15], such as Lindblad equations, which are the most widely used for describing open quantum dynamics. Thermal equilibrium states are replaced by non equilibrium steady states, and dissipative perturbations are considered in addition to the Hamiltonian ones [16]. Extensions of linear response theory have been developed to also account for non-Markovian effects [17]. Different formulations of the response function have enriched the theory [18]. While these formulations produce equivalent results in numerous applications, the existence of different types of response functions provides important theoretical and experimental advantages, allowing one to choose the best one depending on the specific application.

As a result, this perturbative approach has achieved broad applicability and a considerable degree of completeness [19,20]. However, nonlinearities are common [21,22], especially when dissipation, phase transitions, and decoherence phenomena occur. In these situations, even a small perturbation can lead to significant modifications of the state, impairing the applicability of linear response. One possible approach to these issues consists in taking higher order terms in the perturbation expansions, which may provide greater precision, but it is often cumbersome and expensive, and may still fail to capture phenomena like phase transitions and anomalous behaviors. More generally, the perturbations may not be small. In all such cases, the best solution would be an exact theory of response. In recent years, research aiming at this goal has started to emerge. One interesting approach is the use of stochastic methods from classical physics to study quantum dynamics [23,24]. The stochastic reset technique is particularly suitable for representing measurement processes and, in this context, for describing open quantum systems exactly [25]. While promising, this method has the disadvantage of being applicable only to a narrow range of quantum systems, and it is not always possible to move to "classical" stochastic dynamics.

In contrast, an exact response theory already exists in classical statistical mechanics, that sprang from to the fluctuation theorems [6,26,27,28,29], a generalization of the second law of thermodynamics for small systems, allowing the dissipation to become a random variable that can take negative values. A quite general classical exact response theory has been built on the basis of the transient time correlation function (TTCF) [7], and of the fluctuation relations. The key ingredient of this theory is the Dissipation Function, first defined by Evans and Searles as the physically relevant quantity that verifies the fluctuation relation [30,31,32]. Its connection with exact response theory is referred to as the dissipation theorem in Ref. [33]; it concerns systems subjected to time-independent perturbations, such as an external constant field. A dynamical system interpretation has been provided [34]. This theory has offered numerous advantages and new perspectives, and it has been applied in several contexts [35,36]. It has allowed treating hard non-equilibrium problems, at low drivings [37,38], drastically improving the signal-to-noise ratio, and providing a superior method with respect to direct averaging for such calculations. In recent applications, it was demonstrated that this method dramatically improves the accuracy of the results at low shear rates, and that it is suitable to investigate atomistically detailed confined fluids at realistic flow rates [37]. It has also been shown that the TTCF can be used to define a local diffusion coefficient leading to important practical implications for nanoscale and inhomogeneous systems [39]. More recently, the Dissipation Function has been applied to polar molecules in an electric field, yielding excellent results [40]. Machine learning techniques have further advanced the study of the Dissipation Function in nonequilibrium steady states, leading to a more accurate, short-time valid steady-state fluctuation theorem [41]. Moreover, this exact response theory has recently been extended to more complex perturbations, such as stochastic ones [36]. The Dissipation Function and related response theory remain subjects of active research, and further interesting developments are expected in various fields.

Attempts to link quantum response theory with quantum fluctuation theorems have been explored but without advancing beyond perturbative approaches [42]. Therefore, in the present paper, we propose an exact quantum response theory based on the Dissipation Function. In Section 2, we review the fundamental aspects of the classical exact response theory. We formulate a quantum analogue of the Dissipation Function, introducing possible definitions for the corresponding Dissipation Operators. Finally, we derive two exact expressions for the computation of the observables expectation values, based on the new Dissipation Operator, analogously to the corresponding classical response theory. In Section 5 we apply the new expressions to qubit systems, where exact results can be obtained in other ways as well. We then compare the new approach with linear response theory in an application to a spin-1/2 particle in a magnetic field. This example highlights the advantages that an exact response theory can offer. In Section 6, we conduct a study of the Dissipation Operator, providing some results that could be useful in future research. In Section 7, we extend the new method to Lindblad dynamics, and apply it to an open qubit system characterized by decoherence. Under appropriate assumptions, the new expression based on the Dissipation Operator gives exact results for this particular case. In Section 8, we discuss the obtained results and suggest future developments.

2. Classical and Quantum Response Theory

In this section, we review the fundamental aspects of the classical Dissipation Function and its use in response theory [34]. We then propose a formulation for quantum dynamics.

Let us consider a system whose microscopic phase $\Gamma \in \mathcal{M}$ evolves according to the equation of motion $\dot{\Gamma }=G(\Gamma )$. We define the map $S^t:\mathcal{M}\to \mathcal{M}$. $S^t\Gamma$ is the solution at time t of the system with initial condition $\Gamma \in \mathcal{M}$. We assume that the phase space $\mathcal{M}$ is endowed with a probability measure $d{\mu }_0(\Gamma )=f_0(\Gamma )d\Gamma$, of density $f_0$, which evolves according to the generalized Liouville equation ${\partial }_t{f}_t(\Gamma )=-{\nabla }_{\Gamma }·\left(f_{t}G(\Gamma )\right)$. This can be rewritten in terms of the Dissipation Function [34]

$${\Omega }^{f_t}(\Gamma )\equiv -[\Lambda (\Gamma )+G(\Gamma )·{\nabla }_{\Gamma }ln{f}_t(\Gamma )]$$

(1)

as

$${\displaystyle \frac{\partial f_t}{\partial t}}(\Gamma )=f_t{\Omega }^{f_t}(\Gamma ),$$

(2)

where $\Lambda ={\nabla }_{\Gamma }·G$ is the phase space variation rate. The evolution of the ensemble average, defined by:

$${\langle \mathcal{O}\rangle }_t={\int }_{\mathcal{M}}\mathcal{O}(\Gamma )f_t(\Gamma )d\Gamma$$

(3)

can be expressed with respect to the initial distribution $f_0$ as:

$${\langle \mathcal{O}\rangle }_t={\langle \mathcal{O}\rangle }_0+{\int }_0^t{〈(\mathcal{O}\circ S^s){\Omega }^{f_0}〉}_{0}ds.$$

(4)

See, e.g., Refs.[7,33,34,43] for detailed derivations. This result allows us to calculate the system’s response to an external perturbation in an exact, not approximate, way. Here, $f_0$ is the unperturbed distribution, which is usually the equilibrium ensemble for the unperturbed dynamics, while $S^t$ represents the exact (perturbed) dynamics. The origin of the name Dissipation Function comes from the fact that this quantity, for nonequilibrium molecular dynamics (NEMD), corresponds precisely to the energy dissipation rate: product of dissipative force times associated flux. This is clear by comparing expression (4) with those obtained using the TTCF for a thermodynamic system under the influence of an external field $F_{ext}$ [33]:

$${\langle \mathcal{O}\rangle }_t={\langle \mathcal{O}\rangle }_0-{\displaystyle \frac{V}{k_{B}T}}{\int }_0^t{〈(\mathcal{O}\circ S^s)J·F_{ext}〉}_{0}ds.$$

(5)

Another interesting interpretation of $\Omega$ arises in the more general context of dynamical systems. Let ${\Omega }_{0,s}^{f_0}$ denote the time integral of the Dissipation Function between the time instants 0 and s:

$${\Omega }_{0,s}^{f_0}(\Gamma )\equiv {\int }_0^s{\Omega }^{f_0}(S^u\Gamma )du$$

(6)

Then, one gets:

$${〈{\Omega }_{0,s}^{f_0}{\rangle }_0=\langle ln(f_0/f_{-s})〉}_0\ge 0$$

(7)

that is the relative entropy $D(f_0\parallel f_{-s})$[34], or Kullback–Leibler divergence of the distributions $f_0$ and $f_{-s}$. This is of interest, for instance, in large deviation theory, where ${\Omega }^f$ plays the role of the large deviation functional, and has numerous consequences in applications [44]. Expression (4) also offers other advantages. Firstly, it keeps the probability fixed while allowing only the observables to evolve over time. Probability evolution requires the reversed dynamics, which is more difficult to use, whereas observables do not. For equation (4), the dynamics are assumed to be invertible [34], although not necessarily time-reversal invariant, as often required in statistical mechanics. The use of the equilibrium distribution $f_0$ is fortunate because it is known and has the property of smoothing the result, effectively improving the signal to noise ratio [37]. In linear response theory one uses $f_0$ as well, but in that case the result is approximate rather than exact. Additionally, the notion of T-mixing, which is a necessary condition for the Fluctuation Relations to hold in non equilibrium steady states [30], through its connection with the Dissipation Function ($\Omega T-$mixing),

$${〈{\Omega }^{f_0}(\mathcal{O}\circ S^t)〉}_0\to {〈{\Omega }^{f_0}〉}_0{〈\mathcal{O}\circ S^t〉}_0=0\mathrm{for}t\to \infty ,$$

(8)

provides a new approach for describing relaxation towards equilibrium [34].

The objective of this paper is to develop the quantum counterpart of the classical expression (4), since it could allow the extension to quantum mechanics of the benefits of the classical response theory. For sake of simplicity, we focus on finite-dimensional Hilbert spaces, allowing us to concentrate on constructing the new formalism while avoiding the technicalities of infinite dimensional analysis.

In quantum mechanics, a system is completely determined by its quantum state $|\psi \rangle$ which obeys the Schrödinger equation [45]:

$$i\hslash {\displaystyle \frac{\partial }{{\partial }_t}}\psi (x,t)=\mathcal{H}\psi (x,t)$$

(9)

a deterministic evolution of a probabilistic entity. In this case, one assumes that all degrees of freedom of interest are represented in the Hamiltonian $\mathcal{H}$, hence one refers to an isolated system. Observables are then expressed by self-adjoint operators [46]. The analogous notion of classical ensembles is here given by collections of states $|{\psi }_j\rangle$, each taken with probability $p_j$, with j running over a suitable set. Collectively, these states are represented by the density operator

$$\rho \equiv \sum _j{p}_j|{\psi }_j\rangle \langle {\psi }_j|$$

(10)

suitable to treat both pure and mixed states [46]. This operator is self-adjoint, ${\rho }^{†}=\rho$, has a unitary trace, $Tr\left(\rho \right)=1$, and it is semi-definite positive, $\langle \psi \left|\rho \right|\psi \rangle \ge 0\forall \psi$. The expectation value of an observable $\mathcal{O}$ is expressed through the trace operation, $\langle \mathcal{O}\rangle =Tr\left(\rho \mathcal{O}\right)$, that replaces the integral (3) of the classical theory.

There are two equivalent pictures [45]. In Schrödinger picture, observables are time independent and $\rho$ evolves over time according to Von Neumann (Quantum Liouville) equation:

$${\displaystyle \frac{\partial }{\partial t}}{\rho }_t={\displaystyle \frac{1}{i\hslash }}[\mathcal{H},{\rho }_t].$$

(11)

Instead, in the Heisenberg picture, $\rho$ remains fixed and the observables evolve in time according to

$${\displaystyle \frac{d}{dt}}{A}_t={\displaystyle \frac{1}{i\hslash }}[A_t,\mathcal{H}]+{\partial }_{t}A.$$

(12)

If the Hamiltonian is constant, the time-propagator $U\left(t\right)$ is defined as

$$U\left(t\right)\equiv e^{-{\displaystyle \frac{it}{\hslash }}\mathcal{H}},U^{†}\left(t\right)=e^{+{\displaystyle \frac{it}{\hslash }}\mathcal{H}}.$$

(13)

This time propagator is a one-parameter unitary group [46], i.e., it is unitary, $U\left(t\right)U^{†}\left(t\right)=U^{†}\left(t\right)U\left(t\right)=Id$, it satisfies the group property $U(t+s)=U\left(t\right)U\left(s\right)$ and $U\left(0\right)=Id$. Then, if the Hamiltonian is time independent, the solutions of the previous equations can be written as

$$\left(\mathrm{Schr}\text{ö}\mathrm{dinger}\right){\rho }_t=U\left(t\right){\rho }_0{U}^{†}\left(t\right);(\mathrm{Heisenberg})A_t=U^{†}\left(t\right)AU\left(t\right).$$

(14)

These are the main tools we will use in the following.

3. Quantum Dissipation Function for Time Indpendent Perturbations

We formulate the quantum exact response theory based on the Dissipation Function starting from scratch, so it is good to start by considering simple cases. Take a time independent Hamiltonian ${\mathcal{H}}_0$ and an arbitrary initial density operator ${\rho }_0$. Introduce a perturbation, producing a new Hamiltonian:

$$\mathcal{H}={\mathcal{H}}_0+\theta \left(t\right)\lambda {\mathcal{H}}_{ext},$$

(15)

where $\theta \left(t\right)$ is the Heaviside function. ${\mathcal{H}}_0$ represents the equilibrium dynamics, and usually ${\rho }_0$ is the corresponding equilibrium density operator, that is stationary with respect to ${\mathcal{H}}_0$:

$$[{\mathcal{H}}_0,{\rho }_0]=0.$$

(16)

The external perturbation $\lambda {\mathcal{H}}_{ext}$ is assumed to be turned on at time $t=0$, and then kept constant over time, with $\lambda \in \mathbb{R}$. ${\mathcal{H}}_{ext}$ makes ${\rho }_0$ no longer be invariant, and the expectation values of the observables become time-dependent. Trying to translate the classical notion of Dissipation Function in this framework, we find there are some issues to tackle, such as the ordering problem. This requires us to choose a proper form for the Dissipation Function. From the classical expression (2), we obtain:

$${\Omega }^{f_0}(\Gamma )=f_0^{-1}(\Gamma ){\displaystyle \frac{\partial f_0(\Gamma )}{\partial t}}.$$

(17)

that holds only in the part of phase space in which $f_0(\Gamma )>0$. This is reminiscent of the more general notion of absolute continuity of the evolved distributions with respect to the initial [47], and has been called ergodic consistency in molecular dynamics [43,48]. This limitation is not serious in many applications because the initial distribution usually corresponds to an unperturbed equilibrium state and does not vanish anywhere in the phase space of the perturbed evolution, if this satisfies the same constraints of the equilibrium dynamics [34,48].

The quantum mechanical counterpart of the classical Dissipation Function can then be guessed, using the von Neumann equation (11), to be defined as follows.

Definition 1.

Let ${\rho }_0$ be an initial density operator and $\mathcal{H}={\mathcal{H}}_0+\lambda {\mathcal{H}}_{ext}$ a time-independent Hamiltonian operator, the Dissipation Operator can be defined as

$${\Omega }^0\equiv {\displaystyle \frac{1}{i\hslash }}{\rho }_0^{-1}[\mathcal{H},{\rho }_0].$$

(18)

As there are different equivalent definitions of the classical Dissipation Function, our choice is to some extent arbitrary, and other choices can be considered. We are going to show, however, that this form of ${\Omega }^0$ is consistent with linear response for small perturbations, and yields the exact solutions of problems whose exact solution is known. It remains that this definition faces two issues. First, the Dissipation Operator (18) is not self-adjoint, so its expectation value can be a complex number, ${\langle {\Omega }^0\rangle }_t\in \mathbb{C}$. For isolated systems, this may be a challenge, since it does not represent a directly measurable quantity. However, this is not a problem per se in the case of open systems, where relying solely on self-adjoint operators may not be possible or necessary [49,50]. While it would be preferable to have a self-adjoint ${\Omega }^0$, this is not a priority. The second issue is that ${\rho }_0$ is not always invertible, as in the case for pure states. There are, however, ways to handle such a difficulty. One approach is to reduce the dimension of the Hilbert space, if one knows that the evolution remains within a subspace in which ${\rho }_0$ is invertible. Another way is to add arbitrarily small amounts $p_{ϵ}$ to part of the entries of ${\rho }_0$ so that its rank becomes full, and later analyze the results in the $p_{ϵ}\to 0$ limit. In any event, it is the same issue known in classical mechanics, and, analogously, the problem is solved if the initial distribution has support wider than that of the evolved distributions [47,48]. Despite these issues, definition (18) has various advantages. It allows us to derive a quantum expression analogous to equation (4). Moreover, it satisfies two important properties already present in the classical context of the Dissipation Function.

Proposition 1.

The expectation value of the Dissipation Operator calculated with respect to the density operator ${\rho }_0$ is always zero, and its initial time derivative is always positive,

$${\langle {\Omega }^0\rangle }_0=0,{\displaystyle \frac{d}{ds}}{\langle {\Omega }^0\rangle }_s{|}_{s=0}\ge 0.$$

(19)

Proof of Proposition 1.

$${\langle {\Omega }^0\rangle }_0=Tr\left({\rho }_0\left({\displaystyle \frac{i}{\hslash }}{\rho }_0^{-1}[\mathcal{H},{\rho }_0]\right)\right)={\displaystyle \frac{i}{\hslash }}Tr(Id·[\mathcal{H},{\rho }_0])$$

(20)

now we can use first the linear property and then the cyclic property of the trace,

$$\begin{array}{c}\hfill {\displaystyle \frac{i}{\hslash }}Tr(\mathcal{H}{\rho }_0-{\rho }_0\mathcal{H})={\displaystyle \frac{i}{\hslash }}\left(Tr\left(\mathcal{H}{\rho }_0\right)-Tr\left({\rho }_0\mathcal{H}\right)\right)={\displaystyle \frac{i}{\hslash }}\left(Tr\left(\mathcal{H}{\rho }_0\right)-Tr\left(\mathcal{H}{\rho }_0\right)\right)=0\end{array}$$

(21)

We conclude that ${\langle {\Omega }^0\rangle }_0=0$.

We can show the second property by making the incremental limit explicit in the derivative operation,

$$\begin{array}{c}\hfill {\displaystyle \frac{d}{ds}}{\langle {\Omega }^0\rangle }_s{|}_{s=0}=\underset{s\to 0}{lim}{\displaystyle \frac{1}{s}}[{\langle {\Omega }^0\rangle }_s-{\langle {\Omega }^0\rangle }_0]=\underset{s\to 0}{lim}{\displaystyle \frac{1}{s}}\left[Tr\left({\rho }_s{\Omega }^0\right)-Tr\left({\rho }_0{\Omega }^0\right)\right]\\ \hfill =Tr\left(\underset{s\to 0}{lim}{\displaystyle \frac{{\rho }_s-{\rho }_0}{s}}{\Omega }^0\right)=Tr\left({\displaystyle \frac{\partial {\rho }_0}{\partial t}}{\Omega }^0\right)\end{array}$$

(22)

We now use a result that we will show later, eq. (43): we can also use the adjoint of the Dissipation Operator, ${\left({\Omega }^0\right)}^{†}$, to express ${\partial }_t{\rho }_0$,

$${\displaystyle \frac{\partial {\rho }_0}{\partial t}}={\left({\Omega }^0\right)}^{†}{\rho }_0$$

(23)

Substituting this into the last expression, we obtain

$${\displaystyle \frac{d}{ds}}{\langle {\Omega }^0\rangle }_s{|}_{s=0}=Tr\left({\displaystyle \frac{\partial {\rho }_0}{\partial t}}{\Omega }^0\right)=Tr\left({\left({\Omega }^0\right)}^{†}{\rho }_0{\Omega }^0\right)\ge 0$$

(24)

Where the last inequality holds for the semi-positivity of ${\rho }_0$. In fact, let $A\in {\mathbb{C}}^{n\times n}$ be a self-adjoint and semi-definite positive matrix and $C\in {\mathbb{C}}^{n\times m}$ be an arbitrary matrix, then:

$$\mathrm{Tr}\left(C^{†}AC\right)=\sum _{i=1}^m{\left(C^{†}AC\right)}_{ii}=\sum _{i=1}^m{c}_i^{†}A{c}_i\ge 0$$

(25)

where $c_i$ are the columns of C, and the last inequality follows from the definition of semi-definite positive matrix. □

We now derive one of the key result of this paper: an exact expression for the expectation value of observables, based on the new Dissipation Operator (18).

Proposition 2.

Let ${\rho }_0$ be the initial density operator and $\mathcal{H}={\mathcal{H}}_0+\lambda {\mathcal{H}}_{ext}$ the time-independent Hamiltonian operator. The expectation value of any observable $\mathcal{O}$ can be calculated exactly using the Dissipation Operator as

$${\langle \mathcal{O}\rangle }_t={\langle \mathcal{O}\rangle }_0+{\int }_0^t{\langle {\Omega }^0{\mathcal{O}}_s\rangle }_{0}ds$$

(26)

Proof of Proposition 2.

Let’s derive a useful expression for ${\displaystyle \frac{d}{ds}}{\langle \mathcal{O}\rangle }_s$:

$$\begin{array}{c}\hfill {\displaystyle \frac{d}{ds}}{\langle \mathcal{O}\rangle }_s=\underset{h\to 0}{lim}{\displaystyle \frac{1}{h}}\left[{\langle \mathcal{O}\rangle }_{s+h}-{\langle \mathcal{O}\rangle }_s\right]=\underset{h\to 0}{lim}{\displaystyle \frac{1}{h}}\left[Tr\left(\mathcal{O}{\rho }_{s+h}\right)-Tr\left(\mathcal{O}{\rho }_s\right)\right]\\ \hfill =\underset{h\to 0}{lim}{\displaystyle \frac{1}{h}}\left[Tr\left(\mathcal{O}U(s+h){\rho }_0{U}^{†}(s+h)\right)-Tr\left(\mathcal{O}U\left(s\right){\rho }_0{U}^{†}\left(s\right)\right)\right]\\ \hfill =\underset{h\to 0}{lim}{\displaystyle \frac{1}{h}}\left[Tr\left(U^{†}(s+h)\mathcal{O}U(s+h){\rho }_0\right)-Tr\left(U^{†}\left(s\right)\mathcal{O}U\left(s\right){\rho }_0\right)\right].\end{array}$$

(27)

Thanks to the group property of the operator $U\left(t\right)$ for time independent Hamiltonians,

$$U(s+h)=U\left(h\right)U\left(s\right),U^{†}(s+h)=U^{†}\left(s\right)U^{†}\left(h\right),$$

(28)

thanks to the fact that observables and density matrix evolve "at opposite times", and introducing the evolved operators ${\mathcal{O}}_s=U^{†}\left(s\right)\mathcal{O}U\left(s\right)$ we can write:

$$\begin{array}{ccc}& & \underset{h\to 0}{lim}{\displaystyle \frac{1}{h}}\left[Tr\left(U^{†}\left(h\right){\mathcal{O}}_{s}U\left(h\right){\rho }_0\right)-Tr\left({\mathcal{O}}_s{\rho }_0\right)\right]\hfill \end{array}$$

(29)

$$\begin{array}{ccc}& =& \underset{h\to 0}{lim}{\displaystyle \frac{1}{h}}\left[Tr\left({\mathcal{O}}_{s}U\left(h\right){\rho }_0{U}^{†}\left(h\right)\right)-Tr\left({\mathcal{O}}_s{\rho }_0\right)\right]\hfill \end{array}$$

(30)

$$\begin{array}{ccc}& =& \underset{h\to 0}{lim}{\displaystyle \frac{1}{h}}[Tr\left({\mathcal{O}}_s\left(U\left(h\right){\rho }_0{U}^{†}\left(h\right))-{\rho }_0\right)\right]\hfill \end{array}$$

(31)

$$\begin{array}{ccc}& =& Tr\left({\mathcal{O}}_s\underset{h\to 0}{lim}{\displaystyle \frac{{\rho }_h-{\rho }_0}{h}}\right)=Tr\left({\mathcal{O}}_s{\displaystyle \frac{\partial {\rho }_0}{\partial t}}\right)\hfill \end{array}$$

(32)

From the definition (18) of the Dissipation Operator we have:

$${\displaystyle \frac{\partial {\rho }_0}{\partial t}}={\rho }_0{\Omega }^0$$

(33)

and substituting in (32) we obtain:

$$Tr\left({\mathcal{O}}_s{\displaystyle \frac{\partial {\rho }_0}{\partial t}}\right)=Tr\left({\mathcal{O}}_s{\rho }_0{\Omega }^0\right)=Tr\left({\rho }_0{\Omega }^0{\mathcal{O}}_s\right)={\langle {\Omega }^0{\mathcal{O}}_s\rangle }_0.$$

(34)

Then we have:

$${\displaystyle \frac{d}{ds}}{\langle \mathcal{O}\rangle }_s={\langle {\Omega }_0{\mathcal{O}}_s\rangle }_0.$$

(35)

and

$${\langle \mathcal{O}\rangle }_t={\langle \mathcal{O}\rangle }_0+{\int }_0^t{\displaystyle \frac{d}{ds}}{\langle \mathcal{O}\rangle }_{s}ds={\langle \mathcal{O}\rangle }_0+{\int }_0^t{〈{\Omega }^0{\mathcal{O}}_s〉}_{0}ds$$

(36)

□

First, we note the strong similarity with the analogous expression (4) in classical statistical mechanics. The system’s response is expressed in terms of the correlation function of the Dissipation Operator ${\Omega }^0$ and the observable ${\mathcal{O}}_t$, evolved according to the exact dynamics $\mathcal{H}$, calculated with respect to the equilibrium density operator ${\rho }_0$. This use of the initial distribution is common to linear response, but like in the classical case, eq.(26) is not an approximate expression and is not limited to small perturbations.

In general, ${\Omega }^0$ and $\mathcal{O}$ don’t commute, ${\langle {\Omega }^0{\mathcal{O}}_s\rangle }_0\ne {\langle {\mathcal{O}}_s{\Omega }^0\rangle }_0.$ Therefore, we cannot use ${\langle {\mathcal{O}}_s{\Omega }^0\rangle }_0$ in equation (26). In quantum mechanics, other types of correlations are often used, such as Kubo’s canonical correlation or symmetric correlation [3], in order to obtain real numbers as results, but for expression (26), this is automatically obtained. Thus, the meaning of quantum response function can be assigned to the form used in equation (26). Nevertheless, we now consider a symmetrized version of our response formula introducing the following notation:

$${\langle A;B\rangle }_0=Tr\left[{\rho }_0{\displaystyle \frac{AB+BA}{2}}\right].$$

(37)

4. Self Adjoint Quantum Dissipation Operator

For the isolated time independent dynamics considered above, the response of an observable ${\langle \mathcal{O}\rangle }_t$ is computable in the Heisenberg and Schrödinger pictures without any need to introduce the Dissipation Function. In this sense, expression (26) only represents a different formalism for expressing a well-known result. However, the different perspective it offers, it may turn useful in the solution of complex time independent problems, as it happens in classical mechanics, especially in presence of poor signal to noise dynamics [37]. Furthermore, its analogy with the classical expression (4) may prove useful in time dependent situations. Additionally, as in equation (5), there may be opportunities to associate the Dissipation Operator ${\Omega }^0$ with the production of generalized entropy, involving, for instance the notion of Kullback-Leibler divergence (7). These interpretations are supported by the fact that Proposition 1 seems to imply the existence of $h>0$ such that:

$${\langle {\Omega }_{0,s}^0\rangle }_0\ge 0for0<s<h.$$

(38)

Finally, it could be possible to use the $\Omega T-$mixing property [34,43] to study the relaxation of observables toward stationary states. This question would be especially interesting in the context of open quantum systems. Choosing $\mathcal{O}=1$, we find an excellent consistency condition: the probabilistic interpretation of the density operator is preserved using expression (26).

To explore this possibility, let us start observing that the minimum requirements of response theory are satisfied by the notions introduced above. In particular, Proposition 1 implies:

$${\langle 1\rangle }_t={\langle 1\rangle }_0+{\int }_0^t{\langle {\Omega }^{0}1\rangle }_{0}ds=1+0t=1.$$

(39)

and, for a perturbation ${\mathcal{H}}_{ext}$ that commutes with ${\rho }_0$, we have ${\Omega }^0=0$, hence:

$${\langle \mathcal{O}\rangle }_t={\langle \mathcal{O}\rangle }_0+{\int }_0^t{\langle 0\rangle }_{0}ds={\langle \mathcal{O}\rangle }_0$$

(40)

for all observables, as desired.

In case we take $\mathcal{O}={\Omega }^0$, and we wish to obtain a real number, as the dissipation should be, we have however a difficulty: ${\Omega }^0$ is not self adjoint and ${\langle {\Omega }^0\rangle }_t$ may be complex. While this is not necessarily a problem for the theory of open systems, it is interesting to develop a self-adjoint Dissiaption Operator. Therefore, we propose the following symmetrized operator:

$${\tilde{\Omega }}^0\equiv {\displaystyle \frac{{\Omega }^0+{\left({\Omega }^0\right)}^{†}}{2}}.$$

(41)

as Dissipation Function. Then, we note that:

$$\begin{array}{ccc}\hfill {\left({\Omega }^0\right)}^{†}& =& {\left({\displaystyle \frac{1}{i\hslash }}{\rho }_0^{-1}[\mathcal{H},{\rho }_0]\right)}^{†}={\left({\displaystyle \frac{1}{i\hslash }}{\rho }_0^{-1}\mathcal{H}{\rho }_0-{\displaystyle \frac{1}{i\hslash }}\mathcal{H}\right)}^{†}=-{\displaystyle \frac{1}{i\hslash }}{\left({\rho }_0^{-1}\mathcal{H}{\rho }_0\right)}^{†}+{\displaystyle \frac{1}{i\hslash }}{\mathcal{H}}^{†}\hfill \\ & =& -{\displaystyle \frac{1}{i\hslash }}\left({\rho }_0\mathcal{H}{\rho }_0^{-1}-\mathcal{H}\right)=-{\displaystyle \frac{1}{i\hslash }}\left[{\rho }_0,\mathcal{H}\right]{\rho }_0^{-1}={\displaystyle \frac{1}{i\hslash }}\left[\mathcal{H},{\rho }_0\right]{\rho }_0^{-1}\hfill \end{array}$$

(42)

so that:

$${\displaystyle \frac{\partial {\rho }_0}{\partial t}}={\left({\Omega }^0\right)}^{†}{\rho }_0$$

(43)

which leads to the following.

Proposition 3.

By means of the self-adjoint operator (41), eq. (26) for time independent perturbations can be expressed as:

$${\langle \mathcal{O}\rangle }_t={\langle \mathcal{O}\rangle }_0+{\int }_0^t\left({\langle {\tilde{\Omega }}^0;{\mathcal{O}}_s\rangle }_0+{\displaystyle \frac{1}{2}}{\langle [\overline{{\Omega }^0},{\mathcal{O}}_s]\rangle }_0\right)ds,$$

(44)

where $\overline{{\Omega }^0}$ is the anti-self-adjoint operator $\overline{{\Omega }^0}\equiv {\displaystyle \frac{1}{2}}\left({\Omega }^0-{\left({\Omega }^0\right)}^{†}\right).$

Proof. 

From the proof of (26) we know

$${\displaystyle \frac{d}{ds}}{\langle \mathcal{O}\rangle }_s=Tr\left({\mathcal{O}}_s{\displaystyle \frac{\partial {\rho }_0}{\partial t}}\right)$$

(45)

To introduce the Hermitian Dissipation Operator ${\tilde{\Omega }}^0$, we use:

$${\displaystyle \frac{\partial {\rho }_0}{\partial t}}={\rho }_0{\Omega }^0;{\displaystyle \frac{\partial {\rho }_0}{\partial t}}={\left({\Omega }^0\right)}^{†}{\rho }_0.$$

(46)

Therefore, we can write:

$$\begin{array}{ccc}\hfill Tr\left({\mathcal{O}}_s{\displaystyle \frac{\partial {\rho }_0}{\partial t}}\right)& =& {\displaystyle \frac{1}{2}}Tr\left({\mathcal{O}}_s{\displaystyle \frac{\partial {\rho }_0}{\partial t}}\right)+{\displaystyle \frac{1}{2}}Tr\left({\mathcal{O}}_s{\displaystyle \frac{\partial {\rho }_0}{\partial t}}\right)\hfill \\ & =& {\displaystyle \frac{1}{2}}Tr\left({\mathcal{O}}_s{\rho }_0{\Omega }^0\right)+{\displaystyle \frac{1}{2}}Tr\left({\mathcal{O}}_s{\left({\Omega }^0\right)}^{†}{\rho }_0\right)={\displaystyle \frac{1}{2}}Tr\left({\Omega }^0{\mathcal{O}}_s{\rho }_0+{\left({\Omega }^0\right)}^{†}{\rho }_0{\mathcal{O}}_s\right).\hfill \end{array}$$

(47)

We can rewrite this last expression as

$$\begin{array}{ccc}\hfill & & Tr\left({\Omega }^0{\mathcal{O}}_s{\rho }_0+{\left({\Omega }^0\right)}^{†}{\rho }_0{\mathcal{O}}_s\right)=Tr\left({\Omega }^0{\mathcal{O}}_s{\rho }_0+{\left({\Omega }^0\right)}^{†}({\mathcal{O}}_s{\rho }_0+[{\rho }_0,{\mathcal{O}}_s])\right)\hfill \\ & & =Tr\left(({\Omega }^0+{\left({\Omega }^0\right)}^{†}){\mathcal{O}}_s{\rho }_0\right)+Tr\left({\left({\Omega }^0\right)}^{†}[{\rho }_0,{\mathcal{O}}_s]\right)=2Tr\left({\tilde{\Omega }}^0{\mathcal{O}}_s{\rho }_0\right)+Tr\left({\left({\Omega }^0\right)}^{†}[{\rho }_0,{\mathcal{O}}_s]\right).\hfill \end{array}$$

(48)

Where the last equality follows from the definition (41). We would like to eliminate the presence of ${\left({\Omega }^0\right)}^{†}$. For this purpose we note that with similar procedures, but using the commutator property for the first term in the trace, we can equivalently write:

$$Tr\left({\Omega }^0{\mathcal{O}}_s{\rho }_0+{\left({\Omega }^0\right)}^{†}{\rho }_0{\mathcal{O}}_s\right)=Tr\left({\Omega }^0{\rho }_0{\mathcal{O}}_s+{\Omega }^0[{\mathcal{O}}_s,{\rho }_0]+{\left({\Omega }^0\right)}^{†}{\rho }_0{\mathcal{O}}_s\right).$$

(49)

We find

$$Tr\left({\Omega }^0{\mathcal{O}}_s{\rho }_0+{\left({\Omega }^0\right)}^{†}{\rho }_0{\mathcal{O}}_s\right)=2Tr\left({\tilde{\Omega }}^0{\rho }_0{\mathcal{O}}_s\right)+Tr\left({\Omega }^0[{\mathcal{O}}_s,{\rho }_0]\right).$$

(50)

Now, we combine equations (48) and (50). Starting from equation (47), we have

$$Tr\left({\mathcal{O}}_s{\displaystyle \frac{\partial {\rho }_0}{\partial t}}\right)={\displaystyle \frac{1}{4}}\left(Tr\left({\Omega }^0{\mathcal{O}}_s{\rho }_0+{\left({\Omega }^0\right)}^{†}{\rho }_0{\mathcal{O}}_s\right)+Tr\left({\Omega }^0{\mathcal{O}}_s{\rho }_0+{\left({\Omega }^0\right)}^{†}{\rho }_0{\mathcal{O}}_s\right)\right).$$

(51)

We can use equation (48) in place of the first term and equation (50) in place of the second one, obtaining

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{4}}\left(2Tr\left({\tilde{\Omega }}^0{\mathcal{O}}_s{\rho }_0\right)+Tr\left({\left({\Omega }^0\right)}^{†}[{\rho }_0,{\mathcal{O}}_s]\right)+2Tr\left({\tilde{\Omega }}^0{\rho }_0{\mathcal{O}}_s\right)+Tr\left({\Omega }^0[{\mathcal{O}}_s,{\rho }_0]\right)\right)\\ \hfill ={\displaystyle \frac{1}{2}}Tr\left({\tilde{\Omega }}^0{\mathcal{O}}_s{\rho }_0\right)+{\displaystyle \frac{1}{2}}Tr\left({\tilde{\Omega }}^0{\rho }_0{\mathcal{O}}_s\right)+{\displaystyle \frac{1}{4}}\left(Tr\left({\left({\Omega }^0\right)}^{†}[{\rho }_0,{\mathcal{O}}_s]+{\Omega }^0[{\mathcal{O}}_s,{\rho }_0]\right)\right)\\ \hfill ={\displaystyle \frac{1}{2}}Tr\left({\tilde{\Omega }}^0{\mathcal{O}}_s{\rho }_0\right)+{\displaystyle \frac{1}{2}}Tr\left({\mathcal{O}}_s{\tilde{\Omega }}^0{\rho }_0\right)+{\displaystyle \frac{1}{4}}\left(Tr\left(({\Omega }^0-{\left({\Omega }^0\right)}^{†})[{\mathcal{O}}_s,{\rho }_0]\right)\right).\end{array}$$

(52)

We can clearly use the symmetric correlation (37) for terms with ${\tilde{\Omega }}^0$ but we are unable to completely replace ${\Omega }^0$ and ${\left({\Omega }^0\right)}^{†}$ with the hermitian ${\tilde{\Omega }}^0$. We can also observe that ${\Omega }^0-{\left({\Omega }^0\right)}^{†}$ is the anti-self-adjoint part of ${\Omega }^0$:

$$\overline{{\Omega }^0}\equiv {\displaystyle \frac{1}{2}}\left({\Omega }^0-{\left({\Omega }^0\right)}^{†}\right).$$

(53)

In conclusion, we obtain

$${\displaystyle \frac{d}{ds}}{\langle \mathcal{O}\rangle }_s={\langle {\tilde{\Omega }}^0;{\mathcal{O}}_s\rangle }_0+{\displaystyle \frac{1}{2}}Tr\left({\overline{\Omega }}^0[{\mathcal{O}}_s,{\rho }_0]\right)={\langle {\tilde{\Omega }}^0;{\mathcal{O}}_s\rangle }_0+{\displaystyle \frac{1}{2}}{\langle [\overline{{\Omega }^0},{\mathcal{O}}_s]\rangle }_0.$$

(54)

Integrating as already done previously we obtain

$${\langle \mathcal{O}\rangle }_t={\langle \mathcal{O}\rangle }_0+{\int }_0^t\left({\langle {\tilde{\Omega }}^0;{\mathcal{O}}_s\rangle }_0+{\displaystyle \frac{1}{2}}{\langle [\overline{{\Omega }^0},{\mathcal{O}}_s]\rangle }_0\right)ds.$$

(55)

□

This is an alternative form of expression (26), which yields equivalent exact results. Since we are still in the early stages of formulating this theory, it is not clear which form best represents the method based on the Dissipation Operator. Therefore, we present both formulas. Equation (26) is in a simpler form, and its analogy with the classical expression is immediately apparent. On the other hand, the analogy between expression (44) and the classical counterpart is less obvious, but it has the advantage of involving only self-adjoint (and anti-self-adjoint) operators. These operators are commonly used in quantum mechanics and possess well-known properties. Now, $\langle {\tilde{\Omega }}^0\rangle \in \mathbb{R}$, and $\langle {\overline{\Omega }}^0\rangle$ is a pure imaginary number. Additionally, expression (44) is based on the symmetric correlation (37), which shares characteristics much closer to the classical one. For instance, ${\langle {\tilde{\Omega }}^0;{\mathcal{O}}_s\rangle }_0={\langle {\mathcal{O}}_s;{\tilde{\Omega }}^0\rangle }_0$. Finally, with steps similar to those of the proof of Proposition 1, we obtain:

$${\langle {\tilde{\Omega }}^0\rangle }_0={\langle \overline{{\Omega }^0}\rangle }_0=0;{\displaystyle \frac{d}{ds}}{\langle {\tilde{\Omega }}^0\rangle }_s{{|}_{s=0}\ge 0,{\displaystyle \frac{d}{ds}}{\langle \overline{{\Omega }^0}\rangle }_s|}_{s=0}=0,$$

(56)

which makes this new Dissipation Operator $\tilde{\Omega }$ also suitable as a candidate for being associated with the production of generalized entropy.

5. Applications to Qubits and Numerical Tests

In this section, we apply our expressions for the response to perturbations of qubit systems, and we compare with the results obtained using the Heisenberg picture. To perform these comparisons, we calculated the results numerically using a MATLAB code. For mathematical convenience, we set $\hslash =1$ and treat all physical quantities as dimensionless.

A quantum system with only two energy states is referred to as a qubit. Qubits play a fundamental role in quantum computing, as they are the basic units of quantum information. They are described within two-dimensional Hilbert spaces. In this context, all possible physical observables are represented by linear combinations of the identity matrix and the Pauli matrices ${\sigma }_i$, while the density operator is expressed as the Bloch vector $\overrightarrow{a}$. One has:

$$\begin{array}{c}\hfill \mathcal{H}=\overrightarrow{c}·{\overrightarrow{\sigma }}_4,\mathcal{O}=\overrightarrow{q}·{\overrightarrow{\sigma }}_4;\overrightarrow{c},\overrightarrow{q}\in {\mathbb{R}}^4.\\ \hfill \rho ={\displaystyle \frac{1}{2}}(Id+\overrightarrow{a}·\overrightarrow{\sigma }),\overrightarrow{a}\in {\mathbb{R}}^3,\left|\overrightarrow{a}\right|\le 1.\end{array}$$

(57)

where $\overrightarrow{\sigma }=({\sigma }_x,{\sigma }_y,{\sigma }_z)$ and ${\overrightarrow{\sigma }}_4=(Id,\overrightarrow{\sigma })$, and

$${\sigma }_x=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right],{\sigma }_y=\left[\begin{array}{cc}0& -i\\ i& 0\end{array}\right],{\sigma }_z=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right],Id=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right].$$

(58)

We arbitrarily chose the vectors $\overrightarrow{a},\overrightarrow{c}$ and $\overrightarrow{q}$, and then we numerically simulated the subsequent evolution. This process was repeated for several different combinations of Hamiltonian, initial density, observable, and time interval. In all cases, the expressions (26) and (44) yielded equivalent results, coinciding with those of the Heisenberg and Schrödinger pictures, which are exact. Two examples of these tests are shown in Figure 1.

Because it happens that for some simple problems approximate theories yield exact results [51], it is interesting to compare the new approach with linear response theory. Let us consider a spin-${\displaystyle \frac{1}{2}}$ particle in a magnetic field $B_z$ directed along the $z-$axis. The dynamics is described by the Hamiltonian

$${\mathcal{H}}_0={\mu }_B{B}_z{\sigma }_z={\displaystyle \frac{\hslash {\omega }_0}{2}}{\sigma }_z.$$

(59)

The system is initially in equilibrium and described by the density operator

$${\rho }_0=p_{\uparrow }|\uparrow \rangle \langle \uparrow |+p_{\downarrow }|\downarrow \rangle \langle \downarrow |=\left[\begin{array}{cc}{p}_{\uparrow }& 0\\ 0& p_{\downarrow }\end{array}\right]$$

(60)

which is invariant under the unperturbed Hamiltonian, $[{\mathcal{H}}_0,{\rho }_0]=0$. The equilibrium dynamics is then disturbed by the interaction Hamiltonian ${\mathcal{H}}_{ext}$:

$$\mathcal{H}={\mathcal{H}}_0+\lambda {\mathcal{H}}_{ext}={\displaystyle \frac{{\omega }_0}{2}}{\sigma }_z+\lambda ({\displaystyle \frac{{\omega }_x}{2}}{\sigma }_x-{\displaystyle \frac{{\omega }_0}{2}}{\sigma }_z)=(1-\lambda ){\displaystyle \frac{{\omega }_0}{2}}{\sigma }_z+\lambda {\displaystyle \frac{{\omega }_x}{2}}{\sigma }_x$$

(61)

decreasing the z-direction component of the magnetic field, $B_z$, and producing a non-vanishing x-direction component $B_x$.

Linear response [3], considers a perturbed dynamics of the form ${\mathcal{H}}_t=\mathcal{H}-K\left(t\right)\mathcal{A}$, where $K\left(t\right)$ is a time dependent external force applied from the infinite past, $t=-\infty$, when the system was at thermal equilibrium and described by an equilibrium density matrix ${\rho }_e$, i.e. $\rho (-\infty )={\rho }_e$. A is a dynamical quantity conjugate to the applied force K. For this system, the linear response formula is

$$\begin{array}{c}\hfill {\langle B\rangle }_t={\langle B\rangle }_{{\rho }_e}+\Delta \overline{B}\left(t\right),\Delta \overline{B}\left(t\right)={\int }_{-\infty }^{t}d{t}^{\prime }K\left(t^{\prime }\right){\varphi }_{BA}(t-t^{\prime }),{\varphi }_{BA}\left(t\right)={\displaystyle \frac{1}{i\hslash }}{\langle [A,B\left(t\right)]\rangle }_{{\rho }_e};\end{array}$$

(62)

where B is an arbitrary observable and $B\left(t\right)=e^{{\displaystyle \frac{i}{\hslash }}\mathcal{H}t}B{e}^{-{\displaystyle \frac{i}{\hslash }}\mathcal{H}t}$ is the Heisenberg-evolved operator according to the unperturbed dynamics $\mathcal{H}$ (different from the total dynamics ${\mathcal{H}}_t$). The dynamics in (61) is simpler than this, as we are considering a constant disturbance over time. We can assume the disturbance is absent before time $t=0$ and is impulsively turned on at $t=0$. Let us apply a change of notation to match ours:

$$A\to -{\mathcal{H}}_{ext};K\left(t\right)\to \theta \left(t\right)\lambda ;{\mathcal{H}}_t=\mathcal{H}-K\left(t\right)\mathcal{A}\to \mathcal{H}={\mathcal{H}}_0+\lambda {\mathcal{H}}_{ext};B\to \mathcal{O}.$$

(63)

where $\theta \left(t\right)$ is the Heaviside function. Kubo’s formula can now be rewritten for our problem and simplified as

$$\begin{array}{c}\hfill {\langle \Delta \mathcal{O}\rangle }_t={\int }_{-\infty }^{t}d{t}^{\prime }\theta \left(t^{\prime }\right)\lambda \varphi (t-t^{\prime })={\int }_0^{t}d{t}^{\prime }\lambda \varphi (t-t^{\prime })\\ \hfill =-{\int }_t^{0}d{t}^{*}\lambda \varphi \left(t^{*}\right)={\int }_0^{t}d{t}^{\prime }\lambda \varphi \left(t^{\prime }\right).\end{array}$$

(64)

The response function is

$$\varphi \left(t\right)=-{\displaystyle \frac{1}{i\hslash }}{\langle [{\mathcal{H}}_{ext},e^{{\displaystyle \frac{i}{\hslash }}{\mathcal{H}}_{0}t}\mathcal{O}{e}^{-{\displaystyle \frac{i}{\hslash }}{\mathcal{H}}_{0}t}]\rangle }_0=-{\displaystyle \frac{i}{\hslash }}{\langle [U_0^{†}\left(t\right)\mathcal{O}{U}_0\left(t\right),{\mathcal{H}}_{ext}]\rangle }_0,$$

(65)

where we used the subscript 0 in $U_0\left(t\right)$ to distinguish the evolution operator related to the equilibrium dynamics ${\mathcal{H}}_0$ from $U\left(t\right)$ of the total dynamic $\mathcal{H}$: $U_0^{†}\left(t\right)\mathcal{O}{U}_0\left(t\right)\ne {\mathcal{O}}_t=U^{†}\left(t\right)\mathcal{O}U\left(t\right)$. In conclusion, by setting $\hslash =1$, we obtain

$${\langle \mathcal{O}\rangle }_t={\langle \mathcal{O}\rangle }_0-i\lambda {\int }_0^t{\langle [U_0^{†}\left(t\right)\mathcal{O}{U}_0\left(t\right),{\mathcal{H}}_{ext}]\rangle }_{0}d{t}^{\prime }.$$

(66)

In Figure 2, we compare the results with different values of $\lambda$, ${\rho }_0$, time interval and observable $\mathcal{O}$, for the expressions of $\mathcal{O}={\sigma }_y$ given by the linear theory, by the Heisenberg picture and by the Dissipation Function formalism. Qualitatively very similar results (as $\lambda$ varies) were found for other self-adjoint observables, and different initial density operators ${\rho }_0$ (with $[{\mathcal{H}}_0,{\rho }_0]=0)$.

For small $\lambda$ the linear approximation is accurate, and the difference between linear and ${\Omega }^0$ formulas is indistinguishable to the eye. As $\lambda$ increases, the linear approximation becomes increasingly worse, as expected, while the ${\Omega }^0-$response continues to provide results consistent with those of the Heisenberg picture. Moreover, the linear response theory becomes less accurate as time increases, whereas the expression in (26) does not exhibit this flaw. This illustrates the limitations of linear theories even in simple systems. However, for small perturbations, linear response can handle much more complex dynamics than these ones. It remains to be seen how useful the Dissipation Function can be in such situations.

6. Properties of Dissipation and Time Dependent Perturbations

It is clear that the Dissipation Operator plays a very precise role in the temporal evolution of the density operator ${\rho }_t$:

$${\displaystyle \frac{\partial {\rho }_0}{\partial t}}={\rho }_0{\Omega }^0.$$

(67)

However, we can consider any time $t^{*}$ as the initial instant of the dynamics and, consequently, define

$${\Omega }^{t^{*}}={\displaystyle \frac{1}{i\hslash }}{\rho }_{t^{*}}^{-1}[\mathcal{H},{\rho }_{t^{*}}]$$

(68)

which allows us to write the Von Neumann equation (11) as

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial {\rho }_t}{\partial t}}={\rho }_t{\Omega }^t;\rho \left(0\right)={\rho }_0.\end{array}$$

(69)

In general, the formal solution of an equation of this form is expressed using the anti-time ordering operator ${\mathcal{T}}_R$ [52],

$${\rho }_t={\rho }_0{\mathcal{T}}_R\left[e^{{\int }_0^t{\Omega }^{s}ds}\right].$$

(70)

To make use of eq. (70), we need to express ${\Omega }^s$ without explicitly relying on ${\rho }_s$, which is contained in the definition of the Dissipation Operator. This can be done for constant Hamiltonian dynamics.

Proposition 4.

If the Hamiltonian $\mathcal{H}$ is time-independent, ${\Omega }^s$ coincides with the initial Dissipation Operator ${\Omega }^0$ evolved backward in time,

$${\Omega }^s=U\left(s\right){\Omega }^0{U}^{†}\left(s\right).$$

(71)

Proof. 

By definition we have

$${\Omega }^s={\displaystyle \frac{1}{i\hslash }}{\rho }_s^{-1}[\mathcal{H},{\rho }_s].$$

(72)

We can express ${\rho }_s$ using time propagators,

$$\begin{array}{c}\hfill {\Omega }^s={\displaystyle \frac{1}{i\hslash }}{\left(U\left(s\right){\rho }_0{U}^{†}\left(s\right)\right)}^{-1}[\mathcal{H},U\left(s\right){\rho }_0{U}^{†}\left(s\right)]\\ \hfill ={\displaystyle \frac{1}{i\hslash }}{\left(U^{†}\left(s\right)\right)}^{-1}{\rho }_0^{-1}{U}^{-1}\left(s\right)[\mathcal{H},U\left(s\right){\rho }_0{U}^{†}\left(s\right)]\\ \hfill ={\displaystyle \frac{1}{i\hslash }}U\left(s\right){\rho }_0^{-1}{U}^{†}\left(s\right)\left(\mathcal{H}U\left(s\right){\rho }_0{U}^{†}\left(s\right)-U\left(s\right){\rho }_0{U}^{†}\left(s\right)\mathcal{H}\right)\\ \hfill ={\displaystyle \frac{1}{i\hslash }}U\left(s\right){\rho }_0^{-1}{U}^{†}\left(s\right)\mathcal{H}U\left(s\right){\rho }_0{U}^{†}\left(s\right)-{\displaystyle \frac{1}{i\hslash }}U\left(s\right){\rho }_0^{-1}{U}^{†}\left(s\right)U\left(s\right){\rho }_0{U}^{†}\left(s\right)\mathcal{H};\end{array}$$

(73)

this expression can be simplified noting that

$$U\left(s\right){\rho }_0^{-1}{U}^{†}\left(s\right)U\left(s\right){\rho }_0{U}^{†}\left(s\right)=Id,\mathrm{and}[U\left(t\right),\mathcal{H}]=0.$$

(74)

We can then write:

$$\begin{array}{c}\hfill {\Omega }^s={\displaystyle \frac{1}{i\hslash }}U\left(s\right){\rho }_0^{-1}\mathcal{H}{\rho }_0{U}^{†}\left(s\right)-{\displaystyle \frac{1}{i\hslash }}U\left(s\right)U^{†}\left(s\right)\mathcal{H}U\left(s\right)U^{†}\left(s\right)\\ \hfill ={\displaystyle \frac{1}{i\hslash }}U\left(s\right)\left({\rho }_0^{-1}\mathcal{H}{\rho }_0-U^{†}\left(s\right)\mathcal{H}U\left(s\right)\right)U^{†}\left(s\right)\\ \hfill ={\displaystyle \frac{1}{i\hslash }}U\left(s\right)\left({\rho }_0^{-1}\mathcal{H}{\rho }_0-\mathcal{H}\right)U^{†}\left(s\right)\\ \hfill =U\left(s\right)\left({\displaystyle \frac{1}{i\hslash }}{\rho }_0^{-1}[\mathcal{H},{\rho }_0]\right)U^{†}\left(s\right)=U\left(s\right){\Omega }^0{U}^{†}\left(s\right).\end{array}$$

(75)

□

Interestingly, the Dissipation Operator does not evolve in time like standard observables in the Heisenberg representation, $U^{†}\left(t\right)\mathcal{O}U\left(t\right)$. This might seem problematic: we want to consider $\Omega$ as a physical observable, but it evolves like a probability i.e. in reverse time: ${\Omega }^t=U\left(t\right){\Omega }^0{U}^{†}\left(t\right)$. This is not due to the order of its operators in the definition. One could try to redefine the Dissipation Operator by swapping the order of $\mathcal{H},{\rho }_0$ and ${\rho }_0^{-1}$, but it would still evolve in reversed time. The origin of this type of evolution lies in the operators involved: $\mathcal{H}$ and $\rho$. The Hamiltonian operator $\mathcal{H}$ (in this simple case) commutes with time propagators, while $\rho$ and ${\rho }^{-1}$ evolve at reversed times.

If we want a Dissipation Operator that evolves like any other physical observable, we should avoid the use of $\rho$. At the moment, this is not of our concern. Indeed, upon further analysis, this result does not seem to be a disadvantage at all; on the contrary, it looks more a strength. It provides consistent results, it is equivalent with the Heisenberg and Schrödinger pictures for usual observables. Using ${\Omega }^0$, we found

$${\displaystyle \frac{d}{ds}}{\langle \mathcal{O}\rangle }_s=Tr\left({\rho }_0{\Omega }^0{\mathcal{O}}_s\right)=Tr\left(\mathcal{O}{\displaystyle \frac{\partial {\rho }_s}{\partial s}}\right)$$

(76)

which conceptually is reminiscent of the Heisenberg picture, with evolving observables, and of the Schrödinger picture with evolving probabilities. In turn, the equivalence of the two is obtained thanks to the reverse-time evolution of ${\Omega }^0$:

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{d}{ds}}{\langle \mathcal{O}\rangle }_s\right)}^S=Tr\left(\mathcal{O}{\displaystyle \frac{\partial {\rho }_s}{\partial s}}\right)=Tr\left(\mathcal{O}{\rho }_s{\Omega }^s\right)=Tr\left(\mathcal{O}U\left(s\right){\rho }_0{U}^{†}\left(s\right)U\left(s\right){\Omega }^0{U}^{†}\left(s\right)\right)\\ \hfill =Tr\left(U^{†}\left(s\right)\mathcal{O}U\left(s\right){\rho }_0{\Omega }^0\right)=Tr\left({\mathcal{O}}_s{\rho }_0{\Omega }^0\right)=Tr\left({\rho }_0{\Omega }^0{\mathcal{O}}_s\right)={\left({\displaystyle \frac{d}{ds}}{\langle \mathcal{O}\rangle }_s\right)}^H.\end{array}$$

(77)

where the superscript S in the first term and the superscript H in the last one respectively stand for Schrödinger and Heisenberg. Moreover, eq. (71) yields another consistency result:

$${\langle {\Omega }^t\rangle }_t=Tr\left(U\left(t\right){\Omega }^0{U}^{†}\left(t\right)U\left(t\right){\rho }_0{U}^{†}\left(t\right)\right)={\langle {\Omega }^0\rangle }_0=0$$

(78)

which correctly means that every moment can be considered as the initial one. The property of zero-mean, stated in proposition 1, is propagated in time thanks to eq. (71). Finally, in applications of classical dynamical systems, a commonly used quantity is ${\langle {\Omega }^0\rangle }_t$ [35], but here ${\langle {\Omega }^0\rangle }_t\ne {\langle {\Omega }^t\rangle }_0$, unlike what happens to the usual observables.

Equation (70) shows that the solution to the Von Neumann problem (11) can be expressed in several ways: not only in the usual form, eq. (14), but also as in eq. (70); or equivalently, using the adjoint operator ${\Omega }^{†}$ and the time-ordering operator ${\mathcal{T}}_L$, as

$${\rho }_t={\mathcal{T}}_L\left[e^{{\int }_0^t{\left({\Omega }^{†}\right)}^{s}ds}\right]{\rho }_0.$$

(79)

If the solution is unique, expressions (70) and (79) coincide with the usual form (14); they are just different ways of writing the same thing.

Applying our theory to complex, particuarly open quantum systems for which traditional time propagators do not exist, could prove useful. Let us investigate this possibility. Take a time-dependent Hamiltonian $\mathcal{H}\left(t\right)$. The usual time evolution is expressed by [52]:

$$\begin{array}{c}\hfill {\mathcal{U}}^{†}\left(t\right)={\mathcal{T}}_R\left[exp\left(+{\displaystyle \frac{i}{\hslash }}{\int }_0^t\mathcal{H}\left(s\right)ds\right)\right];\mathcal{U}\left(t\right)={\mathcal{T}}_L\left[exp\left(-{\displaystyle \frac{i}{\hslash }}{\int }_0^t\mathcal{H}\left(s\right)ds\right)\right];\\ \hfill {\rho }_t=\mathcal{U}\left(t\right){\rho }_0{\mathcal{U}}^{†}\left(t\right).\phantom{aaaaaaaaaaaaaaaaaaaaaaaaa}\phantom{bbbbbbbbbbbbbbbbbbbbbbbbb}\end{array}$$

(80)

In terms of dissipation operators, one can instead follow different approaches. We choose here self-adjoint dissipation operators, introducing the operator ${\omega }^0\equiv {\hslash }^{-1}{\rho }_0^{-1}\mathcal{H}{\rho }_0$, and then we symmetrize it to make it self-adjoint:

$${\tilde{\omega }}^0={\displaystyle \frac{1}{2\hslash }}\left({\rho }_0^{-1}\mathcal{H}{\rho }_0+{\rho }_0\mathcal{H}{\rho }_0^{-1}\right).$$

(81)

We have:

$$\begin{array}{c}\hfill \overline{{\Omega }^0}={\displaystyle \frac{1}{2}}\left({\Omega }^0-{\left({\Omega }^0\right)}^{†}\right)={\displaystyle \frac{1}{2i\hslash }}\left({\rho }_0^{-1}[\mathcal{H},{\rho }_0]-[\mathcal{H},{\rho }_0]{\rho }_0^{-1}\right)\\ \hfill ={\displaystyle \frac{1}{2i\hslash }}\left({\rho }_0^{-1}\mathcal{H}{\rho }_0-2\mathcal{H}+{\rho }_0\mathcal{H}{\rho }_0^{-1}\right)={\displaystyle \frac{1}{i}}\left({\tilde{\omega }}^0-{\displaystyle \frac{1}{\hslash }}\mathcal{H}\right).\end{array}$$

(82)

Now, we can rewrite the Hamiltonian as $\mathcal{H}\left(s\right)=\hslash ({\tilde{\omega }}^s-i\overline{{\Omega }^s})$ and define the operator ${\widehat{\Omega }}^s\equiv {\tilde{\omega }}^s-i\overline{{\Omega }^s}=\mathcal{H}/\hslash$. This allows us to write:

$$\widehat{\mathcal{U}}\left(t\right)\equiv {\mathcal{T}}_L{e}^{-i{\widehat{\Omega }}_{0,t}},{\widehat{\mathcal{U}}}^{†}\left(t\right)={\mathcal{T}}_R{e}^{+i{\widehat{\Omega }}_{0,t}}={\widehat{\mathcal{U}}}^{-1}\left(t\right);{\mathcal{O}}_t={\widehat{\mathcal{U}}}^{†}\left(t\right)\mathcal{O}\widehat{\mathcal{U}}\left(t\right).$$

(83)

Clearly, for constant Hamiltonians, these ordered exponentials turn mere expoentials, because the two time-dependent operators $\overline{{\Omega }^t}$ and ${\tilde{\omega }}^t$ reduce to the constant operator $\mathcal{H}$. While this is simply a rewriting of the usual operators, it suggests an interesting possibility: it may be possible to construct time propagators for dissipative dynamics by extending the Dissipation Operator formalism to open quantum systems.

7. The Dissipation Operator for Open Quantum Systems

The Lindblad equation [9]:

$$\dot{\rho }={\displaystyle \frac{1}{i\hslash }}[\mathcal{H},\rho ]+\sum _{\alpha }\left(L_{\alpha }\rho L_{\alpha }^{†}-{\displaystyle \frac{1}{2}}\{L_{\alpha }^{†}{L}_{\alpha },\rho \}\right),$$

(84)

with the commutator accounting for the unitary non dissipative evolution and the $L_{\alpha }$ operators for the rest, are commonly used to describe open quantum systems. In particualr, $\mathcal{H}$ is the effective Hamiltonian operator, which may differ from the Hamiltonian of the isolated system due to the interaction with the environment. The second term represents the dissipative part of the evolution, with the Lindblad operators $L_{\alpha }$ describing the interaction between system and external environment. To apply our method to these dynamics, we extend the Dissipation Operator incorporating the dissipative effects introduced by the Lindblad operators. Again, there are several equivalent ways to do it; we adopt the following:

$$D_{\alpha }\equiv {\rho }^{-1}{L}_{\alpha }\rho L_{\alpha }^{†}-L_{\alpha }^{†}{L}_{\alpha }\Rightarrow {\Omega }_L^0\equiv {\Omega }^0+\sum _{\alpha }{D}_{\alpha }^0.$$

(85)

Here, ${\Omega }^0$ concerns as above the Hamiltonian part of the evolution, $D_{\alpha }$ takes into account the dissipative part:

$$D_{\alpha }^{†}={\left({\rho }^{-1}{L}_{\alpha }\rho L_{\alpha }^{†}\right)}^{†}-{\left(L_{\alpha }^{†}{L}_{\alpha }\right)}^{†}={\left(\rho L_{\alpha }^{†}\right)}^{†}{\left({\rho }^{-1}{L}_{\alpha }\right)}^{†}-{\left(L_{\alpha }\right)}^{†}{\left(L_{\alpha }^{†}\right)}^{†}=L_{\alpha }\rho L_{\alpha }^{†}{\rho }^{-1}-L_{\alpha }^{†}{L}_{\alpha }$$

(86)

and linearity implies ${\left({\sum }_{\alpha }{D}_{\alpha }\right)}^{†}={\sum }_{\alpha }{D}_{\alpha }^{†}$. Then, the Lindblad equation can be rewritten as

$${\displaystyle \frac{\partial \rho }{\partial t}}={\displaystyle \frac{1}{2}}\left(\rho {\Omega }_L^t+{\left({\Omega }_L^t\right)}^{†}\rho \right)$$

(87)

Indeed, we can write:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial \rho }{\partial t}}& =& {\displaystyle \frac{1}{2}}\left(\rho {\Omega }_L^t+{\left({\Omega }_L^t\right)}^{†}\rho \right)\hfill \\ & =& {\displaystyle \frac{1}{2}}\left(\rho {\Omega }^t+{\left({\Omega }^t\right)}^{†}\rho \right)+{\displaystyle \frac{1}{2}}\sum _{\alpha }\left(L_{\alpha }\rho L_{\alpha }^{†}-\rho L_{\alpha }^{†}{L}_{\alpha }+L_{\alpha }\rho L_{\alpha }^{†}-L_{\alpha }^{†}{L}_{\alpha }\rho \right)\hfill \\ & =& {\displaystyle \frac{1}{i\hslash }}[\mathcal{H},\rho ]+\sum _{\alpha }\left(L_{\alpha }\rho L_{\alpha }^{†}-{\displaystyle \frac{1}{2}}\{L_{\alpha }^{†}{L}_{\alpha },\rho \}\right).\hfill \end{array}$$

(88)

We note that, for Hamiltonian dynamics, one can express the Von Neumann equation with a linear combination of $\Omega$ and ${\Omega }^{†}$:

$${\displaystyle \frac{\partial \rho }{\partial t}}=a{(\Omega )}^{†}\rho +b\rho \Omega \mathrm{with}a,b\in \mathbb{R}s.t.a+b=1.$$

(89)

Applying definition (85) to the Lindblad equation, we take $a=b=1/2$, consistently with eq. (87). Various practical challanges now arise to explicitly solve the equations of interest. Therefore, to illustrate how our approach works, here we introduce some simplifying hypotheses. In particular, we take:

$${\mathcal{O}}_t=U_L^{\times }\left(t\right)\mathcal{O}{U}_L\left(t\right)=e^{{\displaystyle \frac{i}{\hslash }}(\mathcal{H}+C^{\times })t}\mathcal{O}{e}^{-{\displaystyle \frac{i}{\hslash }}(\mathcal{H}+C)t}.$$

(90)

where the symbol ’×’ indicates that the two evolution operators need not necessarily be adjoints of each other. This hypotesis is suggested by the extension of $\widehat{\Omega }$ in the form:

$${\widehat{\Omega }}_L^s=\left({\tilde{\omega }}^s-i\overline{{\Omega }^s}\right)+\left(g^s-k^s\right)={\displaystyle \frac{1}{\hslash }}\left(\mathcal{H}+C\right),$$

(91)

where $g(L_{\alpha },L_{\alpha }^{†},{\rho }_s)$ and $k(L_{\alpha },L_{\alpha }^{†},{\rho }_s)$ represent the dissipative part of ${\tilde{\omega }}_L$ and ${\overline{\Omega }}_L$ respectively. These are time-dependent operators. The term $k^s$ is determined by the definition of ${\overline{\Omega }}_L$, but the term $g^s$ has no specific constraints. This allows us to define $g^s$ consistently with the correct results. It is reasonable to express these terms through a general operator C, precisely because of this freedom with respect to g. Additionally, assuming that C is time-independent is not unrealistic. Just as the two time-dependent operators $\overline{{\Omega }^t}$ and ${\tilde{\omega }}^t$ reduce to the constant Hamiltonian $\mathcal{H}$, the same could happen for $g^t$ and $k^t$, especially since all the Lindblad operators in equation (84) are time-independent. We stress that in eq. (90), we are not assuming $U_L$ to be a "real" time propagator, in the sense that it would extend the Heisenberg representation, ${\langle \mathcal{O}\rangle }_t=Tr\left(U_L^{\times }\left(t\right)\mathcal{O}{U}_L\left(t\right){\rho }_0\right)$. Instead, we are stating that these operators behave like evolution operators when combined with the generalization of expression (26). This generalization can be derived exactly by repeating the calculations already performed earlier, starting from $d_t{\langle \mathcal{O}\rangle }_t=Tr\left({\mathcal{O}}_t^{*}{\partial }_t{\rho }_0\right)$, where the only difference is represented by eq. (87). This leads to:

$${\langle \mathcal{O}\rangle }_t={\langle \mathcal{O}\rangle }_0+{\displaystyle \frac{1}{2}}{\int }_0^t\left({\langle {\Omega }_L^0{\mathcal{O}}_s^{*}\rangle }_0+{\langle {\mathcal{O}}_s^{*}{\left({\Omega }_L^0\right)}^{†}\rangle }_0\right)ds.$$

(92)

Let us apply expression (92) to a particular Lindblad equation, whose analytical solution is known, suitably defining C and $C^{\times }$, and consequently $U_L\left(t\right)$ and $U_L^{\times }\left(t\right)$. Consider the equation,

$${\displaystyle \frac{\partial \rho }{\partial t}}=-i[\omega {\sigma }_z,\rho ]+\gamma \left({\sigma }^{-}\rho {\sigma }^{+}-{\displaystyle \frac{1}{2}}\{{\sigma }^{+}{\sigma }^{-},\rho \}\right);{\sigma }^{\pm }\equiv {\displaystyle \frac{1}{2}}({\sigma }_x\pm i{\sigma }_y),$$

(93)

which describes a decoherence phenomenon of a qubit open system. This equation is already in the form of (84); we just need to highlight the correspondence with the general notation:

$$\alpha =1;\mathcal{H}=\omega {\sigma }_z;L=\sqrt{\gamma }{\sigma }^{-};L^{†}=\sqrt{\gamma }{\sigma }^{+}.$$

(94)

After writing equation (93) in its single components, the analytical solution can be easily calculated:

$$\rho \left(t\right)=\left[\begin{array}{cc}{\rho }_{00}^0{e}^{-\gamma t}& {\rho }_{01}^0{e}^{-(0.5\gamma +2i\omega )t}\\ {\rho }_{10}^0{e}^{-(0.5\gamma -2i\omega )t}& {\rho }_{11}^0+{\rho }_{00}^0(1-e^{-\gamma t})\end{array}\right]$$

(95)

and then the evolution of observables, ${\langle \mathcal{O}\rangle }_t=Tr\left(\rho \left(t\right)\mathcal{O}\right)$ immediately follows. Comparing this expression with eq. (92) we find the form of C and $C^{\times }$ for this particular case. Consider first diagonal observables:

$$\mathcal{O}=\left[\begin{array}{cc}a& 0\\ 0& b\end{array}\right],a,b\in \mathbb{R}.$$

(96)

Introducing

$$C\equiv -i\gamma {\sigma }^{+}{\sigma }^{-},C^{\times }\equiv +i\gamma {\sigma }^{-}{\sigma }^{+},$$

(97)

formula (92) gives the known analytical solution (95), for any initial condition $\rho \left(0\right)$ and for all real frequencies $\omega ,\gamma$. In Figure 3, this is shown numerically for the projectors ${\pi }_0$ ($a=1,b=0$) and ${\pi }_1$ ($a=0,b=1$), from which all diagonal observables can be obtained. It is important to note that $C^{\times }$ differs from the adjoint of C, $C^{†}=i\gamma {\sigma }^{+}{\sigma }^{-}$, and that the two operators $U_L$ and $U_L^{\times }$ do not form a unitary group. This is consistent with the dissipative nature of the dynamics.

Consider now the following class of observables:

$$\mathcal{O}=\left[\begin{array}{cc}0& c+id\\ c-id& 0\end{array}\right];c,d\in \mathbb{R}.$$

(98)

that we call diag2 operators. We define the evolution of the observables for this second class as:

$${\mathcal{O}}_t^{*}=e^{{\displaystyle \frac{i}{\hslash }}(\mathcal{H}+M^{\times })t}\mathcal{O}{e}^{-{\displaystyle \frac{i}{\hslash }}(\mathcal{H}+M)t};M=-{\displaystyle \frac{i}{4}}\{{\sigma }^{-},{\sigma }^{+}\}=-M^{\times }.$$

(99)

When combined with expression (92), this provides the exact response for all observables of the form (98), with M, $M^{†}=M^{\times }$. In Figure 4, we numerically illustrate this fact for two cases.

Any self-adjoint observable $\mathcal{O}$ in a two-dimensional Hilbert space can be expressed as a linear combination of the two kinds of operators defined by equations (96) and (98). For the specific case of equation (93), the operators $C,C^{\times },M\mathrm{and}{M}^{\times }$ all commute with the Hamiltonian operator $\mathcal{H}$. In conclusion, for this particular Lindblad equation (93), the expression (92) gives exact results when the observables evolve according to:

$${\mathcal{O}}_t^{*}=e^{{\displaystyle \frac{i}{\hslash }}{C}^{\times }t}{\mathcal{O}}_H^1\left(t\right)e^{-{\displaystyle \frac{i}{\hslash }}Ct}+e^{{\displaystyle \frac{i}{\hslash }}{M}^{\times }t}{\mathcal{O}}_H^2\left(t\right)e^{-{\displaystyle \frac{i}{\hslash }}Mt},$$

(100)

where a simple rule applies:

$$M=-{\displaystyle \frac{i}{4}}\{L_{\alpha },L_{\alpha }^{†}\}=-M^{\times };C=-i{L}_{\alpha }^{†}{L}_{\alpha };C^{\times }=i{L}_{\alpha }{L}_{\alpha }^{†}$$

(101)

with the subscript ’H’ indicating evolution in the Heisenberg picture.

How generally rule (101) applies is object of investigation. However, this is just one way to obtain exact results using the Dissipation Operator. Exact results can be obtained equivalently by evolving the observables in the Heisenberg picture and assuming the following time dependence for the Dissipation Operator:

$${\Omega }_L^0\left(t\right)=e^{-\gamma t}{\Omega }_L^0\mathrm{for} \mathrm{diagonal} \mathrm{operators};{\Omega }_L^0\left(t\right)=e^{-{\displaystyle \frac{1}{2}}\gamma t}{\Omega }_L^0\mathrm{for} \mathrm{diag} 2 \mathrm{operators}$$

(102)

Lindblad equations, while widely used, have certain limitations. One issue is that the positivity of the density matrix $\rho$ can be violated, particularly in numerical simulations (1.5.2 in text [9]). Moreover, their validity is restricted to linearly dissipative interactions, thus choosing $L_{\alpha }$ operators can be an ambiguous or approximate process [53,54]. Furthermore, when many energy states are involved, calculating the solution becomes very computationally intensive [55,56]. We conclude observing that Lindblad equations are intimately linked to the concept of quantum entropy [57,58], which is affected by dissipation, of course. Therefore, combining them with the response theory based on the Dissipation Operator may prove useful.

8. Concluding Remarks

In this paper, we have presented a first approach to the formulation of a quantum exact response theory as an extension of the classical theory based on the Dissipation Function. Our future goal is to further develop this approach for studying perturbations in complex open quantum systems. To this end, we have started from constant Hamiltonian dynamics, which suggest at least two quite natural quantum versions of the Dissipation Function, that are equivalent for the cases considered here.

The first quantum Dissipation Operator has been defined by (18), and the corresponding quantum exact response expression has been given by (26), in perfect analogy with the classical formula (4). For constant Hamiltonian dynamics, the correctness of the results can be tested, comparing the expectation values of observables ${\langle \mathcal{O}\rangle }_t$ with the corresponding Schrödinger and Heisenberg pictures.

We then noted that, although it is not strictly required for open dynamics, it is still interesting to use self-adjoint operators. Therefore, we proposed a self-adjoint Dissipation Operator, defined by eq. (44), that is equivalent to the first. In our derivation, this requires the introduction of an anti-self-adjoint operator. In this formulation, symmetric quantum correlations appear, with properties quite similar to the classical case.

We applied these expressions to qubit systems and obtained numerically correct results. A comparison with linear response theory revealed that as the perturbation value and time increase, the linear response naturally worsens, whereas the expressions (26) and (44) continue to yield exact results. We also examined the role of the Dissipation Operator in the time evolution of the density matrix $\rho$. In particular, we found that for constant Hamiltonian dynamics, the operator ${\Omega }^t$ evolves in reverse time, as described in equation (71). This is not what observables do, but it is consistent with the classical theory, where the Dissipation Function is at once an observable and the generator of the evolution of probabilities. As far as response is concerned, this poses no problem.

Thinking preferable to keep the Dissipation Operator ${\Omega }^0$ fixed in time, as in the classical theory, we introduced time propagators based on the dissipation operators, as given in eq. (83). We did it consistently, rewriting the Hamiltonian in terms of dissipation operators, which may seem trivial but suggests an interesting insight: by extending this approach, it may be possible to properly treat observables of dissipative dynamics. Building on this idea, we extended the method to Lindblad equations, incorporating the dissipative part (described by $L_{\alpha }$ and $L_{\alpha }^{†}$) in the definition of the Dissipation Operator, leading to the quantum exact response expression (92). We applied this approach to a specific Lindblad equation, comparing the results obtained through our new expression with the known analytical solution. With the aid of assumptions suggested by both previous studies and by inspection, we obtained exact results with expression (92).

While this approach is still in its infancy, and it is open to different approaches, it offers an original perspective on quantum response theory and open quantum systems. Its main strength lies in its analogy with the classical exact response theory, which is robust and general [31,33,35,36,37]. Moreover, the computational advantages recently offered by classical theory across various applications [37,59], along with the development of new methods for investigating the Dissipation Function [41], provide further compelling evidence of the potential held by its proposed quantum analogue. A promising avenue for future research is to investigate whether $\Omega$ can be associated with entropy production rates, and with the mathematical tools of large deviation theory, such as the Kulback-Leibler divergence. Future studies intend to clarify this, also in connection with the application to different equations for open quantum systems. Finally, extension to infinite-dimensional Hilbert spaces should be considered.