Introduction —

The dynamics of open quantum systems arise from the interplay between coherent evolution and environmental interactions [1], [2]. Beyond inducing dissipation and decoherence, these interactions can renormalize the energy structure of the system, giving rise to phenomena such as the Lamb [3], [4] and AC Stark shifts [5]. The Lamb shift, a prototypical renormalization effect caused by electromagnetic vacuum fluctuations, exemplifies how environmental couplings modify intrinsic system properties, with far-reaching implications for both theory and applications.

Within the Lindblad framework, such shifts emerge naturally from tracing out environmental degrees of freedom under weak-coupling and Markovian assumptions [1]. Extending this understanding to strong-coupling and non-Markovian regimes, however, remains an open challenge [6]–[12].

A powerful approach to non-Markovian dynamics is provided by time-convolutionless (TCL) master equations [1], [13], which describe exact reduced dynamics through a time-local generator \(\mathcal{L}_{t}\). The generator can be written in a generalized Lindblad form [14], [15]: \[\begin{align} \label{eq:no95gauge} \mathcal{L}_t[X] = -i [K(t), X] + \mathcal{D}_t[X], \end{align}\tag{1}\] where \(K(t)\) denotes the effective Hamiltonian encoding coherent energy shifts, and \(\mathcal{D}_t[X]\) represents the dissipative part, which may include negative rates —a possible signature of non-Markovian behavior [16], [17].

However, the decomposition of \(\mathcal{L}_t\) into its coherent and dissipative components is not unique due to the intrinsic symmetry of the generator [1], [18]. Resolving this ambiguity is essential for a consistent quantification of coherent energy shifts beyond the Markovian regime.

Hayden and Sorce recently proposed a foundational solution to this issue [19], defining the effective Hamiltonian by minimizing the dissipative contribution. This “minimal dissipation” approach isolates the coherent component of the dynamics and yields a dissipator with traceless Lindblad operators, recovering the standard weak-coupling limit where Hamiltonian renormalization is well understood. The physical relevance of this effective Hamiltonian has been experimentally confirmed in trapped-ion platforms operating in a highly non-Markovian regime Colla2024exp?, where the resulting level shifts are typically time-dependent.

Beyond its foundational importance, the effective Hamiltonian also provides a key conceptual link to quantum thermodynamics in the strong-coupling and non-Markovian domains [20]. In quantum thermal machines, it captures energy renormalization effects arising from system–bath interactions, which crucially influence performance metrics such as efficiency and work extraction [21]. For time-dependent effective renormalization, \(K(t)\) can describe work exchange between system and bath, revealing novel thermodynamic behaviors and challenging conventional energy-flow paradigms [20], [22], colla2024FA?.

Despite its significance, determining \(K(t)\) is often analytically or numerically demanding, especially in strongly coupled or high-dimensional systems, as time-local generators are rarely accessible in closed form. Exact expressions exist only for specific integrable models such as the Jaynes–Cummings [23], Fano–Anderson [21], [22], colla2024FA?, and Caldeira–Leggett [24] Hamiltonians, linearly coupled Gaussian systems [25], and pure dephasing models [26], [27]. In general, \(\mathcal{L}_t\) may be obtained numerically using methods such as hierarchical equations of motion (HEOM) [28] or pseudomode approaches [29], [30], or experimentally by quantum process tomography. In perturbative regimes, TCL expansions in the coupling strength provide an effective alternative.

In this Letter, we derive an operational expression for the effective Hamiltonian \(K(t)\) directly from the time-local generator \(\mathcal{L}_t\), enabling straightforward extraction of the coherent part of the dynamics once \(\mathcal{L}_t\) is known. We further present a recursive perturbative expansion for \(K(t)\) in powers of the system–bath coupling, expressed in terms of interaction operators and bath correlation functions. This construction allows systematic analysis of Hamiltonian renormalization across coupling regimes. We apply the method to paradigmatic spin models, illustrating how environmental characteristics affect the structure of the renormalized Hamiltonian.

Effective Hamiltonian –

The procedure recently proposed by Hayden and Sorce to uniquely fix the effective Hamiltonian is based on projecting the generator \(\mathcal{L}\) into the subspace of superoperators that act as commutators with Hermitian operators, i.e. \(\mathcal{H} = -i[H, \cdot]\), where \(H = H^\dagger\) [19]. This projection is performed using a scalar product defined by:

\[\begin{align} \label{eq:sp95Haar} \langle \mathcal{L}_1, \mathcal{L}_2 \rangle := \overline{\bra{\psi} \overline{\mathcal{L}_1[\ket{\phi}\bra{\phi}] \mathcal{L}_2[\ket{\phi}\bra{\phi}] } \ket{\psi}}, \end{align}\tag{2}\] where the overline denotes averaging over Haar-random pure states \(\ket{\phi}\) and \(\ket{\psi}\). This approach minimizes the dissipator assuming maximum uncertainty over input and output states (i.e. uniformly sampling over pure states), yielding a canonical effective Hamiltonian \(K\) for any linear, Hermiticity-preserving, trace-annihilating generator \(\mathcal{L}\in \mathfrak{htp}(\mathscr{H})\) (where \(\mathscr{H}\) denotes the Hilbert space of the reduced system).

The physical interpretation of this procedure is better understood via the average fidelity. The fidelity between a noiseless von Neumann generator \(\mathcal{H} = -i[H, \cdot]\), associated with some Hamiltonian \(H\), and the generator of a noisy channel described by \(\mathcal{L}\) is defined as: \[\begin{align} \overline{\mathcal{F}}(\mathcal{H}, \mathcal{L}) = \int d\psi\, \mathrm{Tr}\left\{\mathcal{H}[\ket{\psi}\bra{\psi}]\,\mathcal{L}[\ket{\psi}\bra{\psi}]\right\}. \end{align}\] This quantity measures how “faithfully” the noisy dynamics generated by \(\mathcal{L}\) approximates that of the noiseless \(\mathcal{H}\) when applied to randomly chosen pure states, and is connected to the scalar product 2 between the generators simply through a scaling with the Hilbert space dimensions \(d = \dim(\mathscr{H})\), i.e. \(\braket{\mathcal{H},\mathcal{L}} = \frac{1}{d}\overline{\mathcal{F}}(\mathcal{H}, \mathcal{L})\).

Expanding \(H\) in a basis of traceless Hermitian operators \(\{H_j\}_j\), orthonormal under the Hilbert-Schmidt product, the effective Hamiltonian is expressed as (see Appendix 1): \[\begin{align} \label{eq:Ks-fidelity} K = \frac{d+1}{2} \sum_j \overline{\mathcal{F}}(\mathcal{H}_j, \mathcal{L}) H_j, \end{align}\tag{3}\] where \(\mathcal{H}_{j} =-i [H_{j},\cdot]\) This expression for \(K\) establishes a direct connection between the average fidelity and the effective Hamiltonian. It reveals that the effective Hamiltonian is a weighted sum of basis operators \(H_j\), with weights determined by the average fidelity between \(\mathcal{L}\) and the generators corresponding to each basis operator. It quantifies how well different components of the dynamics preserve quantum information on average and clarifies the equal apriori probability – or maximum uncertainty – principle embedded in the choice of the norm induced by the scalar product 2 .

Equation 3 provides physical insight into \(K\), but it is not practical for calculations. Therefore, we provide a general expression for the effective Hamiltonian associated with the canonical splitting under 2 of any given generator \(\mathcal{L}\in \mathfrak{htp}(\mathscr{H})\), which constitutes our first result.

The above equation can be validated by substituting \(\mathcal{L}\) with a generic canonical generator satisfying minimal dissipation, i.e eq. 1 with traceless Lindblad operators. We further provide a detailed derivation of eq. ?? starting from the definition of the scalar product in eq. 2 (see Appendix 2). Using the operator basis \(F_\alpha = |j\rangle\langle k|\), where \(|j\rangle\) form an orthonormal basis on \(\mathscr{H}\), this becomes: \[\begin{align} \label{eq:K-jk} K = \frac{1}{2i d} \sum_{j,k=1}^d [|j\rangle\langle k|,\,\mathcal{L}[|k\rangle\langle j|]]. \end{align}\tag{4}\] These new expressions are applicable to a wide range of quantum systems and coupling regimes, and offer a more practical approach for computing \(K\) whenever the exact generator is explicitly known. However, in most cases, an explicit form of \(\mathcal{L}\) is only accessible via perturbative treatment. In these cases, the form of the generator is typically given by \(\mathcal{L}[X] = \sum_{ij} \omega_{ij} V_{i} X W_{j}\), with \(V_{i}\) and \(W_{j}\) satisfying Hermiticity-preservation and trace-annihilation conditions. This leads to an effective Hamiltonian of the following form: \[\begin{align} \label{eq:K95rec} K = \frac{1}{2i d} \sum_{ij} \omega_{ij} \left( \mathrm{Tr}(V_i) W_j - \mathrm{Tr}(W_j) V_i \right), \end{align}\tag{5}\] where the result from [19] is recovered when \(\omega_{ij}\) is diagonal (pseudo-Kraus representation). In what follows, we provide a technique to recover the contributions to the effective Hamiltonian at each order of the perturbation expansion.

Systematic perturbation expansion —

The generator \(\mathcal{L}_t\) of an open quantum system typically arises from the microscopic unitary evolution of a closed bipartite system, where the system \(S\) interacts with an environment \(E\) (the irrelevant degrees of freedom that are traced out). The total Hamiltonian governing this evolution is given by: \[\begin{align} H = H_S + H_E + \lambda H_I, \end{align}\] where \(H_S\) and \(H_E\) are the free Hamiltonians of the system and environment, respectively, and \(\lambda\) is the coupling parameter for the interaction \(H_I\). For simplicity, we consider a factorized interaction Hamiltonian, \(H_I = A \otimes B\), where \(A\) and \(B\) act on the system and environment, respectively. Extensions to sums of such terms are straightforward. We further assume the total initial state to factorize as \(\rho_{SE}(0) = \rho_S \otimes \rho_E\), since initial system–environment correlations do not affect the coherent part of the system dynamics, i.e., the effective Hamiltonian \(K\) [31]. Their only effect is the addition of an inhomogeneous term to the homogeneous map describing the evolution of the initially uncorrelated state, while the homogeneous part of the dynamics characterized by the generator in eq. 1 remains unchanged. See Appendix 3 for a detailed discussion.

Finding the exact effective Hamiltonian \(K\) for the reduced system \(S\) is generally challenging without an explicit expression for the time local generator \(\mathcal{L}_{t}\) associated to the reduced dynamics of the system only. Here, we provide a systematic perturbation expansion for \(K\) in powers of the coupling strength \(\lambda\). This constitutes the second main result of this work. Below, we outline the strategy; technical details are provided in [32].

The interaction picture is used to simplify calculations, where the interaction Hamiltonian becomes time-dependent: \(\tilde{H}_t = A(t) \otimes B(t)\). The total unitary evolution can then be written as: \[\begin{align} \mathcal{U}_{t}[\rho_{SE}(0)] = \mathcal{T}\left( e^{-i \int_0^t d\tau (\tilde{H}\tau^L - \tilde{H}\tau^R)} \right) [\rho_{SE}(0)], \end{align}\] where \(\mathcal{T}\) denotes the time ordering and \(\tilde{H}_\tau^{L/R}\) represents the operator acting on \(\rho\) from the left and from the right [33]–[35].

The reduced dynamics of the system is obtained by tracing out the environment degrees of freedom, yielding to \(\Phi_t[\rho_S] = \mathrm{Tr}_E \{\mathcal{U}_t[\rho_S \otimes \rho_E]\}\). Expanding the unitary evolution in a Taylor series, we express \(\Phi_{t}\) as:

\[\begin{align} \Phi_t = 1 + \sum_{n=1}^\infty (-i)^n \lambda^n \mu_n, \end{align}\] where each term \(\mu_n\) depends on time-ordered integrals over system and bath operators. These contributions can be decomposed as: \[\begin{align} \mu_n^k & = \int_0^t d\boldsymbol{\tau}_1^k d\boldsymbol{s}_1^{n-k} D(\boldsymbol{\tau}_1^k, \boldsymbol{s}_1^{n-k}) A^L(\boldsymbol{\tau}_1^k) A^R(\boldsymbol{s}_1^{n-k}), \end{align}\] where we have defined \[\begin{align} {A^{L/R}}(\boldsymbol{t}_1^k)[\,\cdot\,] &:= A^{L/R}_{t_1}\circ{A^{L/R}}_{t_2}\circ...\circ{A^{L/R}}_{t_k} [\,\cdot\,] \; , \end{align}\] and with shorthand integral notation \(\int_{0}^{t}d\boldsymbol{\tau}_{1}^{k} d\boldsymbol{s}_{1}^{n-k}\) indicating integration over \(k\) times \({\tau_1, ..., \tau_k}\) and \((n-k)\) times \({s_1, ..., s_{n-k}}\). Here, \(D(\boldsymbol{\tau}_1^k, \boldsymbol{s}_1^{n-k})\) encodes environmental correlations via \(n\)-point bath correlation functions: \[\begin{align} D(\boldsymbol{\tau}_1^k, \boldsymbol{s}_1^{n-k}) = & \mathrm{Tr} \left\{B^R(\boldsymbol{s}_1^{n-k}) B^L(\boldsymbol{\tau}_1^k) \rho_E \right\} \nonumber\; \theta_{\boldsymbol{\tau}_1^k} \theta_{\boldsymbol{s}_1^{n-k}}, \end{align}\] where \(\theta_{\boldsymbol{\tau}_{1}^{n}}\) ensures time ordering within each set of variables \(\{\tau_{j}\}_{j\in[1,n]}\), i.e. \(\theta_{\boldsymbol{\tau}_{1}^{n}}\) is equal to 1 if \(t_{1}> t_{2}>\dots\tau_{n}\) and zero otherwise.

Following the steps in [32], [36] the generator can be expressed as: \[\begin{align} \mathcal{L}_{t} = \sum_{n=0}^\infty \lambda^n \mathcal{L}_n, \end{align}\] where each term \(\mathcal{L}_n\) depends on combinations of \(\mu_n\) terms and their time derivatives. (For a detailed derivation the perturbative series see [32]) Applying eq. 5 for each order of \(\mathcal{L}_n\), we obtain a perturbative expansion for the effective Hamiltonian: \[\begin{align} K(t) = \sum_{n=0}^{\infty} \lambda^n K_n. \end{align}\] With notation \(\braket{X}_{1/d}:= \mathrm{Tr}\{X\}/d\), and defining the operator \({\mathbb{A}}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right) := \langle A(\boldsymbol{s}_{1}^{n-k})^{\dagger}\rangle_{1/d} A(\boldsymbol{\tau}_{1}^{k})\), the \(n\)th-order contribution to \(K(t)\) is given by:

where the functions \(\mathcal{D}\) are bath \(n\)-ordered cumulants and are linear combinations of products of the \(n\)-point correlation functions (see companion paper [32]) which can be recovered via the following recursion formula: \[\begin{align} \label{eq:Dcal-rec} \mathcal{D}(\boldsymbol{\tau}_{1}^{k},\boldsymbol{s}_{1}^{n-k}) &= \dot{D}(\boldsymbol{\tau}_{1}^{k},\boldsymbol{s}_{1}^{n-k}) \nonumber\\ & \quad -\sum_{l=0}^{k}\sum_{r=0}^{n-k} \mathcal{D}(\boldsymbol{\tau}_{1}^{l},\boldsymbol{s}_{1}^{r}){D}(\boldsymbol{\tau}_{l+1}^{k},\boldsymbol{s}_{r+1}^{n-k}) \end{align}\tag{7}\] where \[\begin{align} \dot{D}(\boldsymbol{\tau}_{1}^{k},\boldsymbol{s}_{1}^{n-k})=D(\boldsymbol{\tau}_{1}^{k},\boldsymbol{s}_{1}^{n-k}) \big( \delta_{\tau_1,t}+\delta_{s_{1},t} \big) \; , \end{align}\] and we have taken the convention \(\int_0^td\tau \delta_{\tau, t} f(\tau) = f(t)\) for any function of \(f\). The time ordering functions \(\theta\) inside each term \(D\) are responsible for the time ordering of all operators and functions in the integrals, and will give rise to a number of different terms that are ordered in “chunks”, similarly to the ordered cumulants of van Kampen [37], [38].

These expressions allow systematic computation of \(K(t)\) at any order in \(\lambda\), revealing how system-bath interactions renormalize energy levels. Notably, eq. 6 provides explicit information about the structure of \(K(t)\) in terms of system operators (\(A\)) and their interaction-picture dynamics and the \(n\)-point correlation functions of the bath.

Spin systems —

To elucidate the structure of the effective Hamiltonian, we apply our results to open spin systems. Without loss of generality, we take \(H_S = \frac{\omega}{2}\sigma_z\). The recursive series for \(K\) here derived reveals whether the interaction with a thermal bath generates only a shift of the energy levels of the free Hamiltonian or also a rotation of its eigenvectors.

Suppose that \(A=\sigma_z\), resulting in a situation of pure dephasing for the spin, where \(A(t)\equiv A\) \(\forall t\). Using 6 , we find (see Appendix 4) that the contribution of the odd terms of the series commute with the free Hamiltonian, and the even terms are vanishing. This leads to an effective Hamiltonian for pure dephasing of the form \[\begin{align} \label{ex:K95s95diag} K(t) = \frac{\omega_r(t)}{2}\sigma_z \;, \end{align}\tag{8}\] where only the odd ordered cumulants of the bath contribute to the time-dependent energy shift.

This result agrees with time-dependent shifts obtained in [39] (dispersive spin oscillator) and in [11] (pure qubit-qubit dephasing). It furthermore explains the absence of renormalization found in [39] (displaced spin oscillator), and in the well-known dephasing spin-boson model [40]: \[\begin{align} H = \frac{\omega}{2}\sigma_z + \sum_j \omega_j a^\dagger_j a_j + \sum_j g_j \sigma_z \otimes (a^\dagger_j + a_j) \end{align}\] with an initial state of the environment that is diagonal in the Fock basis (including, in particular, a thermal state), for which the exact TCL master equation (already exact at second order [27]) is well known to be \[\begin{align} \dot{\rho}_S(t) = -i \left[\frac{\omega}{2}\sigma_z, \rho_S(t)\right] + \gamma(t) [\sigma_z \rho_S(t)\sigma_z - \rho_S(t)]\;, \end{align}\] where \(\gamma(t)\) is the decoherence rate.

This drastic difference between the nature of even and odd contributions of the series is not just limited to dephasing models. Let us consider an unbiased spin model, where the interaction term is described by an operator that is orthogonal, under the Hilbert-Schmidt scalar product, to the Hamiltonian term, such as \(A=\sigma_x\).

In this case, the odd terms of the series do not commute with the free Hamiltonian, i.e., \([K_{2m+1},H_S]\neq 0\), inducing a rotation of the eigenvectors of the Hamiltonian (see Appendix 4). In contrast, the even terms do commute with \(H_S\), contributing to a shift in the energy spectrum. Thus, when odd orders vanish, the effective Hamiltonian only results in a renormalization of the spin frequency, as expressed in eq. (8 ). Notably, this is the case for a spin-boson model at finite temperature: \[\begin{align} H = \frac{\omega}{2} \sigma_z + \sum_j \omega_j a^\dagger_j a_j + \sum_j g_j \sigma_x \otimes (a^\dagger_j + a_j) \;. \end{align}\] This result proves analytically what was conjectured in [41] due to numerical evidence.

Commonly studied scenarios are those in which the dynamics is characterized by two-point correlation functions, as in the case of linearly coupled, zero-mean Gaussian baths. The present analysis clarifies why such models exhibit simpler renormalization effects: in the case of pure dephasing, renormalization is entirely absent, while for orthogonal couplings, the renormalized Hamiltonian remains diagonal in the original eigenbasis. In contrast, when the bath distribution exhibits time asymmetry (i.e., when odd-order cumulants are nonzero), more intricate behaviors emerge, such as rotations of the Hamiltonian eigenbasis.

Structure of the Hamiltonian —

The examples above show that the effective Hamiltonian exhibits qualitatively distinct behaviors depending on the interaction type and environmental structure. More generally, the perturbative expansion reveals a fundamental symmetry difference between even and odd orders. This feature is not limited to spin models: it arise from decomposing the ordered cumulants of the bath into real and imaginary parts, \(\mathcal{D}=\mathcal{D}_{\mathcal{R}}+i\mathcal{D}_{\mathcal{I}}\), which are respectively symmetric and antisymmetric under time ordering reversal. A similar decomposition holds for the system contribution, where \({\mathbb{A}}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right)\) splits into \({\mathbb{A}_{\mathcal{R}}}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right)\) and \({\mathbb{A}_{\mathcal{I}}}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right)\), which respectively are symmetric and antisymmetric under time ordering reversal.

Making these decompositions explicit in eq. 6 we find that even-order terms \(K_{2m}\) couple bath and system contributions of opposing time symmetry classes, namely \[\begin{align} \label{eq:Kn-even-fd} K_{2m}\! =\!\!\! \sum_{k=0}^{2m} & (-)^{m+k+1} \!\!\int_0^t\!\!\!\! d \boldsymbol{\tau}_{1}^{k} d \boldsymbol{s}_{1}^{n-k}\!\! \left[ \mathcal{D}_{\mathcal{R}}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right){\mathbb{A}_\mathcal{I}}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right)\right.\nonumber\\ &\left.+ \mathcal{D}_{\mathcal{I}}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right){\mathbb{A}_\mathcal{R}}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right)\right] \; , \end{align}\tag{9}\] ensuring an overall antisymmetric structure under time-ordering reversal. Conversely, the odd terms \(K_{2m+1}\) couple terms of the same symmetry class and are overall symmetric \[\begin{align} \label{eq:Kn-odd-fd} K_{2m\!+\!1} \!=\!\!\!\! \sum_{k=0}^{2m+1}\!\!& (-)^{m+k}\!\!\! \int_0^t\!\!\!\!\! d \boldsymbol{\tau}_{1}^{k} d \boldsymbol{s}_{1}^{n-k}\!\! \left[ \mathcal{D}_{\mathcal{I}}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right)\!{\mathbb{A}_\mathcal{I}}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right)\right.\nonumber\\ &\left.- \mathcal{D}_{\mathcal{R}}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right){\mathbb{A}_\mathcal{R}}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right)\right] \;, \end{align}\tag{10}\] highlighting the fundamentally distinct symmetry properties and structure of even and odd orders in the effective Hamiltonian.

From this perspective, we can clearly interpret the terms up to second order. The first order contribution arises only if the interaction-picture bath operator has a non-zero average with respect to the bath’s initial state – often implying the presence of non-zero initial coherences in the environment colla2024FA? – and corresponds to a driving of the system resulting from a displacement of the bath state relative to the interaction Hamiltonian.

The second order contribution, instead, can be understood via linear response theory [42] as arising from fluctuations and dissipation in the environment. In this framework, dissipation is associated with the imaginary part of the correlation function \(\chi(t, s) := i \langle [{B}(t),{B}(s)]\rangle\), i.e. the response function, which is antisymmetric under time reversal. This component reflects the causal response of the environment and accounts for the irreversible energy exchange associated with it. Diffusion, on the other hand, is associated to the real part of the correlation function \(S(t, s) := \frac{1}{2} \langle \{{B}(t),{B}(s)\}\rangle - \langle {B}(t) \rangle\langle {B}(s) \rangle\), which is symmetric under time reversal. This symmetry represents the environmental fluctuations that drive random perturbations in the system state. These random perturbations are responsible for processes such as decoherence, and the spreading of the system’s distribution over time. In terms of these quantities, the second order contribution to \(K\) can be written as \[\begin{align} K_2 = & \int_0^t d\tau \Big[\frac{1}{2i}S(t,\tau) [A(t),A(\tau)] \\ &- \frac{\chi(t,\tau)}{4}\left(\{A(t),A(\tau)\} - \langle\{A(t),A(\tau)\}\rangle_{1/d}\right) \Big] \; . \nonumber \end{align}\]

These insights suggest that distinct environmental features underpin the structural differences observed between even and odd orders of the renormalized Hamiltonian. Although our explicit analysis has been limited to second-order terms, it is reasonable to expect higher-order terms to similarly reflect more intricate environmental correlations.

Conclusions – We introduced a systematic framework to obtain the effective Hamiltonian governing coherent dynamics in non-Markovian open quantum systems. By employing the minimal dissipation principle, we uniquely isolated the coherent component of the time-local generator. We provided an operational expression enabling the direct extraction of the effective Hamiltonian from the generator, as well as a recursive perturbative expansion in terms of bath correlation functions and system-bath interactions.

Our perturbative analysis clarified how environmental correlations and their interplay with system operators lead to energy renormalization, demonstrating clear structural distinctions between even and odd orders of the expansion. These distinctions are used to explain recently observed behaviors in spin systems, including eigenbasis rotations and time-dependent energy shifts. Moreover, our perturbative series can serve as a complementary tool alongside established numerical approaches such as the HEOM, providing analytical guidance and a means to interpret numerical results. This synergy is illustrated in our spin system example, where the perturbative insights help corroborate and elucidate the behavior observed in simulations. Together, these advances have direct relevance for strong-coupling quantum thermodynamics and the engineering of quantum thermal machines.

A.C. acknowledges support from MUR via the PRIN 2022 Project “Quantum Reservoir Computing (QuReCo)” (contract n. 2022FEXLYB). A.C. and H.P.B. acknowledge financial support from the European Union’s Framework Programme for Research and Innovation Horizon 2020 (2014-2020) under the Marie Skłodowska-Curie Grant Agreement No. 847471. G.G. acknowledges financial support by the DAAD, the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation, project numbers 447948357 and 440958198), the SinoGerman Center for Research Promotion (Project M-0294), the German Ministry of Education and Research (Project QuKuK, BMBF Grant No. 16KIS1618K), and the Stiftung Innovation in der Hochschullehre and the Spanish MICIN (project PID2022-141283NB-I00) with the support of FEDER funds.

1 \(K\) in terms of generator fidelity↩︎

The scalar product defined in [19] for two generators reads \[\label{eq:sp95haar} \braket{\mathcal{L}_1,\mathcal{L}_2} = \int\! d \phi \! \int\! d\psi \bra{\phi}\mathcal{L}_1\left[\ket{\psi}\!\!\bra{\psi}\right] \mathcal{L}_2 \left[\ket{\psi}\!\!\bra{\psi}\right] \ket{\phi} \;,\tag{11}\] where \(\psi\) and \(\phi\) are Haar-random generated pure states. The average fidelity [43] between generators of a gate and a noisy channel defined in the main text is then directly proportional to the scalar product 11 of the generators: \[\label{sp95fidelity} \overline{\mathcal{F}}(\mathcal{H}, \mathcal{L})= \int d\psi \mathrm{Tr}\left\{\mathcal{H}\left[\ket{\psi}\bra{\psi}\right] \mathcal{L} \left[\ket{\psi}\bra{\psi}\right] \right\} = d \braket{\mathcal{H},\mathcal{L}}\;,\tag{12}\] since \(\int d \phi \bra{\phi} X \ket{\phi} = \mathrm{Tr}\{X\}/d\). Then, in [19] one chooses a basis of traceless Hamiltonians \(\{H_j\}_j\) such that their associated generators \(\mathcal{H}_j\) form an orthonormal basis for \(\mathfrak{ham}(\mathscr{H})\) under 11 , and find the effective Hamiltonian to be given by \(K = \sum_j H_j \braket{\mathcal{H}_j,\mathcal{L}}\). The Hamiltonian elements \(H_j\) are orthogonal under the Hilbert-Schmidt scalar product but not orthonormal, as \(\mathrm{Tr}\{H_j^2\}= d(d+1)/2\). By defining a normalized basis for Hilbert-Schmidt \[\begin{align} \tilde{H}_j= \sqrt{\frac{2}{d(d+1)}}H_j \end{align}\] we obtain the equation for the effective Hamiltonian in terms of generator fidelity i.e. eq. (4).

2 Derivation of eq. ??↩︎

The derivation is given as in [19] and [20] by projecting the generator onto the subspace \(\mathfrak{ham}(\mathscr{H})\) of superoperators acting as \(-i\) times a commutator with a traceless Hamiltonian, \[\mathcal{H}= \Pi(\mathcal{L}) = \sum_j \mathcal{H}_j \langle \mathcal{H}_j , \mathcal{L}\rangle \; ,\] with \(\{\mathcal{H}_j\}_{j=1}^{d^2-1}\) an orthonormal basis of \(\mathfrak{ham}(\mathscr{H})\), through the scalar product 11 on \(\mathfrak{htp}(\mathscr{H})\). Given the Hamiltonians \(H_j\) associated to the commutators in each \(\mathcal{H}_{j}\), which obey \(\mathrm{Tr}\{ H_i H_j\} = \frac{d(d+1)}{2} \delta_{ij}\), the effective Hamiltonian \(K\) associated to \(\mathcal{H}\) reads \(K = \sum_j H_j \langle \mathcal{H}_j , \mathcal{L}\rangle\) [20]. Using the second moment Haar integral and the cyclic properties of the trace, the coefficients \(\langle \mathcal{H}_j , \mathcal{L}\rangle\) are given simply by \[\langle \mathcal{H}_j , \mathcal{L}\rangle = -\frac{i}{d} \mathrm{Tr}\{H_j\overline{\big[\ket{\phi}\bra{\phi}, \mathcal{L}[ \ket{\phi}\bra{\phi}]\big]}\} \; .\] Given that for any operator \(X\) it holds that \(\sum_j H_j \mathrm{Tr}\{ H_j X\} = \frac{d (d+1)}{2 } X\), we find \[\begin{align} K =& \frac{d+1}{2i} \overline{\big[\ket{\phi}\bra{\phi}, \mathcal{L}[\ket{\phi}\bra{\phi}]\big]} \nonumber \\ =& (d+1) \text{Im}\Big(\overline{\ket{\phi}\bra{\phi}\mathcal{L}[\ket{\phi}\bra{\phi}]}\Big) \; . \end{align}\] To evaluate the fourth moment Haar integral for an unknown generator, consider an orthonormal basis \(\{\ket{i}\}_{i=1}^{d}\) of \(\mathscr{H}\), and take a specific element \(\ket{\bar{k}} =: \ket{\phi_0}\) as reference state such that \(\ket{\phi} = U\ket{\bar{k}} = \sum_i u_{i\bar{k}} \ket{i}\) with \(U\) a Haar random unitary transformation. Then \[\ket{\phi}\bra{\phi}\mathcal{L}[\ket{\phi}\bra{\phi}] = \sum_{ijkl} u_{i\bar{k}}u_{k\bar{k}}u^*_{j\bar{k}}u^*_{l\bar{k}} \ket{i}\!\!\bra{j} \mathcal{L}[\ket{k}\!\!\bra{l}] \; .\] The Haar average therefore acts only on the unitary transformation elements, and the fourth moment is given in terms of the Weingarten functions \(\text{Wg}\), [44]: \[\begin{align} \langle u_{i\bar{k}}u_{k\bar{k}}u^*_{j\bar{k}}u^*_{l\bar{k}}\rangle =& (\delta_{ij}\delta_{kl}\delta_{\bar{k}\bar{k}} +\delta_{il}\delta_{jk}\delta_{\bar{k}\bar{k}})(\text{Wg}[1] + \text{Wg}[2])\nonumber \\ =& \frac{1}{d(d+1)} (\delta_{ij}\delta_{kl}+ \delta_{il}\delta_{jk}) \; , \end{align}\] which gives \[\begin{align} \overline{\ket{\phi}\!\!\bra{\phi}\mathcal{L}[\ket{\phi}\!\!\bra{\phi}]} = \frac{1}{d(d+1)} \sum_{jk}&\Large( \ket{j}\!\!\bra{j} \mathcal{L}[\ket{k}\!\!\bra{k}] \nonumber \\ & + \ket{j}\!\!\bra{k} \mathcal{L}[\ket{k}\!\!\bra{j}] \Large)\; . \end{align}\] The first term is Hermitian and drops out of the imaginary part, so that \[\begin{align} K =& \frac{1}{d} \sum_{j,k=1}^d \text{Im}\Big( \ket{j}\bra{k}\mathcal{L}[\ket{k}\!\!\bra{j}]\Big) \nonumber \\ =& \frac{1}{2i d} \sum_{j,k=1}^d \big[\ket{j}\!\!\bra{k}, \mathcal{L}[\ket{k}\!\!\bra{j}]\big] \; . \label{eq:Ks95jk} \end{align}\tag{13}\] Since the choice of the basis \(\ket{i}\) is arbitrary, one can perform any unitary transformation of the basis \(\{\ket{j}\!\!\bra{k}\}_{j,k}\) onto any basis of operators \(F_{\alpha}\) that is orthonormal with respect to the Hilbert-Schmidt product of operators, see e.g. [45], to find eq. ?? .

Basis choice: generators of \(\mathfrak{su}(d)\)↩︎

Let us choose the orthonormal basis \(F_0=\mathbb{1}/\sqrt{d}\) and \(F_i=\sigma_i/\sqrt{d}\) for \(i\neq 0\) with \(\sigma_i\) generators of \(\mathfrak{su}(d)\), with algebraic relations \[\sigma_j \sigma_k = \delta_{jk} +\sum_{l=1}^{d^2-1} (i f_{jkl}+ d_{jkl})\sigma_l\] where \(f_{jkl}\) are the completely antisymmetric structure constants of \(\mathfrak{su}(d)\), and \(d_{jkl}\) completely symmetric. Then the contribution from \(F_0\) drops out and the effective Hamiltonian is given in terms of the generators only: \[\label{Ks95sigma} K = \frac{1}{2i d} \sum_{j=1}^{d^2-1} \left[\sigma_j, \mathcal{L}[\sigma_j] \right] \; .\tag{14}\] This expression might be particularly useful for systems with low-dimensional Hilbert spaces. If we also write the generator in terms of this basis, \(\mathcal{L}=\sum_{\alpha \beta} l_{\alpha \beta} F_{\alpha} \cdot F_{\beta}\) we get simply \[K = \frac{1}{2i d} \sum_{j} \left(l_{0j} - l_{j0} \right) \sigma_j \; .\]

3 Initial correlations↩︎

Following [31], initial system–environment correlations do not modify the Hamiltonian part of the reduced dynamics. The composite initial state can be written as \[\begin{align} \rho_{SE}(0) = \rho_S \otimes \rho_E + \chi, \end{align}\] where \(\rho_S = \mathrm{Tr}_E\{\rho_{SE}\}\), \(\rho_E = \mathrm{Tr}_S\{\rho_{SE}\}\), and the correlation operator \(\chi = \rho_{SE} - \rho_S \otimes \rho_E\) satisfies \(\mathrm{Tr}_E\{\chi\} = \mathrm{Tr}_S\{\chi\} = 0\). The presence of \(\chi\) introduces an inhomogeneous contribution to the reduced dynamics, leading to an affine time-local generator that can be uniquely extended to a linear map on the space of trace-one operators, without affecting the Hamiltonian term. Explicitly \[\begin{align} \label{eq:corrgenerator} \mathcal{L}^{\chi}_t [ X ] = \mathcal{L}_t[X] + \mathcal{J}^{\chi}_t \mathrm{Tr}\{ X\} \; , \end{align}\tag{15}\] where \(\mathcal{L}_t\) is the generator for initially uncorrelated states and \[\begin{align} \mathcal{J}_t^\chi = \dot{I}_t^\chi - \mathcal{L}_t[I_t^\chi] \; , \quad I_t^{\chi} = \mathrm{Tr}_E \{ U_t \chi U_t^{\dagger} \} \; . \end{align}\] The effective Hamiltonian in the presence of initial correlations can then be identified by applying eq. ?? to the generator in eq. 15 . Importantly, this yields the same Hamiltonian part \(K(t)\) as in the uncorrelated case. Thus, initial correlations affect only the dissipative (non-Hamiltonian) part of the evolution, leaving the effective Hamiltonian invariant.

4 Contributions to \(K\) in spin systems↩︎

4.1 Pure decoherence↩︎

Suppose that \(A=\sigma_z\), so that \(A(t)\equiv A\) \(\forall t\), namely a situation of pure decoherence of the spin. Then we have \({\mathbb{A}}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right) = \langle \sigma_z^{n-k}\rangle_{1/d} \sigma_z^{k}\), which is easily found depending on the parity of \(n\) and \(k\): \[{\mathbb{A}}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right) = \begin{cases} \sigma_z & \text{(n odd, k odd)} \\ \mathbb{1} & \text{(n even, k even)} \\ 0 & \text{elsewhere} \end{cases}\; .\] This means that there are no even order terms in the effective Hamiltonian – as the final Hamiltonian is always traceless and any term proportional to the identity will vanish – and that any odd order term will give a contribution proportional to \(\sigma_z\). Thus, we conclude that for a spin under pure decoherence, the effective Hamiltonian is in general of the form \(K(t) = \frac{\omega(t)}{2}\sigma_z\).

4.2 Unbiased spin system↩︎

Suppose the interaction term to be orthogonal to \(\sigma_z\). For example, take \(A=\sigma_x\). Then \(A(t) = \sigma_+ e^{i\omega t} + \sigma_- e^{-i\omega t}\), and the terms \(A(\boldsymbol{t}_1^m)\) have a different structure depending on the parity of \(m\). For odd \(m\): \[A(\boldsymbol{t}_1^m) = \sigma_+ e^{i\phi(\boldsymbol{t}_1^m)} + \sigma_- e^{-i\phi(\boldsymbol{t}_1^m)}\;,\] while for even \(m\): \[\begin{align} A(\boldsymbol{t}_1^m) &= \ket{1}\bra{1} e^{i\phi(\boldsymbol{t}_1^m)} + \ket{0}\bra{0} e^{-i\phi(\boldsymbol{t}_1^m)}\nonumber\\ &={1{\boldsymbol{l}}}\, \cos(\phi(\boldsymbol{t}_{1}^{m})) +i\sigma_z\, \sin(\phi(\boldsymbol{t}_{1}^{m})) \;, \end{align}\] where the phases are defined as \(\phi(\boldsymbol{t}_1^m):= - \omega \sum_{j=1}^m (-1)^j t_j\). This leads to \[{\mathbb{A}}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right) = \begin{cases} {\mathbb{A}_1}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right)& \text{(n odd, k odd)} \\ {\mathbb{A}_2}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right) & \text{(n even, k even)} \\ 0 & \text{elsewhere} \end{cases}\; ,\] with \[\begin{align} {\mathbb{A}_1}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right) =& \cos(\phi(\boldsymbol{s}_1^{n-k}))\left[\sigma_+ e^{i\phi(\boldsymbol{\tau}_1^k)} + \sigma_- e^{-i\phi(\boldsymbol{\tau}_1^k)}\right] \\ {\mathbb{A}_2}\!\left(\boldsymbol{\tau}_{1}^{k}\!, \boldsymbol{s}_{1}^{n-k}\right) =& \cos(\phi(\boldsymbol{s}_1^{n-k}))\left[\mathbb{1} + i \sin(\phi(\boldsymbol{\tau}_1^{k}))\sigma_z \right] \end{align}\]

From the above it is clear that even orders give contributions on the diagonal of \(K\), thus maintaining a structure proportional to \(\sigma_z\), while odd orders give only contributions on off-diagonals of the effective Hamiltonian, thus changing its eigenbasis.

References↩︎