## Driving of an open quantum system at finite temperature across first- and second-order quantum phase transitions

Felipe Matus ∗ and Pavel Cejnar †

Institute of Particle and Nuclear Physics, Faculty of Mathematics and Physics, Charles University, V Holeˇ soviˇk´ch c a 2, 180 00 Prague, Czechia

(Dated: December 18, 2024)

An open fully connected system of qubits at nonzero temperature is driven within a finite time interval along various paths in the space of its control parameters. The driving leads across finitesize precursors of firstand second-order quantum phase transition from factorized to entangled ground-state phases, aiming at the preparation of the complex ground state of the system at the final parameter point with maximal fidelity. During the drive, the system is coupled to a heat bath at a constant temperature, the dynamics being determined in a nonperturbative way by the method of Hierarchical Equations of Motion. It is shown that the presence of the heat bath in combination with specific patterns of avoided crossings affecting the ground and excited states in the parameter region around the quantum phase transition may considerably improve the fidelity of preparation of the target ground state.

## 1. INTRODUCTION

One of the aims of arising quantum information science is to develop efficient and reliable techniques for preparing general correlated states of complex quantum systems [1, 2]. Acommon strategy is based on driving a suitable system involving complex interactions among its elementary constituents along a selected path in the space of control parameters [3-6]. The path starts at an initial point corresponding to an uncorrelated, easy-to-prepare ground state, and terminates at a final point, where the ground state contains the desired complex quantum correlations. Reproducing the target ground state with high fidelity and in moderate time is a task that has direct applications in quantum control, quantum computation and information processing [7].

driving increases the fidelity of the final-state preparation by repeated measurement-like interactions of the driven system with a suitable ancilla [21-24].

In the last decades different strategies have been proposed to accomplish this task. One can mention, among the most striking approaches, the adiabatic, transitionless, geometric, and decoherence-assisted types of driving. The adiabatic driving represents a direct application of the quantum adiabatic theorem, which states that for a sufficiently slow variation of control parameters, the system follows the instantaneous (adiabatic) ground state in the parameter space [8-10]. The transitionless (also called counter-diabatic) driving emulates the adiabatic evolution along a given parameter path with the aid of additional terms in the Hamiltonian (which usually are non-local) that compensate non-adiabatic effects of the finite-speed driving [11-14]. The geometric driving makes use of the geodesic path between the initial and final parameter points with respect to the Provost-Vallee metric on the ground-state manifold in attempt to maximize the overlap of the evolving state with the adiabatic ground state [15-20]. Finally, the decoherence-assisted

∗ matus@ipnp.mff.cuni.cz

† cejnar@ipnp.mff.cuni.cz

Some of these driving strategies were recently tested by us in systems composed of a single or multiple interacting qubits [24-26]. The multi-qubit environment was implemented through a model from the Lipkin-MeshkovGlick (LMG) family [27], which is numerically treatable and available to experiments, but simultaneously shows complex phenomena like ground-state and excited-state quantum phase transitions [28-32]. The ground-state quantum phase transitions (QPTs), as sudden changes (non-analytic in the limit of infinite system size) of the ground-state energy and wave function with Hamiltonian control parameters [33], are of a straightforward importance for the driving problem if the initial and target states belong to different ground-state phases. Crossing of the QPT point usually sets the most stringent bounds on the total driving times and/or the final fidelity achieved [34-36]. On the other hand, excited-state quantum phase transitions (ESQPTs), which represent an extension of the QPT to the excited domain [37-44], become relevant as soon as quantal or thermal fluctuations create nonvanishing populations of states in the ESQPT domain during the driving process. This typically applies to driving through a wider QPT region, where ESQPT structures appear at low energies and have considerable effect on dynamics of the ground-state occupation.

Our previous studies [24-26] were focused on fully coherent and decoherence-assisted types of driving in strictly isolated systems only. This means that the system was initiated at exactly zero temperature and excited states of the intermediate Hamiltonians were populated solely by quantum transitions induced by nonadiabaticity of the driving. However, in realistic situations, the driven system is likely to interact-at least to a limited extent-with its surrounding environment. This leads to thermal population of excited states during the course of the driving process. If the driving leads across the QPT and ESQPT parameter regions, thermal

noise can excite the system due to reduced energy gaps and relax it afterward. Thus, the impact of temperature on the final fidelity may become relatively strong [45, 46] To analyze such effects is the main purpose of the present paper.

Open systems [47-49] are most commonly studied within the Lindblad formalism, i.e., assuming that the system is weakly coupled to a Markovian reservoir and considering the Born and rotating wave approximations [50-52]. Although the equations obtained in this way are numerically simple to solve, their validity is subject to too stringent conditions and the results are not guaranteed to be accurate [53-59]. In this work, we use the method of Hierarchical Equations of Motion (HEOM) [60-64], which is based on a numerically exact solution of quantum dynamics through the influence functional formalism [65-67]. This approach is applicable to a large span of driving times and allows for narrow energy gaps. The method has already been used to study quantum dissipative dynamics under timedependent driving fields [68], entanglement dynamics [69] as well as more complex dynamics [70-73].

The structure of this paper is as follows: In Sec. 2 we outline the model, along with its quantum critical properties, and describe the quantities that characterize the environment and its interaction with the system. In Sec. 3 we introduce the driving procedures used to evolve the system and briefly describe the HEOM method and its parameters. In Sec. 4 we present and discuss the main results of our numerical calculations. In Sec. 5, we give a brief summary and conclusion.

## 2. MODEL

We consider an interacting system of qubits (hereafter denoted by S) coupled to an external environment (denoted by E). The total Hamiltonian is decomposed into three parts:

$$\hat { H } = \hat { H } _ { \text{S} } ( \boldsymbol \Lambda ) + \hat { H } _ { \text{E} } + \hat { H } _ { \text{I} }, \quad \ \ ( 1 ) \quad \stackrel { \lambda } { \text{int} }.$$

where ˆ H S ( Λ ) is the Hamiltonian of qubits, with Λ denoting the set of control parameters that will be later subject to the driving procedure, ˆ H E is the Hamiltonian of the environment (also called the bath), and ˆ H I stands for the system-environment interaction. The forms of these Hamiltonians are given below.

Note that in this paper, the overall energy scale is given by an implicit constant ε , whose choice is arbitrary. So all components of the Hamiltonian (1) and in general all energies E are considered dimensionless, expressed in units of ε . Since we additionally set ℏ = k B = 1, the temperature T is given in units of ε and the time t in units of 1 /ε .

## 2.1. Qubit system

The system S consists of N &gt; 1 mutually interacting qubits, equivalent to spin1 2 particles, whose dynamics is described within the LMG framework [27]. The Hamiltonian is written in terms of collective spin (quasispin) operators

$$\hat { J } _ { \alpha } = \frac { 1 } { 2 } \sum _ { i = 1 } ^ { N } \hat { \sigma } _ { \alpha } ^ { ( i ) }, \quad \alpha = x, y, z, \quad \quad ( 2 )$$

where ˆ σ ( ) i α are Pauli matrices acting on the i th-qubit space. The operators ˆ J α satisfy the usual SU(2) commutation relations.

We use the LMG Hamiltonian used in our previous studies [24, 25],

$$\hat { H } _ { \text{S} } ( \boldsymbol \Lambda ) = \hat { J } _ { z } - \frac { 1 } { N } \left [ \lambda \hat { J } _ { x } ^ { 2 } & \ + \ \chi \left \{ \hat { J } _ { x }, \hat { J } _ { z } + \frac { N } { 2 } \right \} \\ & \ + \ \chi ^ { 2 } ( \hat { J } _ { z } + \frac { N } { 2 } ) ^ { 2 } \right ], \quad ( 3 )$$

where { A,B } = AB + BA . The first term ˆ J z of this Hamiltonian describes a system of noninteracting qubits, while the terms in square brackets represent interactions connecting any qubit with all the other qubits of the system. The dimensionless interaction strengths λ and χ , jointly denoted as Λ ≡ { Λ µ } µ =1 2 , , are considered as the system control parameters. These will be subject (see Sec. 3.2 below) to an externally driven time variation, Λ ( ) t = ( λ t , χ t ( ) ( )), which will make the system Hamiltonian time-dependent: ˆ H S ( ) t = ˆ ( H S Λ ( )). t As seen in Eq. (3), the above-introduced energy constant ε , in units of which the Hamiltonian is expressed, coincides with the energy difference between the up and down qubit states in absence of interaction ( λ = χ = 0).

The total quasispin ˆ J 2 = ˆ J 2 x + ˆ J 2 y + ˆ J 2 z is conserved by ˆ H S ( Λ ), so the full 2 N -Hilbert space of qubits splits into a sum of (2 j + 1)-dimensional subspaces (almost all of them appearing in numerous replicas) with different values of the quantum number j [74]. We restrict here to the unique ( N + 1)-dimensional subspace with j = N 2 , which is fully symmetric under the exchange of qubits, and assume that this subspace is invariant under the full Hamiltonian (1) since also the interaction term ˆ H I conserves the exchange symmetry. The energy spectrum of ˆ H S ( Λ ) is expressed in an ordered form: E 0 ( Λ ) ≤ E 1 ( Λ ) ≤ · · · ≤ E N ( Λ ).

For χ = 0, another quantity, which is conserved by the system Hamiltonian, is the parity

$$\hat { \mathcal { P } } = ( - 1 ) ^ { \hat { J } _ { z } - \frac { N } { 2 } }.$$

However, in the following we will assume that this symmetry is not necessarily preserved by the interaction Hamiltonian ˆ . H I

The ground state of Hamiltonian (3) exhibits QPTs of first and second order. These can be described in terms of two order parameters associated with the ground-state expectation values ⟨ ˆ J z + N 2 ⟩ and ⟨ ˆ J x ⟩ . For N →∞ , the plane of control parameters λ × χ splits in three groundstate phases: phase I with ⟨ ˆ J z + N 2 ⟩ = ⟨ ˆ J x ⟩ = 0 in the domain λ &lt; λ c ( χ ), where

$$\lambda _ { c } ( \chi ) = 1 - \frac { \chi ^ { 2 } } { 1 - \chi ^ { 2 } }, \quad \quad ( 5 ) \quad \stackrel { \Re } { \ } \\. \hat { \ } \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \vec { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat {. } \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }.. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }... \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \ \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \hat { \ }. \end{cases}$$

phase II with ⟨ ˆ J z + N 2 ⟩ &gt; 0 and ⟨ ˆ J x ⟩ &gt; 0 in the domain λ &gt; λ c ( χ ), χ &gt; 0, and phase III with ⟨ ˆ J z + N 2 ⟩ &gt; 0 and ⟨ ˆ J x ⟩ &lt; 0 in the domain λ &gt; λ c ( χ ), χ &lt; 0. The N →∞ form of the ground state in phase I is fully factorized, corresponding to separated qubits in the down spin projection states, whereas in phases II and III all qubits in the ground-state wave function are mutually entangled. Since phases II and III are symmetric with respect to the χ ↔-χ inversion, we consider here only the case χ ≥ 0 with phases I and II.

̸

The QPTs between these phases are of the first order for χ = 0 and of the second order for χ = 0 [38]. In the first-order QPT, the order parameters change discontinuously, while the ground-state energy E 0 ( Λ ) shows a discontinuity of the first derivative and the ground-state energy gap ∆ 10 ( Λ ) = E 1 ( Λ ) -E 0 ( Λ ) vanishes exponentially with N →∞ . In contrast, the second-order QPT comes through a discontinuity of the first derivative of the order parameter, which is connected with a discontinuity of the second derivative of the ground-state energy and a polynomially vanishing ground-state energy gap. In the second-order QPT, the parity of the ground state is spontaneously broken, which means that for λ &gt; λ c (0) and χ = 0 the N →∞ ground state becomes a degenerate parity doublet. A more elaborate analysis of the QPTs and the ground-state geometry of Hamiltonian (3) can be found in [25, 75]. FiniteN precursors of the QPTs of the first and second order are observed in the evolution of the lowest energy levels in panels (a) and (b) of Fig. 1. This figure depicts lower parts of energy spectra of our LMG Hamiltonian for two particular paths in the λ × χ plane, which will be used (although with a different value of N ) in the driving protocols to be introduced in Sec. 3.

The ESQPT nonanalyticities, affecting in the N →∞ limit higher energy levels, are also present in the spectrum of Hamiltonian (3). Their finiteN precursors are also shown in Fig. 1 (see the parts of the spectra with the colored background). Associated with the first-order QPT, see Fig. 1(a), two chains of avoided level crossings of excited states issue from the λ ≈ λ c ( χ ) avoided crossing of the ground state to both λ &lt; λ c and λ &gt; λ c directions. In addition, many other avoided crossings appear in a finite region above the ESQPTs. The second-order QPT in Fig.1(b) is accompanied by a single chain of avoided crossings in the direction λ ≳ λ c (0). These structures, which are generic accompaniments of the first- and second-order QPTs [40, 44], play an important role in

Figure 1. Energy spectra of the LMG Hamiltonian (3) for two cuts through the λ × χ parameter plane: (a) λ ∈ [0 , 0 5], . χ = 4 8 . λ and (b) λ ∈ [0 , 4], χ = 0. We set N = 40. The vertical dotted line segments demarcate N →∞ critical points λ c ( χ ) of (a) the first-order QPT and (b) the second-order QPT. In panel (b) the positive and negative parity levels are drawn by full and dashed curves, respectively, visualizing the spontaneous parity breaking at the second-order QPT. All level crossings are avoided. Their patterns (highlighted by the background color) around both ground-state QPTs represent precursors of ESQPTs.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00635v4-figure-a9840735-214e-407a-9bfd-08f2b132e19e.jpg)

driving-induced dynamics of the system at nonzero temperature.

## 2.2. Environment and qubit-environment coupling

We consider the prototypical environment consisting of a collection of a large number of harmonic oscillators in thermal equilibrium at temperature T = 1 /β . The Hamiltonian (in units of ε ) of this heat bath is

$$\hat { H } _ { \text{E} } = \sum _ { k } \omega _ { k } \hat { b } _ { k } ^ { \dagger } \hat { b } _ { k },$$

where ˆ b † k and ˆ b k denote the bosonic creation and annihilation operators and ω k are frequencies of individual modes (phonon energies). To describe the qubitenvironment interaction, we follow the Caldeira-Leggett approach [66, 67], where the qubit system is assumed to be coupled linearly to the position of all the environmental oscillators. Therefore, the interaction Hamiltonian (again in units of ε ) is expressed as

$$\hat { H } _ { 1 } = \frac { \dot { \, 1 } } { \sqrt { N } } \, \hat { Q } \otimes \underbrace { \sum _ { k } g _ { k } ( \hat { b } _ { k } ^ { \dagger } + \hat { b } _ { k } ) } _ { \hat { x } }, \, \quad \quad ( 7 )$$

where ˆ is a weighted position of the oscillators, with X g k denoting non-negative coupling strengths to individual

modes, and ˆ is the system-environment coupling operQ ator acting in the space of qubits. The operator ˆ Q is naturally taken as a linear combination of the quasispin components ˆ , for instance J α

$$\hat { Q } = \sin \theta \, \hat { J } _ { x } + \cos \theta \, \hat { J } _ { z }. \quad \quad ( 8 ) \quad \text{ w}$$

We use two choices of the angle θ , namely θ = 0 and θ = π 2 . While in the former case the environment just modulates energies of the unperturbed qubit states, in the latter case it induces transitions between these states. Note that with θ = 0 and χ = 0 the total Hamiltonian (1) conserves the parity (4), so for χ = 0 the choice θ = π 2 captures the generic situation in which the qubitenvironment interaction violates non-fundamental symmetries of the qubit system. The factor in front of the expression (7) ensures proper scaling of the S-E interaction energy with an increasing number of qubits N .

The total Hamiltonian (1) with individual terms fixed by Eqs. (3), (6) and (7) suffers from the drawback that for infinite sizes of both the qubit and environment subsystems and for sufficiently large coupling strengths g k its ground state energy may cease having a lower bound [76]. In case of a finite N , the lower bound exists, but the whole system with increasing g k 's may exhibit a sharp crossover from a normal phase, characterized by low average occupancies ⟨ ˆ ˆ b † k b k ⟩ of all bath modes, to a 'superradiant' phase, in which these occupancies get macroscopically large [77-79]. This transition is experimentally verified in various laboratory realizations of the Dicketype Hamiltonians in which the bilinear coupling of the form (7) is also present, see, e.g., [80, 81].

To avoid the ground-state instability and the transition to the superradiant-like phase, the system Hamiltonian (3) is often renormalized by adding an extra term, which counterbalances the influence of the environment [66, 67]. This leads to a tranformation

$$\hat { H } _ { S } ( \boldsymbol \Lambda ) \longmapsto \hat { H } _ { S } ^ { \prime } ( \boldsymbol \Lambda, q ) = \hat { H } _ { S } ( \boldsymbol \Lambda ) + r \frac { q } { N } \hat { Q } ^ { 2 }, \quad ( 9 ) \quad \text{$\underset { \boldsymbol \tau } { \text{ratio} } }$$

↦

where variable r will be explained below and

$$q = \sum _ { k } \frac { g _ { k } ^ { 2 } } { \omega _ { k } } \quad \quad \ ( 1 0 ) \quad \ ( 1 }$$

plays the role of an effective S-E coupling strength. With the replacement (9), the total Hamiltonian (1) is guaranteed to have a lower bound of energy [82].

In this work, the qubit system is always considered finite, so the inclusion of the counterterm is not necessary. Nevertheless, as we wish to estimate the effect of renormalization in the driving process, we include the counterterm to some of our simulations. Whether the renormalization is on or off, respectively, is distinguished by values r = 1 or r = 0 of the auxiliary parameter in Eq. (9). In particular, the option r = 1 is used along with ˆ = Q ˆ , J x when the coupling term (7) can be interpreted as the electric dipole interaction of the qubits (atoms carrying the qubit states) with ω k modes of a quantized electromagnetic field and the counterterm stands for the A 2 term (squared vector potential) of the atomfield Hamiltonian [76, 81]. Hence for θ = π 2 in Eq. (8) we use both r = 0 and r = 1 options, while for θ = 0 we use only r = 0. With the choice r = 0, we intentionally allow for the situation in which the increasing effective strength (10) of the S-E interaction induces a crossover effect in the ground state of the whole system to a superradiant-like phase, which may have a large impact on results of the driving process.

Thermal fluctuations of the environment are described by its self-correlation function

$$\stackrel { \text{by its set-correlation function} } { \text{bit-} } \\ \text{ym-} \quad C ( t ) & = \text{tr} \left [ \hat { X } ( t ) \hat { X } ( 0 ) \, \hat { \rho } _ { E } ^ { ( \beta ) } \right ] \quad ( 1 1 ) \\ \text{rac-} \quad & = \sum _ { k } g _ { k } ^ { 2 } \left [ \left \{ 2 n _ { E } ( \omega _ { k } ) + 1 \right \} \cos ( \omega _ { k } t ) - i \sin ( \omega _ { k } t ) \right ] \\ \text{ixed} \quad & = \frac { 1 } { \pi } \int _ { 0 } ^ { \infty } d \omega \, J ( \omega ) \left [ \coth \left ( \frac { \beta \omega } { 2 } \right ) \cos ( \omega t ) - i \sin ( \omega t ) \right ], \\ \text{k} \text{ its window } \, \stackrel { \text{when} } { \wedge } \, \hat { ( \beta ) } \text{ is a human-identity, even for the union}$$

where ˆ ρ ( β ) E is a thermal density operator of the environment and ˆ ( X t ) = e + ˆ iH E t ˆ Xe -iH ˆ E t represents time evolution of the weighted oscillator coordinate from Eq. (7). In the second line we use the Bose distribution n E ( ω k ) = ( e βω k -1) -1 and in the third line we introduce the spectral density

$$J ( \omega ) = \pi \sum _ { k } g _ { k } ^ { 2 } \delta ( \omega - \omega _ { k } ), \text{ \quad \ \ } ( 1 2 )$$

which fully determines the correlation function [47, 48]. For discrete modes, the spectral density consists of a sequence of δ -peaks, but in the limit of a continuous distribution of bath mode energies (infinite number of bosonic modes) it becomes a smooth function of ω . While the real part of the correlation function (11) describes thermal fluctuations and decoherence, the imaginary part captures the effects of dissipation. For more detailed explanation we suggest Refs. [47-49, 60].

It is assumed that the dissipative effect of the bath decays as Im C t ( ) ∼ e -γt , where γ -1 = τ E is the dissipation time scale. By imposing this condition one finds (12) in the Drude-Lorentz form [60, 64, 69]:

$$J ( \omega ) = 2 q \, \frac { \gamma \omega } { \gamma ^ { 2 } + \omega ^ { 2 } }, \text{ \quad \ \ } ( 1 3 )$$

where γ represents the width of the bath mode distribution. With this expression, it is possible to write the discreet sum in (10) as the integral q = π -1 ∫ ∞ 0 dωJ ω /ω ( ) .

The use of this formalism in numerical calculations based on the HEOM method will be further commented in Sec. 3.4.

## 3. DRIVEN DYNAMICS

In the driving procedure, the total Hamiltonian (1) is made time dependent by imposing an external time

dependence Λ ( ) of control parameters of the system of t qubits. Thus, ˆ ( H t ) = ˆ H S ( Λ ( ) t ) + ˆ H E + ˆ . H I Below we describe the evolution induced by this Hamiltonian and identify the target states of the qubit system.

## 3.1. Initial state and overall evolution

We assume that the system S and the environment E are at the initial time t = 0 prepared with the same inverse temperature β in a factorized initial state:

$$\hat { \rho } ( 0 ) = \underbrace { \frac { e ^ { - \beta \hat { H } _ { S } ( 0 ) } } { \hat { \rho } _ { S } ( 0 ) } } _ { \hat { \rho } _ { E } ( 0 ) } \otimes \underbrace { \frac { e ^ { - \beta \hat { H } _ { E } } } { \hat { \rho } _ { E } ( 0 ) = \hat { \rho } _ { E } ^ { ( \beta ) } } } _ { \hat { \rho } _ { E } ( 0 ) = \hat { \rho } _ { E } ^ { ( \beta ) } }. \quad \ \ ( 1 4 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
 \hat { \hat { \hat { t } } } } _ { \hat { \hat { t } } }$$

Here, ˆ H S (0) is the system Hamiltonian at t = 0, and ˆ (0) and ˆ ρ S ρ E (0) are thermal density operators of S and E, respectively. This corresponds to the situation when in times t &lt; 0 both S and E have thermalized separately through interactions with a common thermal reservoir but with no interaction between each other. We note that the thermal initial state of S in Eq. (14) is probably a more realistic choice than a pure zero-temperature state ˆ (0) = ρ S | E 0 (0) ⟩⟨ E 0 (0) . | Moreover, it will be shown that the nonzero temperature of S has a non-trivial effect on the final fidelity of the driving procedure.

The interaction between S and E is initiated at t = 0, indicating the start of the driving procedure on S, see Sec. 3.2. One can assume that the environment E is connected with the instrumentation needed to perform the drive. Since then, the S-E interaction completely dominates over the interaction with the other reservoir so that the latter can be neglected. The evolved state ˆ( ) ρ t of the composite S-E system for t &gt; 0 satisfies the quantum Liouville-von Neumann equation

$$\frac { d } { d t } \hat { \rho } ( t ) = - i \left [ \hat { H } ( t ), \hat { \rho } ( t ) \right ]. \quad \quad ( 1 5 ) \quad \text{$$

It is very probable that because of the interaction between S and E, the total density operator ˆ( ) is no more ρ t factorized. The density operator of the qubit system is then extracted by the common procedure involving the partial trace over the Hilbert space of the environment:

$$\hat { \rho } _ { \text{S} } ( t ) = \text{tr} _ { \text{E} } \hat { \rho } ( t ). \text{ \quad \ \ } ( 1 6 )$$

We outline the method how this evolution is determined in Sec. 3.4. Before, we focus on the externally driven variation of the system Hamiltonian ˆ H S ( Λ ( ) t ) and on the definition of the state of the qubit system which is the target of our driving procedure.

## 3.2. Driving protocols

The time-dependent system Hamiltonian ˆ H S ( Λ ( ) t ) follows from a predefined time dependence of control parameters Λ of Hamiltonian (3). The function Λ ( ) t is

Figure 2. Top panel: Relevant part of the parameter space of the qubit Hamiltonian ˆ H S ( Λ ). The dashed curve represents the N →∞ QPT separatrix between the ground-state phases I and II. The tilted (solid red) line segment represents the driving path across the first-order QPT. The horizontal (solid blue) line segment represents the driving path across the second-order QPT. The gray background encodes the energy gap ∆ 10 ( Λ ) for N = 10. Lower panels: The planar speed u (solid line) and geometric speed v (dotted line) for driving protocols A and B along both paths for N = 10 and t F = 1. The upper and lower rows (with red and blue curves) correspond to the paths across the first- and second-order QPTs, respectively (overflow parts of the dependencies are displayed in the insets). The vertical dashed line segments mark the crossing of the N →∞ QPT separatrix from the upper panel.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00635v4-figure-35b6e3a3-eb48-449b-9bc7-5cd69c69ca9a.jpg)

determined by the chosen path in the parameter space and by the speed with which the system is driven along this path. In our case, the path will always be a straight line connecting the initial and final parameter points Λ I and Λ F , respectively. This line is parametrized as

$$\Lambda ( t ) = \Lambda _ { I } + s ( t ) \left ( \Lambda _ { F } - \Lambda _ { I } \right ), \text{ \quad \ \ } \left ( 1 7 \right )$$

where s t ( ) is a variable satisfying s (0) = 0 and s t ( F ) = 1, with t F denoting the final time, i.e., the total time of the driving. This ensures that Λ (0) = Λ I and Λ ( t F ) = Λ F . Moreover, we assume a continuous and monotonous variation of s t ( ) between its initial and final values 0 and 1, and consider the following two options for its time dependence.

Protocol A: The first option is the simplest linear dependence s t ( ) = t/t F . This ensures that the speed on the

parameter plane

$$u ( t ) = \sqrt { \sum _ { \mu, \nu } \delta _ { \mu \nu } \frac { d \Lambda ^ { \mu } } { d t } \frac { d \Lambda ^ { \nu } } { d t } }, \quad \quad ( 1 8 ) \quad \begin{smallmatrix} ( \lambda, \\ \text{dri} \\ \text{qu} \\ \text{in} \end{smallmatrix} }$$

hereafter referred to as the 'planar speed', is kept constant, equal to u = | Λ F -Λ I | /t F .

Protocol B: The second option is a more complicated time dependence of parameter s t ( ), that keeps constant the so-called 'geometric speed'

$$v ( t ) = \sqrt { \sum _ { \mu, \nu } g _ { \mu \nu } ( \Lambda ( t ) ) \frac { d \Lambda ^ { \mu } } { d t } \frac { d \Lambda ^ { \nu } } { d t } }. \quad \quad ( 1 9 ) \quad$$

Here, g µν ( Λ ) is the metric tensor defined on the manifold of ground states | E 0 ( Λ ) ⟩ of Hamiltonian ˆ H S ( Λ ) via the distance element dl 2 = 1 -|⟨ E 0 ( Λ ) | E 0 ( Λ + d Λ ) ⟩| 2 [15]. Assuming a non-degenerate spectrum for finite N , we can write

$$g _ { \mu \nu } = \text{Re} \sum _ { n > 0 } \frac { \langle E _ { 0 } | \frac { \partial } { \partial \Lambda ^ { \mu } } \hat { H } _ { \text{S} } | E _ { n } \rangle \langle E _ { n } | \frac { \partial } { \partial \Lambda ^ { \nu } } \hat { H } _ { \text{S} } | E _ { 0 } \rangle } { ( E _ { n } - E _ { 0 } ) ^ { 2 } }, \ \ ( 2 0 ) \quad \text{$\text{$\text{$\text{$\text{$th}$} }$}$} \\. \quad. \. \. \. \. \. \. \. \. \. \. \.$$

where for brevity we skip all dependencies on Λ .

The metric (20) and the geometric speed (19) reflect variations of the ground-state wave function with control parameters Λ . In protocol B, the geometric speed is kept constant, equal to v = ∫ Λ F Λ I √ g µν d Λ µ d Λ µ /t F , while the planar speed (18) varies in time, taking smaller values in the parameter domains where g µν ( Λ ) is large, i.e., particularly in regions with a small energy gap ∆ 10 ( Λ ). For drivings of isolated systems initiated in the ground state | E 0 ( Λ I ) ⟩ , the v =const. strategy improves the final fidelity for the target state | E 0 ( Λ F ) ⟩ [25, 26] (see also [83] for a closely related method). We will show that this advantage lasts even at some nonzero (sufficiently low) temperatures. In principle, one can go beyond protocol B by using a curved path that minimizes the geometric length from Λ I to Λ F , but this approach was shown to yield unsatisfactory results [25] and we do not apply it here.

The initial and final points in the Λ = ( λ, χ ) plane define two driving paths:

$$\Lambda _ { \mathbf I } & = ( 0, 0 ), & ( 2$$

$$\Lambda _ { F } = \begin{cases} ( 0. 2 5, 1.$$

This means that the driving always starts at Λ I corresponding to the factorized ground state of the system of independent qubits. In contrast, both final points Λ F in Eq. (22) are associated with an entangled ground state of the system of strongly interacting qubits. For the first choice of the final point, the driving trajectory crosses the first-order QPT, the parity (4) being broken. For the second choice the trajectory crosses the second-order QPT along the parity-conserving line with χ = 0. The crossing of the respective critical point (its N →∞ realization) corresponds to s ≈ 0 57 for the first-order QPT . path and s = 0 5 for the second-order QPT path. .

In Fig. 2 we plot the two driving paths in the plane ( λ, χ ) and the speeds (18) and (19) along both paths for driving protocols A and B for the system with N = 10 qubits. We note the reduction of the planar speed u t ( ) in both firstand second-order QPT regions for the v =const. driving procedure B. This slowdown is correlated with a reduced ground-state energy gap ∆ 10 , the minimal gap for N = 10 being however shifted with respect to the N →∞ QPT separatrix.

## 3.3. Target states and fidelity

The target state for any of the above-described driving procedures is the ground state of the qubit system at the final parameter point Λ = Λ F . For the path across the first-order QPT with Λ F = (0 25 1 2), . , . the final ground state is given by a single eigenvector | E 0 ( Λ F ) ⟩ , which is energetically well separated from the other eigenvectors. So, if ˆ ( ρ S t F ) is the density operator of the qubit system at the end of the driving procedure, see Eq. (16), then the fidelity of the target-state preparation is

$$\mathcal { F } _ { 1 } = \langle E _ { 0 } ( \bold$$

This formula expresses the probability of finding the target ground state | E 0 ( Λ F ) ⟩ in the statistical ensemble associated with the final density operator ˆ ( ρ S t F ).

For the χ = 0 path across the second-order QPT with Λ F = (2 0), the final ground state in the , N →∞ limit is a degenerate doublet of positive- and negative-parity states. Although the degeneracy of the positive-parity eigenstate | E 0 ( Λ F ) ⟩ and the negative-parity eigenstate | E 1 ( Λ F ) ⟩ is not exact for finite N , we calculate the final fidelity by summing the overlaps with both these states:

$$\begin{array} {$$

where ˆ P gs ( Λ F ) is the projector to the subspace spanned by vectors | E 0 ( Λ F ) ⟩ and | E 1 ( Λ F ) ⟩ . This formula expresses the probability that any state randomly drawn from the statistical ensemble associated with ˆ ( ρ S t F ) lies in the quasidegenerate ground-state subspace at Λ = Λ F .

## 3.4. Notes on the HEOM calculations

The evolution of the qubit density operator ˆ ( ) from ρ S t Eq. (16) is determined by the HEOM method [60-64], which is capable to yield exact (nonperturbative) numerical description of the system dynamics. In particular, it goes beyond the approximations involved in the common Lindblad formalism, which strictly relies on the assumption that environment-induced processes on the system are Markovian. The method originates from the influence functional formalism of Feynman and Vernon [65], which arises from the integration over all degrees of freedom

of the environment and determines the system evolution ˆ ( ) on the basis of the system self-Hamiltonian ρ S t ˆ H S , the coupling operator ˆ Q from Eq. (7) and the environment self-correlation function C t ( ) from Eq. (11).

We do not immerse here deeper into the description of the HEOM method but just present its very rough outline and introduce the parameters involved. The evolution ˆ ( ) ρ S t is described by a set of coupled first-order differential equations [60], which can be efficiently solved in terms of some auxiliary density matrices. These matrices do not represent real physical states but serve just as a numerical tool to account for the S-E correlations depending on the coupling strength q , temperature T and dissipation time scale τ E . When q is small and τ E much lower than the internal dynamical time scale of the system, the dynamics is approximately Markovian and only a few auxiliary density matrices are needed. In contrast, for large q and τ E the dynamics is non-Markovian and a large number of matrices is required to account for the long-lasting S-E correlations. The auxiliary matrices are grouped under a certain hierarchy of layers enumerated by 0 1 , , . . . , L , which are connected with the inclusion of time correlations of an increasing order. The depth L of this hierarchy must be chosen large enough to ensure the convergence. Details can be found, e.g., in Ref. [61].

In order to solve the set of HEOM equations, the selfcorrelation function (11) is expanded as a series of exponentials. This can be done analytically, in general, for a spectral density J ω ( ) that admits a finite number of poles in the complex plane and that goes to zero as ω → ∞ . For the Drude-Lorentz distribution (13) the series takes the form

$$C ( t ) = c _ { 0 } \, e ^ { - \gamma t } +$$

with coefficients c 0 ∈ C and c 1 , c 2 , . . . ∈ R depending on parameters q and γ , and on the temperature T . The explicit expressions are given, e.g., in Ref. [64]. In numerical calculations the sum in Eq. (25) can be truncated at its M th term. For low temperatures, the convergence is slow and M must be relatively large to reach an acceptable precision of calculations. On the other hand, for high temperatures even very small M may lead to reasonable results. In any case, to compensate this truncation, a Lindbladian term that speeds up the convergence (socalled terminator) is included [53].

During the last decades, several open-source libraries have been developed to implement general calculations based on the HEOM formalism. Here we use the version found in the Python library QuTiP-BoFiN, see Ref. [64].

## 4. RESULTS

In our calculations we use the value γ = 10 for the width of the Drude-Lorentz distribution (13) in units of ε . This means that the typical time scale of the environment is smaller (10 times for λ = χ = 0) than that

̸

of the qubit system. We consider three values (again in units of ε ) of the effective system-environment coupling strengths: q = 0 (no coupling), q = 0 1 . (medium coupling) and q = 1 (strong coupling). The above nomenclature concerning the coupling strength should not be taken literally as the effect of the coupling is expected to depend on the driving path. In particular, for the path across the first-order QPT both q = 0 values represent relatively strong coupling regimes (the case q = 1 may even be called an ultrastrong coupling), while for the path across the second-order QPT the case q = 0 1 corre-. sponds to a weak coupling regime. In most calculations, we set the size of the qubit system to a relatively low value, namely N = 10. In that case, the width γ of the Drude-Lorentz distribution covers the whole spectrum of the qubit system with λ = χ = 0.

At low temperatures, values of the truncation parameter up to M = 18 are needed to ensure the convergence of the HEOM calculations, while at high temperatures, the value M = 5 is found to be sufficient. For all temperatures, the number of layers is set to L = 3, the contributions of the L ≥ 4 layers being shown to yield only negligible corrections. This holds even at the low temperatures, when a short period of Re C t ( ) &lt; 0 induces a strong non-Markovian behavior.

At very high temperatures, T ≳ 25, all states of the qubit system become almost equally populated. Thus, independently of the coupling strength q and the driving time t F , the fidelity for drivings across the first- and second-order QPT approaches the values

$$\mathcal { F } _ { 1 } \underset { T \rightarrow \infty } { \longrightarrow } \frac { 1 } { N + 1 }, \quad \mathcal { F } _ { 2 } \underset { T \rightarrow \infty } { \longrightarrow } \frac { 2 } { N + 1 }, \quad ( 2 6 )$$

where ( N +1) is the dimension of the qubit Hilbert space. Therefore, we consider the fidelity for different driving times t F as a function of temperature T ≤ 25 ≈ 10 1 4 . .

Results for the firstand second-order QPT driving paths are described separately. A short comparison of the first-order QPT results obtained by the HEOM calculations with those based on the Lindblad formalism is presented afterwards.

## 4.1. First-order QPT

̸

Results of calculations of the target state fidelity F 1 from Eq. (23) for drivings across the first-order QPT are summarized in Fig. 3. The two rows of the figure correspond to driving protocols A and B (see Sec. 3.2), columns in each row present results for the three values of the effective system-environment coupling strength q given above. For both q = 0 values (the medium and strong coupling cases) we distinguish the two choices of the coupling operator ˆ from Eq. (8) by the angle Q θ . The inclusion of the counterterm in (9) is indicated by r = 1 ( r = 0 implies no addition of it). Individual curves in each plot correspond to distinct values of the driving time

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00635v4-figure-845f8cfe-c35c-4927-944e-b072f3182496.jpg)

10

Figure 3. Final fidelity of the N = 10 qubit system for drivings across the first-order QPT precursor. Results for driving protocols A and B are shown in the top and bottom rows, respectively. The effective coupling q from (10), the angle θ in (8) and the parameter r in (9) are indicated above or within each subplot. Each color and marker type represents a distinct value of the driving time t F , as indicated in the bar on the right (the markers are used just to label the curves and do not correspond to all calculated points). The range of temperature is log 10 T ∈ -[ 0 8 . , 1 4]. .

t F . These are distinguished by colors and markers that will be used also in some forthcoming figures.

Figure 3 contains rather complex information about the dependence of the fidelity on numerous variables and conditions. We start to disentangle this information by picking up some basic trends:

For q = 1, we observe a systematic increase of fidelity (compared to the corresponding r = 0 cases) for long times t F and a decrease for shorter times. This holds for both driving types. For q = 0 1 the . effect of the counterterm depends nonsystematically on the driving type and time.

- (a) Quite expectedly, slow drivings (large times t F ) lead in general to higher fidelity than fast drivings (small t F ). The difference is substantial for low temperatures T and weakens for higher T .
- (b) In absence of coupling to the environment ( q = 0), protocol B leads to higher fidelity than A.

̸

- (c) In the case of protocol A and sufficiently long times t F , the onset of interaction with the environment ( q = 0) leads to a considerable increase of fidelity. The effect is present for both medium and strong couplings, being more pronounced for θ = π 2 and r = 0. For protocol B the effect is absent, so the above-mentioned generally better performance of B does not hold for q &gt; 0.
- (d) Whereas for low and medium temperature T the fidelity depends (more or less sensitively) on all the above-mentioned variables and conditions, for high temperatures we observe convergence of all curves to the uniform limit (26).
- (e) We observe nonmonotonous dependencies of some of the F 1 curves on T for both driving types.
- (f) Effects caused by the inclusion of the counterterm ( r = 1) in Eq. (9) depend on the coupling strength:

If one seeks a way how to increase the fidelity of the target ground-state preparation, the above observations can be translated to the following practical instructions:

First, based on item (a): Cool down the system and set a large driving time. This strategy is just a straightforward extension of the common notion of adiabaticity from isolated to open systems.

Second, based on item (b): If the system is nearly isolated, use the driving protocol B instead of A Indeed, for . q = 0 the two slowest drivings with protocol B reach almost the perfect fidelity F 1 ≈ 1, while with protocol A the fidelity barely approaches to 0.54 and 0.12 for the same driving times. This instruction, which is consistent with the results of Ref. [25], follows from the reduction of the planar speed u t ( ) in the QPT region with a small ground-state energy gap, see Fig. 2, which suppresses transitions to higher-energy levels.

The third instruction is an alternative to the second one and issues from item (c): Stick with driving protocol A but let the system interact with the environment. We see that the slowest driving of type A with the S-E coupling strength q = 0 1 reaches the fidelity . F 1 ≈ 0 99 . at low temperature, comparable to driving B with q = 0, but in this case the effect survives even to medium temperatures. The increase of fidelity is present even for the strong coupling q = 1 , although in a less distinct form and also depending on the coupling operator ˆ . Q For the

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00635v4-figure-d159489a-c373-4b24-8d22-96eb83b2c4f9.jpg)

F

Figure 4. Evolving occupation probabilities P n ( ) of instantaneous energy eigenstates ( t n = 0 1 , , . . . ) of the N = 10 qubit system for driving protocol A along the path crossing the first-order QPT with t F = 10 2 8 . and T = 10 -0 04 . . The choices of q , θ and r are indicated. The bottom row depicts the whole low-lying part of the spectrum, with P n ( ) t encoded in the color of the respective level (see the color bar). The top row shows the occupation probability of the ground state. The vertical dashed lines mark the avoided crossings of the ground and excited states, which are correlated with changes of the occupation probabilities.

strong coupling, the counterterm in Eq. (9) decreases the enhancement for fast drivings but increases it for slower drivings. The increase of the target ground-state population in these cases is due to the dissipation and ongoing thermalization of the system during the drive. The same mechanism apparently does not work for the driving protocol B. In this case, the fidelity with q = 0 is already high and the onset of nonzero coupling washes out the advantage of the v = const . strategy because the S-E interaction populates excited states of the system (which makes the ground-state metric less relevant than in the q = 0 case).

To unfold the mechanism of the environment-induced increase of fidelity, we show in Fig. 4 an example of instantaneous populations of the low-lying states in the spectrum of the qubit system during the driving procedure A for some particular temperature and driving time. The occupation probability of the n th state is given by

The last and maybe most surprising strategy for increasing the target ground-state fidelity comes from the observation mentioned in the above item (e): Increase the temperature to the value for which F 1 is maximal. Indeed, many curves in Fig. 3, particularly those corresponding to smaller driving times, show a clear maximum at a certain optimal temperature T opt &gt; 0. The increase of fidelity reached in this way is not huge (and sometimes no increase is even present), but in some cases it is not negligible, like for driving A with t F = 10 2 8 . in the strong coupling case with θ = 0, when one goes from F 1 ≈ 0 2 at . low T to F 1 ≈ 0 6 at . T = T opt . While for driving protocol A this effect is present for all values of q , for protocol B it appears only for q &gt; 0.

$$P _ { n } ( t ) = \langle E _ { n } ( \boldsymbol \Lambda ( t ) ) | \hat { \rho } _ { \mathrm S } ( t ) | E _ { n } ( \boldsymbol \Lambda ( t ) ) \rangle. \quad ( 2 7 ) \quad \text{ditte}$$

So P 0 ( ) t is the instantaneous occupation of the ground state, P 1 ( ) t is the instantaneous occupation of the first excited state, and so on. For q = 0 only transitions of the Landau-Zener type modify instantaneous populations of individual states, which are initially set to the thermal values P n (0) ∝ e -βE n ( Λ I ) . We observe that in this case the ground-state occupation probability P 0 ( ) t sharply drops at the QPT-related ground-state avoided level crossing. However, for q &gt; 0, the dissipation and thermalization processes induced by interactions with the environment gradually push the ground-state occupation probability to higher values after the QPT crossing, in some cases going even above the initial population P 0 (0). This holds for both r = 0 and 1 cases.

The nature of this effect is studied in Fig. 5. It compares evolving occupation probabilities (27) of individual energy eigenstates of the qubit system for protocol A with a particular final time t F , the driving being performed at different temperatures T (this time we only show results without the counterterm). We again (as in Fig. 4) observe changes of level populations induced by transitions of the Landau-Zener type and those caused by interactions with the environment. In the top row of plots, which all correspond to T = 0, the depopulation of the ground state in the QPT region via the Landau-Zener mechanism plays an important role for the final fidelity. An interesting exception is the case with q = 1 and θ = π 2 , where energy transfers from the environment excite the system even far before the QPT. The bottom row of plots corresponds to the same driving at higher temperatures. We notice that in this case, the thermal population of excited states can sometimes increase the population of the ground state right after the QPT. This happens via the same LandauZener mechanism, which now has a partly positive effect.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00635v4-figure-42398e67-b3fa-4f95-9181-ce51b9f7b1f0.jpg)

F

Figure 5. Occupation probabilities P n ( ) t of instantaneous low-energy eigenstates of the N = 10 qubit system for driving protocol A along the path across the first-order QPT with t F = 10 2 8 . at temperatures T specified in each plot. The values of q and θ are indicated, r = 0. All plots in the top row correspond to zero temperature. In the bottom row, the temperature is set to the respective optimal value T opt for the q = 0 and q = 1, θ = 0 plots. In the other cases, for which the dependence of F 1 on T in Fig. 3 is monotonous, T is set to the same value as in the latter case.

This effect is even enhanced by transfers of populations to lower excited states when passing the ESQPT-related sequence of avoided-crossings before the QPT. In particular, for q = 0 and for q = 1, θ = 0 we set T to its respective optimal value T opt and observe an improvement of the final fidelity. However, because of complex interplay between all effects involved in the evolution of individual populations, the increase of temperature does not always lead to an advantage, as seen in Fig. 3. Indeed, for the other choices of q and ˆ in Fig. 5, for which θ F 1 in Fig. 3 monotonously decrease with T , we detect a lowering of the ground-state population at the end of the driving procedure in comparison with the T = 0 case.

Motivated by the experimental work in Ref. [46], which analyzed how the maximal temperature needed to reach fidelity at or above a certain limiting value scales with the size of the system, we further investigate the dependence of the above-discussed temperature-driven effect on the number N of qubits. In Fig. 6 we show the N -dependence of the optimal temperature T opt and of the locally maximal value of fidelity max F 1 at this temperature. We use protocol A with different driving times t F and set q = 0 for simplicity. We see that the optimal temperature grows approximately linearly with N , hence T opt = aN with the coefficient a depending on t F . On the other hand, the local maximum of fidelity at T opt drops algebraically, max F 1 = bN -κ with b depending on t F and κ ≈ 1. It is obvious that this finding strongly reduces the applicability of the last fidelity-increase strategy in systems with a larger number of qubits.

In any case, nonmonotonous dependencies of fidelity on temperature have already been detected experimentally in Ref. [45]. Although the mechanism of the thermally

Figure 6. Optimal temperature T opt (left panel) and the fidelity max F 1 obtained for this temperature (right panel) as functions of the size parameter N for driving protocol A across the first-order QPT with q = 0. The markers represent numerical data for different driving times t F (we use the same marker colours and types as in Fig. 3), while the dashed lines are linear fits.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00635v4-figure-0f2cbc16-db86-4f33-af8d-37b1f8433ebe.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00635v4-figure-49e17243-4c77-4251-85ca-4b47411f3f51.jpg)

assisted quantum annealing discussed in that study is basically the same as the mechanism of the present effect, some aspects of both systems are different. In particular, the setup of Ref. [45] involved a single level crossing including only the ground and first excited states, while the higher energy eigenstates were just spectators. In that situation, it was possible to derive an analytic estimate of the optimal temperature T opt . In our case, however, that estimate does not work because of the more complicated ESQPT-related structure of avoided crossings.

Let us finally note that our present results (see particularly the T = 0 row of Fig. 5) can be compared with the previous studies [51, 84] of driving of a dissipative singlequbit (two-level) model through a single avoided level crossing at zero temperature. The system-environment coupling operators employed in these studies were identical with our choice, namely ˆ Q ∝ ˆ σ z and ˆ Q ∝ ˆ σ x . Although some of our findings are compatible with these older results (in particular, the increase of the groundstate population with the onset of dissipation [51]), our present study demonstrates that the N &gt; 1 qubit system yields more complex dependencies than the simplest case of N = 1.

## 4.2. Second-order QPT

Figure 7 displays the target ground-state fidelity F 2 from Eq. (24) for the path across the second-order QPT. The plots in this figure show all dependencies as in the previous case of the first-order QPT (Sec. 4.1), i.e., the dependence on the temperature T , the S-E coupling strength q , the angle θ of the coupling operator (8), the driving time t F , the parameter r in (9), and on the driving protocols A and B. Since the minimal ground-state energy gap ∆ 10 along second-order QPT path decreases with N algebraically, the driving times t F needed to

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00635v4-figure-f4353440-79ec-4bef-81e8-29cbba69943f.jpg)

10

Figure 7. Same as in Fig. 3 but for the path across the second-order QPT and smaller driving times.

reach values of fidelity comparable to the first-order QPT path (for which the gap closes exponentially) are much shorter than in that case. So we move from the time scale t F ∈ [10 0 4 . , 10 3 6 . ] used for drivings across the firstorder QPT to t F ∈ [10 0 , 10 1 6 . ] used for drivings across the second-order QPT.

T opt and a roughly algebraic decrease of max F 2 with the size N of the qubit system.

As in the case of the first-order QPT drivings, we try to itemize the immediate observations following from Fig. 7:

- (a) Trivially, large t F and small T imply higher fidelity than small t F or large T .
- (b) For zero or weak coupling q the driving protocol B still yields higher fidelity than A, but the effect is visible only for large t F .
- (c) The effect of environment-induced increase of fidelity is more or less gone. Interactions with the environment start playing some role only for strong coupling (as expected) and they mostly decrease the fidelity.
- (d) The high-temperature limit (26) of fidelity is valid.
- (e) Nonmonotonous dependencies of some fidelity curves on the temperature are again observed.
- (f) The counterterm in the Hamiltonian improves the fidelity for the strong coupling (compared to the r = 0 case).
- (g) For the strong coupling, q = 1, we observe a strong dependence of F 2 on the coupling operator ˆ = Q ˆ J x (i.e., θ = π 2 ) and ˆ J z (i.e., θ = 0).

The peak of some fidelity curves, see item (e), has been subject to the same analysis as in the first-order QPT case, the result being summarized in Fig. 8. We again observe the linear increase of the optimal temperature

In the case of the second-order QPT path, special attention needs to be paid to the choice of the systemenvironment coupling operator ˆ in Eq. (8). Q As pointed out in item (g), this choice makes a large difference in the values of F 2 obtained for the strong coupling. This effect is due to the difference between both choices in terms of the parity conservation. We know that for q = 0 the driving from the positive-parity initial state along the χ = 0 path can excite the system only to states with positive parity. However, for q = 0 the parity may be violated by system-environment interaction. This happens for θ = π 2 , so during the drive the environment induces strong transitions to the first excited state with negative parity. This population is then counted in the summed final fidelity F 2 in Eq. (24), so for slow-enough driving (when transitions to higher excited states can be neglected) this fidelity may be relatively high.

̸

In contrast, for θ = 0 the parity is conserved even by the system-environment interaction, so the strongest transitions during the drive lead to the second excited state with positive parity, which does not contribute to F 2 . This explains why slow drivings in the strongcoupling columns of Fig. 7 yield higher fidelity for θ = π 2 than for θ = 0. In the latter case all F 2 curves converge to a narrow band, washing out the difference between slow and fast driving. It seems that for θ = 0 the environmentinduced transitions and the transitions of the LandauZener type have similar effects which complement each other as the driving time t F varies.

## 4.3. Comparison with the Lindblad method

In this section we compare the above HEOM results with simulations using the Lindblad method. The use of

Figure 8. Same as Fig. 6 but for the second-order QPT path. The dashed lines represent dependencies T opt = aN and max F 2 ≈ b/N -c/N 2 with fitted coefficients a, b, c .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00635v4-figure-d1e3e8be-04b3-4f5f-8ced-c1d404632f97.jpg)

the Lindblad formalism for time-dependent Hamiltonians is discussed, e.g., in Refs. [50-52]. In our calculations we apply the Lindblad formula

$$\frac { d } { d t } \hat { \rho } _ { S } ( t ) = - i [ \hat { H } _ { S } ( t ) + \hat { H } _ { \text{L} } ( t ), \hat { \rho } _ { S } ( t ) ] + \mathcal { D } _ { t } \hat { \rho } _ { S } ( t ) \quad ( 2 8 ) \quad t = 2 \text{ \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

with the dissipator,

$$\mathcal { D } _ { t } \hat { \rho } _ { S } ( t ) = & \sum _ { \epsilon } \frac { \Gamma ( \epsilon ) } { N } \left [ \hat { S } ( \epsilon ) \hat { \rho } _ { S } ( t ) \hat { S } ^ { \dagger } ( \epsilon ) - \frac { 1 } { 2 } \left \{ \hat { S } ^ { \dagger } ( \epsilon ) \hat { S } ( \epsilon ), \hat { \rho } _ { S } ( t ) \right \} \right ]. \quad \underset { ( 2 9 ) } { \text{Indb} } \quad \text{In F} \\ & ( 2 9 ) \quad \text{by bo}$$

Here the sum goes over all possible differences ϵ = E m ( Λ ( )) t -E n ( Λ ( )) between instantaneous system ent ergies and

$$\hat { S } ( \epsilon ) \ = \sum _ { \substack { m, n \\ E _ { m } - E _ { n } = \epsilon \\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \colon \\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cd$$

$$\Gamma ( \epsilon ) \ = \ 2 \left [ 1 + n _ { E } ( \epsilon ) \right ] \left [ J ( \epsilon ) \Theta ( \epsilon ) - J ( - \epsilon ) \Theta ( - \epsilon ) \right ], \ ( 3 1 ) \quad \text{sauce} \\ \text{and} \, 1.$$

where Θ denotes the Heaviside step function. The function (31) is positive and satisfies Γ(0) = 4 qT/γ for the Drude-Lorentz form (13). Negative or positive values of ϵ indicate excitation or relaxation of the system, respectively, while ϵ = 0 indicates dephasing [47]. The hermitian term

$$\hat { H } _ { \text{L} } ( t ) = \frac { 1 } { N } \sum _ { \epsilon } \left [ \frac { c _ { 0 } \epsilon } { \gamma ^ { 2 } + \epsilon ^ { 2 } } + \sum _ { k = 1 } ^ { \infty } \frac { c _ { k } \epsilon } { ( 2 \pi k T ) ^ { 2 } + \epsilon ^ { 2 } } \right ] \hat { S } ^ { \dagger } ( \epsilon ) \hat { S } ( \epsilon ) \quad \text{This} \quad \\ \text{(32)} \quad \text{in gen} \quad \\ \text{in } (98)$ represents the s o called Lamb shift  Haimiltonian \quad \text{the e}$}$$

in (28) represents the so-called Lamb shift Hamiltonian, with coefficients c 0 , c 1 , . . . introduced in (25).

In the comparison of the Lindblad calculations with the HEOM results, we consider only the path across the first-order QPT with the driving protocol A (see Sec. 4.1). The S-E coupling operator (8) has θ = 0 and the counterterm in (9) is not included ( r = 0). If the dynamics for N = 10 qubits is modeled using the two-state approximation, the adiabatic timescale is of the order

Figure 9. Fidelity for driving A across the first-order QPT (with θ = 0 and r = 0) for the N = 10 qubit system obtained by the HEOM method (full curves) and by the Lindblad method according to Ref. [51] (dashed curves). The driving times t F and coupling strengths q are indicated, the range of temperature is log 10 T ∈ -[ 0 8 . , 1 4]. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00635v4-figure-1aa698c7-b8a2-460f-a7d1-145261f01e30.jpg)

t = 2 | d ∆ 10 /ds /π | ∆ 2 10 ∼ 10 , 3 where all terms are evaluated at the crossing [45]. The method of Ref. [50, 51] is limited to nearly-adiabatic drivings, so only the drivings with t F ≫ 10 3 can be expected to agree with the Lindblad calculations.

In Fig. 9 we show the fidelity F 1 from Eq. (23) obtained by both the HEOM and Lindblad methods for driving times t F = 10 2 8 . and 10 3 6 . . We see, in accord with the above expectation, that for the larger t F the discrepancy between both methods affects only the drivings at lower temperatures ( T ≲ γ/ π 2 , see Ref. [63]) and shows up only for the strong S-E coupling. For this driving time, the Lindblad approach predicts the final state of the system that essentially coincides with the canonical equilibrium state ˆ ρ ( β ) S ( Λ F ) ∝ e -βH ˆ S (Λ F ) . The discord of the HEOM calculations with this prediction for the stronger coupling and low temperatures agrees with the results of Refs. [5456], where deviations of the actual final state from the canonical state were also observed. Note that a realistic approximation for the equilibrium state for all coupling strengths was derived in Refs. [58, 59]. For the shorter t F in Fig. 9, the HEOM and Lindblad results differ for both coupling strengths and all temperatures below the highT limit (26).

Therefore, we can conclude that non-Markovian effects, properly captured in our HEOM calculations, play in general an important role, which is not reproduced by the simpler Lindblad calculations.

## 5. CONCLUSION

In this paper we studied effects of thermalization and dissipation in quantum driven dynamics across finite-size precursors of QPTs of the first and second order. Using a fully connected system of qubits coupled to an external

environment at a nonzero temperature, we tested various strategies to maximize the fidelity of finite-time preparation of the pure target state associated with the entangled ground state of the qubit system at the final parameter point behind the QPT.

In agreement with our previous studies of driving in the same system without environment-induced effects [25, 26] we found that the final fidelity for a fixed finite value of the driving time can be increased when one applies the driving with constant geometric instead of planar speed. This improvement is attributed to the suppression of exciting transitions in the QPT region due to the reduced planar speed in that region.

However, thermalization and dissipation processes can strongly modify results of the driving procedure, in some cases being able to considerably improve its fidelity. One of the effects results from the environment-induced thermalization of the qubit system in the course of the driving procedure, which under some circumstances increases population of the target ground state. Another effect follows from an initial thermalization of the system at the start of the driving procedure, which for a certain optimal temperature leads to an efficient transfer of populations from excited states to the ground state in the QPT and ESQPT regions. Both effects are present for drivings across the first- as well as second-order QPT. The

- [1] L.K. Grover, Synthesis of Quantum Superpositions by Quantum Computation , Phys. Rev. Lett. 85 , 1334 (2000).
- [2] M. Plesch and C. Brukner, ˇ Quantum-state preparation with universal gate decompositions , Phys. Rev. A 83 , 032302 (2011).
- [3] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, Quantum Computation by Adiabatic Evolution , arXiv: quant-ph/0001106 (2000).
- [4] D. Aharonov and A. Ta-Shma, Adiabatic Quantum State Generation , SIAM J. Comput. 37, 47 (2007).
- [5] D. Aharonov, W. van Dam, J. Kempe, Z. Landau, S. Lloyd, and O. Regev, Adiabatic Quantum Computation is Equivalent to Standard Quantum Computation , SIAM J. Comput. 37, 166 (2007).
- [6] T. Albash and D. A. Lidar, Adiabatic quantum computation , Rev. Mod. Phys. 90 , 015002 (2018).
- [7] M. A. Nielsen and I. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, Cambridge, UK, 2010).
- [8] T. Kato, On the Adiabatic Theorem of Quantum Mechanics , J. Phys. Soc. Jpn. 5 , 435 (1950).
- [9] S. Jansen, M.-B. Ruskai, R. Seiler, Bounds for the adiabatic approximation with applications to quantum computation , J. Math. Phys. (N.Y.) 48, 102111 (2007).
- [10] A. Elgart, and G.A. Hagedorn, A note on the switching adiabatic theorem , J. Math. Phys. (N.Y.) 53, 102202 (2012).
- [11] M. Demirplak and S. A. Rice, Adiabatic Population Transfer with Control Fields , J. Chem. Phys. 107 , 9937 (2003).

latter effect was nevertheless found to weaken with an increasing size of the system, which limits its application in driving procedures involving large numbers of qubits.

A by-product of our analysis is a comparison of the sophisticated HEOM calculations of the systemenvironment dynamics with much simpler and more popular Lindblad calculations. We have demonstrated that in the present setup of driving across the finite-size QPT precursors the Lindblad calculations, which completely disregard non-Markovian effects, do not satisfactorily reproduce the more precise HEOM calculations for stronger system-environment couplings or for faster drivings.

Finally, it needs to be stressed that all the above results were obtained in the framework of our strongly simplified model of the interacting qubit system and its environment. Nevertheless, we consider the effects analyzed here as sufficiently robust, at least on the qualitative level, to play important roles in realistic situations, whose quantitative analysis may require additional calculations with modified model assumptions.

## ACKNOWLEDGEMENTS

We thank Jakub Dolejˇ s´ ı for his help in the initial stages of this work within his diploma thesis [85].

- [12] M. V. Berry, Transitionless quantum driving , J. Phys. A 42 , 365303 (2009).
- [13] A. del Campo, Shortcuts to Adiabaticity by Counterdiabatic Driving , Phys. Rev. Lett. 111 , 100502 (2013).
- [14] D. Gu´ ery-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Mart´ ınez-Garaot, and J.G. Muga, Shortcuts to adiabaticity: Concepts, methods, and applications , Rev. Mod. Phys. 91 , 045001 (2019).
- [15] J. P. Provost and G. Vallee, Riemannian structure on manifolds of quantum states , Commun. Math. Phys. 76 , 289 (1980).
- [16] J. Anandan, Y. Aharonov, Geometry of quantum evolution , Phys. Rev. Lett. 65 , 1697 (1990).
- [17] A. Miyake, M. Wadati, Geometric strategy for the optimal quantum search , Phys. Rev. A 64 , 042317 (2001).
- [18] A. T. Rezakhani, W.-J. Kuo, A. Hamma, D. A. Lidar, and P. Zanardi, Quantum Adiabatic Brachistochrone , Phys. Rev. Lett. 103 , 080502 (2009).
- [19] M. Kolodrubetz, D. Sels, P. Mehta, and A. Polkovnikov, Geometry and non-adiabatic response in quantum and classical systems , Phys. Rep. 697 , 1 (2017).
- [20] M. Bukov, D. Sels, and A. Polkovnikov, Geometric Speed Limit of Accessible Many-Body State Preparation , Phys. Rev. X 9 , 011034 (2019).
- [21] A.M. Childs, E. Deotto, E. Farhi, J.Goldstone, S. Gutmann, and A.J. Landahl, Quantum search by measurement , Phys. Rev. A 66 , 032314 (2002).
- [22] L. Roa, A. Delgado, M.L. Ladr´nde Guevara, and A.B. o Klimov, Measurement-driven quantum evolution , Phys. Rev. A 73 , 012322 (2006).

- [23] S. Hacohen-Gourgy, L.P. Garc´ ıa-Pintos, L.S. Martin, J. Dressel, and I. Siddiqi, Incoherent Qubit Control Using the Quantum Zeno Effect , Phys. Rev. Lett. 120 , 020505 (2018).
- [24] P. Cejnar, P. Str´nsk´, a y J. Stˇ releˇek, c and F. Matus, Decoherence-assisted quantum driving , Phys. Rev. A 107 , L030603 (2023).
- [25] F. Matus, J. Stˇeleˇek, P. Str´nsk´, and P. Cejnar, r c a y Search for optimal driving in finite quantum systems with precursors of criticality , Phys. Rev. A 107 , 012216 (2023).
- [26] F. Matus, J. Stˇ releˇek, c and P. Cejnar, Analytic approach to the Landau-Zener problem in bounded parameter space , J. Phys. A: Math. Theor. 56 , 235303 (2023).
- [27] H. J. Lipkin, N. Meshkov, and A. J. Glick, Validity of many-body approximation methods for a solvable model: (I). Exact solutions and perturbation theory , Nucl. Phys. 62 , 188 (1965).
- [28] R. Gilmore and D. H. Feng, Phase transitions in nuclear matter described by pseudospin Hamiltonians , Nucl. Phys. A 301 , 189 (1978).
- [29] R. Or´s, S. Dusuel, and J. Vidal, u Equivalence of Critical Scaling Laws for Many-Body Entanglement in the LipkinMeshkov-Glick Model , Phys. Rev. Lett. 101 , 025701 (2008).
- [30] T. Zibold, E. Nicklas, C. Gross, and M.K. Oberthaler, Classical Bifurcation at the Transition from Rabi to Josephson Dynamics , Phys. Rev. Lett. 105 , 204101 (2010).
- [31] S. Puri, C. K. Andersen, A. L. Grimsmo, and A. Blais, Quantum annealing with all-to-all connected nonlinear oscillators , Nat. Commun. 8 , 15785 (2017).
- [32] M. J. Cervia, A. B. Balantekin, S. N. Coppersmith, C. W. Johnson, P. J. Love, C. Poole, K. Robbins, and M. Saffman, Lipkin model on a quantum computer , Phys. Rev. C 104 , 024305 (2021).
- [33] S. Sachdev, Quantum Phase Transitions (Cambridge University Press, Cambridge, UK, 1999).
- [34] J. I. Latorre and R. Or´s, u Adiabatic quantum computation and quantum phase transitions , Phys. Rev. A 69 , 062302 (2004).
- [35] W.H. Zurek, U. Dorner, and P. Zoller, Dynamics of a Quantum Phase Transition , Phys. Rev. Lett. 95 , 105701 (2005).
- [36] R. Sch¨tzhold and G. Schaller, u Adiabatic quantum algorithms as quantum phase transitions: First versus second order , Phys. Rev. A 74 , 060304(R) (2006).
- [37] P. Cejnar, M. Macek, S. Heinze, J. Jolie, and J. Dobeˇ s, Monodromy and excited-state quantum phase transitions in integrable systems: collective vibrations of nuclei , J. Phys. A: Math. Gen. 39 , L515 (2006).
- [38] P. Ribeiro, J. Vidal, and R. Mosseri, Thermodynamical Limit of the Lipkin-Meshkov-Glick Model , Phys. Rev. Lett. 99 , 050402 (2007).
- [39] M. Caprio, P. Cejnar, and F. Iachello, Excited state quantum phase transitions in many-body systems , Ann. Phys. 323 , 1106 (2008).
- [40] P. Cejnar, P. Str´nsk´, a y Impact of quantum phase transitions on excited-level dynamics , Phys. Rev. E 78 , 031130 (2008).
- [41] P. P´ erez-Fern´ndez, a A. Rela˜o, n J. M. Arias, J. Dukelsky, and J. E. Garc´ ıa-Ramos, Decoherence due to an excited-state quantum phase transition in a two-level boson model , Phys. Rev. A 80 , 032111 (2009).
- [42] L. F. Santos, M. T´vora and F. P´ erez-Bernal, a Excitedstate quantum phase transitions in many-body systems with infinite-range interaction: Localization, dynamics, and bifurcation , Phys. Rev. A 94 , 012113 (2016).
- [43] W. Kopylov, G. Schaller, and T. Brandes, Nonadiabatic dynamics of the excited states for the Lipkin-MeshkovGlick model , Phys. Rev. E 96 , 012153 (2017).
- [44] P. Cejnar, P. Str´nsk´, M. Macek, M. Kloc, a y Excited-state quantum phase transitions , J. Phys. A: Math. Theo. 54 , 133001 (2021).
- [45] N. Dickson, M. Johnson, M. Amin et al , Thermally assisted quantum annealing of a 16-qubit problem , Nat. Commun. 4 , 1903 (2013).
- [46] T. Albash, V. Martin-Mayor and I. Hen, Temperature Scaling Law for Quantum Annealing Optimizers , Phys. Rev. Lett. 119 , 110502 (2017).
- [47] H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Clarendon Press, Oxford, 2002).
- [48] U. Weiss, Quantum Dissipative Systems (World Scientific, 4th ed., Singapore, 2012).
- [49] M. Schlosshauer, Decoherence and the Quantum-toclassical transition (Springer, 1st ed., Berlin, 2010).
- [50] T. Albash, S. Boixo, D. A. Lidar, and P. Zanardi, Quantum adiabatic Markovian master equations , New J. Phys. 14 , 123016 (2012).
- [51] M. Yamaguchi, T. Yuge and T. Ogawa, Markovian quantum master equation beyond adiabatic regime , Phys. Rev. E 95 , 012136 (2017).
- [52] R. Dann, A. Levy, and R. Kosloff, Time-dependent Markovian quantum master equation , Phys. Rev. A 98 , 052129 (2018).
- [53] A. Ishizaki1 and Y., Tanimura, Quantum Dynamics of System Strongly Coupled to Low-Temperature Colored Noise Bath: Reduced Hierarchy Equations Approach , J. Phys. Soc. Jpn. 74 , 3131-3134 (2005).
- [54] A. G. Dijkstra and Y. Tanimura, System Bath Correlations and the Nonlinear Response of Qubits , J. Phys. Soc. Jpn. 81 , 063301 (2012).
- [55] J. Iles-Smith, N. Lambert, and A. Nazir, Environmental dynamics, correlations, and the emergence of noncanonical equilibrium states in open quantum systems , Phys. Rev. A 90 , 032114 (2014).
- [56] J. Iles-Smith, A.G. Dijkstra, N. Lambert, and A. Nazir, Energy transfer in structured and unstructured environments: Master equations beyond the Born-Markov approximations , J. Chem. Phys. 144 , 044110 (2016).
- [57] J. W. Abbott, Quantum Dynamics of Bath Influenced Excitonic Energy Transfer in Photosynthetic PigmentProtein Complexes Creators , Master Thesis (2020)
- [58] N. Anto-Sztrikacs, A. Nazir, and D. Segal, Effective Hamiltonian theory of open quantum systems at strong coupling , PRX Quantum 4 , 020307 (2023)
- [59] N. Anto-Sztrikacs, B. Min, M. Brenes and D. Segal, Effective Hamiltonian theory: An approximation to the equilibrium state of open quantum systems , Phys. Rev. B 108 , 115437 (2023)
- [60] Y. Tanimura and R. Kubo, Time Evolution of a Quantum System in Contact with a Nearly Gaussian-Markoffian Noise Bath , J. Phys. Soc. Jpn. 58 , 101 (1989).
- [61] Y. Tanimura, Stochastic Liouville, Langevin, Fokker-Planck, and Master Equation Approaches to Quantum Dissipative Systems , J. Phys. Soc. Jpn. 75 , 082001 (2006).

- [62] A. Fruchtman, N. Lambert, and E.M. Gauger, When do perturbative approaches accurately capture the dynamics of complex quantum systems? , Sci. Rep. 6 , 28204 (2016).
- [63] Y. Tanimura, Numerically 'exact' approach to open quantum dynamics: The hierarchical equations of motion (HEOM) , J. Chem. Phys. 153 , 020901 (2020).
- [64] N. Lambert, T. Raheja, S. Cross, P. Menczel, S. Ahmed, A. Pitchford, D. Burgarth, and F. Nori, QuTiPBoFiN: A bosonic and fermionic numerical hierarchicalequations-of-motion library with applications in lightharvesting, quantum control, and single-molecule electronics , Phys. Rev. Research 5 , 013181 (2023).
- [65] R. P. Feynman and F. L. Vernon Jr., The theory of a general quantum system interacting with a linear dissipative system , Ann. Phys. (USA) 24 , 118 (1963).
- [66] A. O. Caldeira and A. J. Leggett, Influence of Dissipation on Quantum Tunneling in Macroscopic Systems , Phys. Rev. Lett. 46 , 211 (1981).
- [67] A. O. Caldeira and A. J. Leggett, Quantum tunnelling in a dissipative system , Ann. Phys. (NY) 149 , 374 (1983).
- [68] J. Xu, R.-X. Xu, and Y.J. Yan, Exact quantum dissipative dynamics under external time-dependent driving fields , New J. Phys. 11 , 105037 (2009).
- [69] J. Ma, Z. Sun, X. Wang, and F. Nori, Entanglement dynamics of two qubits in a common bath , Phys. Rev. A 85 , 062323 (2012) .
- [70] A. Ishizaki and G. R. Fleming, Theoretical examination of quantum coherence in a photosynthetic system at physiological temperature , PNAS 106 , 17255 (2009).
- [71] L. Chen, R. Zheng, Q. Shi, and Y. Yan, Optical line shapes of molecular aggregates: Hierarchical equations of motion method , J. Chem. Phys. 131 , 094502 (2009).
- [72] J. Str¨mpfer and K. Schulten, u The effect of correlated bath fluctuations on exciton transfer , J. Chem. Phys. 134 , 095102 (2011).
- [73] N. Lambert, S. Ahmed, M. Cirio, and F. Nori, Modelling the ultra-strongly coupled spin-boson model with unphys-
- ical modes , Nat. Commun. 10 , 3721 (2019).
- [74] P. Cejnar and P. Str´nsk´, a y Quantum phase transitions in the collective degrees of freedom: nuclei and other manybody systems , Phys. Scr. 91 , 083006 (2016).
- [75] J. Stˇeleˇek and P. Cejnar, r c Quantum geometry in manybody systems with precursors of criticality , arXiv: quantph/2408.00635 (2024).
- [76] K. Rza˙ zewski, K. W´dkiewicz, and W. Zakowicz, o ˙ Phase Transitions, Two-Level Atoms, and the A 2 Term , Phys. Rev. Lett. 35 , 432 (1975).
- [77] K. Hepp and E. H. Lieb, On the superradiant phase transition for molecules in a quantized radiation field: the Dicke maser model , Ann. Phys. 76 , 360 (1973).
- [78] Y. K. Wang and F. T. Hioe, Phase Transition in the Dicke Model of Superradiance , Phys. Rev. A 7 , 831 (1973).
- [79] C. Emary and T. Brandes, Quantum Chaos Triggered by Precursors of a Quantum Phase Transition: The Dicke Model , Phys. Rev. Lett. 90 , 044101 (2003).
- [80] K. Baumann, C. Guerlin, F. Brennecke and T. Esslinger, Dicke quantum phase transition with a superfluid gas in an optical cavity , Nature 464 , 1301 (2010).
- [81] P. Nataf and C. Ciuti, No-go theorem for superradiant quantum phase transitions in cavity QED and counterexample in circuit QED , Nat Commun 1 , 72 (2010).
- [82] W. Kopylov and G. Schaller, Polaron-transformed dissipative Lipkin-Meshkov-Glick model , Phys. Rev. A 100 , 063815 (2019).
- [83] J. Roland and N. J. Cerf, Quantum search by local adiabatic evolution , Phys. Rev. A 65 , 042308 (2002).
- [84] M. Wubs, K. Saito, S. Kohler, P. H¨ anggi and Y. Kayanuma, Gauging a Quantum Heat Bath with Dissipative Landau-Zener Transitions , Phys. Rev. Lett. 97 , 200404 (2006).
- [85] J. Dolejˇ ´, Diploma Thesis (Faculty of Mathematics and sı Physics, Charles University, 2020).
Felipe Matus Pavel Cejnar
Driving of an open quantum system at finite temperature across first- and second-order quantum phase transitions
An open fully connected system of qubits at nonzero temperature is driven
within a finite time interval along various paths in the space of its control
parameters. The driving leads across finite-size precursors of first- and
second-order quantum phase transition from factorized to entangled ground-state
phases, aiming at the preparation of the complex ground state of the system at
the final parameter point with maximal fidelity. During the drive, the system
is coupled to a heat bath at a constant temperature, the dynamics being
determined in a nonperturbative way by the method of Hierarchical Equations of
Motion. It is shown that the presence of the heat bath in combination with
specific patterns of avoided crossings affecting the ground and excited states
in the parameter region around the quantum phase transition may considerably
improve the fidelity of preparation of the target ground state.