Trace preserving quantum dynamics using a novel reparametrization-neutral summation-by-parts difference operator

.


Introduction
The dynamical evolution of dissipative quantum systems constitutes the focus of a broad and active research community, ranging from quantum information science to condensed matter physics and even to high energy nuclear physics. Take the physics of impurities for example, where a heavy particle is introduced into a reservoir of light particles. The interactions of the impurity with the surrounding environment change its properties, leading to a wealth of phenomenologically relevant quantum phenomena, as well as decoherence.
In the context of condensed matter physics such impurities have been discussed as Bose polarons (see e.g. [1]). On the other hand in the context of heavy-ion collisions it is the behavior of heavy quark-antiquark pairs, so called heavy quarkonium and their interaction with the hot collisions remnants that are of interest (for a recent review see [2]). Investigating the physics of the small subsystem inevitably leads one to consider dissipative dynamics and to the framework of open quantum systems (for an excellent introduction see [3]).
It is common to consider a quantum system, which consists of a small subsystem (S), immersed in a large environment (E). Its state is encoded in the eponymous state-vectors |ψ that are elements of the full Hilbert space H. The subsystem and environment degrees of freedom, described by |ψ S and |ψ E respectively, can be found in separate subspaces H = H S ⊗H E . The symbol ⊗ denotes the direct product. The overall system is closed and its physics described by the hermitian Hamilton operatorĤ tot acting on state vectors in the full Hilbert space. It can be decomposed into the following sumĤ tot =Ĥ S ⊗ I E + I S ⊗Ĥ E +Ĥ int . (1) The first termĤ S refers to the Hamiltonian governing the dynamics of the subsystem. It leaves the environment unchanged, as indicated by the direct product with the identity operator I E on the environment subspace of the full Hilbert space.Ĥ E describes the environment degrees of freedom and the interactions between the two are contained inĤ int . The latter explicitly couples the subsystem and environment subspace of the full Hilbert space. Many relevant properties of such a system are captured by its density matrix operatorρ where we denote the adjoint of a state vector by |ψ i † = ψ i |. The density matrix thus represents the outer product of all states realized in the system under consideration, weighted by their probability p i . In the so-called Schrödinger pictureρ tot evolves according to the von-Neumann equation where [A, B] = AB − BA denotes the commutator. The hermiticity ofĤ tot translates into unitary time evolution for the total density matrix ρ tot (t) = U (t, 0)ρ tot (0)U † (t, 0) implemented via the time evolution operator U (t, 0) = T exp[−i dtH int (t)]. Here the adjoint of the operator is denoted by the ( †) symbol U † (t, 0) = U −1 (t, 0) = U (0, t) and the symbol T refers to time ordering, relevant for explicitly time dependent Hamilton operators (for a more detailed exposition we refer the reader to [4]).
Measurable properties of the quantum system are encoded in hermitian operators (Â † =Â), so called observables. Their expectation value, representing the mean value of the associated physical property, obtained over repeated experiments, may be computed via trace over the density matrix Â = Tr [ρÂ]. For observables related to the small subsystem, we now wish to simplify the description. Instead of having to carry out the trace over the full Hilbert space every time we compute an expectation value, we carry out the trace over the environment degrees of freedom in the full Hilbert space a priori. This leads us to the reduced density matrix where we have expressed the trace as sum over the inner products with the n-th environment state vector ψ E n |ρ|ψ E n . Hence, in case thatÂ =Â S ⊗ I E , we have Â (t) = Tr[ρ(t)Â] = Tr S [ρ S (t)Â S ]. We thus need an evolution equation for the reduced density matrixρ S .
In general, the time scales of the evolution in the subsystem and the environment may be of the same order and the time evolution ofρ S remains genuinely non-Markovian, i.e. memory effects play a role. On the other hand if there exists a separation of timescales between the medium and the subsystem, i.e. if the subsystem takes longer to relax than the environment correlations take to decay, one may encounter approximate Markovian dynamics, in which memory effects can be neglected. In that case it has been shown that the most general equation of motion forρ S can be written in terms of the (Gorini-Kossakowski-Sudarshan)-Lindblad equation [5,6] This equation describes the in-general non-unitary, i.e. dissipative, time evolution of the subsystem S through its reduced density matrixρ S . I.e. all operators here act only on states |ψ S ∈ H S in the subsystem subspace of the full Hilbert space. The influence of the environment E manifests itself in the presence of Lindblad operatorsL k , damping rates γ k and possible modifications of the subsystem HamiltonianH S =Ĥ S (explicit examples of these quantities will be derived in Section 2). Formulating the dissipative dynamics in terms of a Lindblad equation is advantageous, as it can be proven that this equation preserves the main physical properties of the reduced density matrix, i.e. positivity, hermiticity and unitarity The preservation of unit trace in particular is important, as the probability interpretation of the density matrix rests upon it. The reason is that understanding the energy, momentum and particle exchange between the subsystem and environment is of central interest in our study. Thus it is paramount to ensure that the formulation of the problem does not introduce artificial loss channels. These may e.g. deplete the probabilitiesp l in Eq. (4) beyond the true effect induced by the presence of the environment. In general Eq. (5) can be expressed in the coordinate space basis of the Hilbert space, where the matrix elements of the reduced density matrix ρ(x 1 , x 2 , . . . , y 2 , y 1 , t) = x 1 , x 2 , . . . |ρ S (t)| . . . , y 2 , y 1 ∈ C form a complex function with a dependence on twice the number of coordinates x i , y i ∈ R 3 as are particles present in the system. The ensuing partial differential equation remains linear in the density matrix, but contains both spatially varying and complex valued coefficients and derivative terms for each coordinate = F x 1 , x 2 , . . . , y 2 , y 1 , ∇ x1 , ∇ x2 , . . . , ∇ y2 , ∇ y1 , t ρ(x 1 , x 2 , . . . , y 2 , y 1 , t).
The computational challenge, which we address in this paper, lies in discretizing Eq. (7) and implementing it with a stable and accurate numerical procedure, which in addition guarantees that the properties in Eq. (6) are preserved. Using an arbitrary complex test function f (x 1 , x 2 , . . .) ∈ C and denoting by δ (3) (x 1 − y 1 ) the three-dimensional delta function, these properties can be formulated in terms of the matrix elements as follows.
Unit trace : The numerical treatment of initial-boundary value problems (IBVPs), among them the Navier-Stokes and Schrödinger-like equations, such as Eq. (7), has seen significant progress over the past decade with the development and refinement of summation-by-parts (SBP) difference operators (for reviews see e.g. [7,8,9]). As these operators build upon the finite difference approach (although they can be formulated for many other schemes, see [10,11,12,13,14,15,16,17,18]) they are straight forward to implement and their numerical evaluation cost is low. The fact that they mimic the integration by parts property of the continuum theory facilitates proofs of stability, e.g. when deploying SBP operators in time stepping approaches for computational fluid dynamics [19]. After the development of SBP operators for first derivatives, higher derivative approximations [20,21] have been derived. More recently the SBP technique has also been applied to derivatives in time direction [9,22,23]. While in this study only periodic boundary conditions will be deployed, the SBP operators can easily accommodate non-trivial boundary conditions (in the weak sense) via the Simultaneous Approximation Term (SAT) technique [24].
To make the paper self-contained, we provide a brief introduction to SBP operators and recommend [7,8] for extensive reviews. Let the domain [x L , x R ] be discretized with N +1 equidistant grid points x i = x L +i∆x, i = 0, . . . , N , where ∆x = (x R − x L )/N . Denote by u(t) = [u 0 , . . . , u N ] the vector containing the function u(t, x) evaluated at spatial grid points at time t. The approximation of the spatial derivative is given by where u x contains the analytical derivative evaluated on the grid. For two functions u, v defined on the grid, we have where the matrix H is diagonal, positive definite and defines an inner product and a corresponding norm. Furthermore, the differentiation operator D satisfies the SBP property where In the second order case, H is the composite trapezoidal rule and D is the standard stencil for the symmetric central difference in the interior and appropriate forward and backward stencils at the boundaries: In the periodic case, the operators simplify to The SBP property already provides a crucial ingredient in the formulation of stable approximations of IBVPs, such as Eq. (7). In order to also preserve the trace of a relevant class of Lindblad equations, we will show that another continuum property of derivatives needs to be fulfilled: reparametrization neutrality.
Reparametrization neutrality refers to the fact that among a set of derivatives with respect to different variables, we may freely change to derivatives expressed in linear combinations of these variables. As a concrete example take x and y and the corresponding derivatives d dx and d dy . Considering instead z = x − y and z = x + y we may reexpress For the derivation of the trace conservation in the Lindblad equation we will need to switch from expressions in x and y to expressions in z and z , which, when discretized with the conventional symmetric finite difference operator, turns out to be impossible. Hence we set out to define a novel finite difference operator, which besides the summation-by-parts property also remains neutral under reparametrization. We proceed in Section 2 by formulating an explicit expression for Eq. (7) for the dissipative dynamics of a heavy quarkonium particle, interacting with a hot environment. We will discuss the preservation of the defining properties of the density matrix in the continuum and pinpoint where it fails after discretization. Using this insight we will in Section 3 construct a novel reparametrization neutral SBP operator, which, as we show, retains the continuum properties in the discretized evolution equations. In Section 4 we will present numerical results from the simulation of quarkonium dissipative dynamics, showcasing the successful preservation of the continuum properties of the density matrix. We close with a brief conclusion and outlook in Section 5

Quarkonium Lindblad master equation
As a concrete example of a Lindblad equation describing dissipative dynamics of a phenomenologically relevant system, we present the case of heavy quark-anti-quark bound states, so called heavy quarkonium at high temperature. The dissipative dynamics of these bound states immersed in a thermal medium play an important role in our understanding of heavy-ion collisions carried out e.g. at the Large Hadron Colldier at the CERN Laboratory. In such collisions, nuclei of heavy atoms are smashed into each other at ultra-relativistic momenta, so that the protons and neutrons making up the nuclei become compressed and heated to temperatures beyond 200, 000× the temperatures present in the core of the sun. In turn they melt into their microscopic constituents, the light quarks and gluons, which form a so called quark-gluon-plasma (QGP). Some of the kinetic energy in the original projectiles is converted into additional particles, such as heavy quark-anti-quark pairs, which may form quarkonium bound states that find themselves immersed in the approximately locally thermal QGP. Using the quantum field theory quantum-chromo-dynamics (QCD) one wishes to understand how the interaction between the quarkonium and the medium affects their binding properties. If we translate these in-medium modifications into changes in the quarkonium particle production rate, we may use their measured yields in heavy-ion collisions to deduce the properties of the QGP created therein.
In [25,26] the Lindblad equation for a single quarkonium particle at high temperature was derived. Here the term high refers to the fact that the strong interactions between quarks and gluons become weaker at high energy scales and thus a weak-coupling expansion could be deployed. The explicit expressions for the subsystem Hamiltonian and the Lindblad operators in terms of the three-dimensional coordinates of the quark and anti-quark x Q and xQ and their momenta p Q ,pQ read The Lindblad operators in this case depend on a continuous momentum variable k and a discrete "color" index a related to the triple valued charge of strongly interacting particles. I.e. each "entry" of the Hamiltonian and the Lindblad operators corresponds to a 6 × 6 matrix built from direct products of the 3 × 3 Gell-Mann matrices denoted by t a .
The form ofL k,a is intuitively understandable: there are two terms, one for the interaction of the medium with the quark-, another term for the interaction with the anti-quark constituent making up the quarkonium particle. The relative minus sign between the two terms in the second line of Eq. (12) indicates that we are dealing with a system consisting of a particle and anti-particle. All aspects of the dynamics are captured in two quantities, the real-valued potential V (x) and the dissipation kernel D(x). These two quantities are intimately related to the real-and imaginary part of the proper real-time heavy-quark potential computed perturbatively in [27,28,29] and in numerical simulations of the strong interactions (so called lattice QCD) in [30,31,32,33]. It has been shown that the effects of dissipation are encoded in the momentum k dependent terms in the parentheses, without which the dissipationless evolution of the recoilless-limit is obtained.
The description of the two-body quarkonium system may be simplified by going over to relative and center of mass coordinates for the quark-anti-quark degrees of freedom Tracing out also the center of mass coordinate, it was shown in [34] that the following relative coordinate operators ensuẽ L rel k,a = D (k) In the following we will further simplify the description by considering the quarkonium particle to be at rest P = 0 and neglecting the explicit matrix structure of the Lindblad operators and the Hamiltonian, essentially setting t a to unity and retaining only a single color component. Evaluating the ensuing Lindblad equation in the coordinate space basis for the relative coordinates, we obtain the following partial differential equation with spatially and potentially temporally varying coefficient terms The vectorial derivatives acting on the coordinate vectors x and y are defined as ∇ x = ( ∂ ∂x 1 , ∂ ∂x 2 , ∂ ∂x 3 ) and ∇ y = ( ∂ ∂y 1 , ∂ ∂y 2 , ∂ ∂y 3 ) respectively and their components denoted by ∇ i x and ∇ i y . The scalar function F 1 , the vectorial F 2 and the tensorial F ij 3 have been introduced to conveniently summarize the contributions arising from the dissipation kernel and the function A(x) = D(x)/8T 2 . Up to this point no approximation beyond the weak coupling expansion and time coarse graining enter our description. While our long-term goal is to solve the full three-dimensional dynamics of quarkonium based on Eq. (16), we restrict ourselves in this paper to the one-dimensional case. It already presents us with many of the relevant technical challenges, while requiring significantly less computational resources for its implementation. It furthermore simplifies the presentation, without compromising with the fundamental computational development.
In one dimension, derivatives in the x and y coordinate do not carry indices anymore and Eq. (16) reduces to Here the physics of the dissipation kernel D(x) enters via three real-valued scalar functions, which, keeping in mind that A(x) = D(x)/8T 2 , read

Conservation of defining properties in one-dimension
In order for positivity and hermiticity of the density matrix to be conserved in the Lindblad formalism, the function D(x) in momentum space needs to be positive and real. As the dissipation kernel is supplied as external input, this property can be explicitly checked for and we will make sure it is fulfilled in the simulations that follow. Here we focus on the preservation of the trace in Eq. (20), which presents the central challenge in its discretization. In the functional language, the trace over quantum states translates into an integration of ρ rel (x, y, t) over x and y in the presence of a delta function δ(x − y).
Some of the terms in the trace over the first, second, fourth and fifth line of Eq. (20) vanish identically, due to the arguments of the V and F functions being evaluated at x = y, e.g.
In addition, the following term from lines four and five vanishes This can be seen by inspecting the properties of the function D(x) and the definition of F 2 , in which the first and third derivative of D(x) enters (remember A(x) = D(x)/8T 2 ). The derivation of the Lindblad equation from QCD in [25,26] leads to a function D(x) that possesses a maximum around the origin and which is furthermore symmetric around the origin. Thus, when we take the first and third derivative of D(x) at the origin, both contributions vanish, i.e. F 2 (0) = 0. (We ensure that this property is respected in our numerical simulations, see Eq. (56).) For some terms in the trace over Eq. (20) we need to apply integration by parts. Take e.g. the derivatives in the first line In order to show that these two contributions cancel each other, we need to be able to transform the double derivative in x into a corresponding derivative in y, which is possible due to the delta-function. We get from twice integrating by parts and the fact that we consider periodic boundary conditions. Now we exploit the anti symmetry of the argument of the delta function to replace ∂ ∂x by − ∂ ∂y twice. (One may use a regularized form of the delta-function here to make the operation mathematically well defined.) Applying integration by parts twice again, we thus obtain which hence cancels the term in Eq. (25). Similarly, by application of integration by parts in the trace over terms in line six and seven of Eq. (20), we find that the following combination vanishes The need to perform integration by parts in the continuum theory tells us that the difference operators, to be deployed for discretization of Eq. (20), need to fulfil the summation-by-parts property.
For the other terms to vanish, additional continuum properties of derivatives are required. Let us start with the remaining F 3 terms in the last two lines of Eq. (20) This expression will not vanish by itself but instead produces a remnant term, which in turn will cancel with the A dependent contributions in the trace over Eq. (20). At first sight moving around the derivatives on the first F 3 term (denoted as T 61 ) appears to involve the application of the product rule which would lead to subsequent numerical complications, as splitting would be required [35]. Note however the particular form of this coefficient. It depends only on the sum z = x + y of the coordinates, while the delta function only depends on z = x − y, which invites us to introduce Let us have a look at how the reparametrization property of the differentials of the two sets of coordinates in the continuum can be used to rewrite the mixed derivative term in Eq. (27). We obtain The result of Eq. (30) can be directly applied to the mixed derivative term T 61 in Eq. (27): In Eq. (32), we carried out one integration by parts in both x and y. Since δ(x − y) only depends on z and not z , the derivative ∂/∂z after integration by parts acts solely on the F 3 term.
Let us suggestively rewrite the derivatives over z in terms of a variable ξ, which, due to the presence of the delta function, we will later identify with x or y Note that the factor 2 from Eq. (32) has disappeared due to the application of the chain rule. Using the definition in Eq. (19), we can turn these two F 3 terms into derivatives acting on the A function Comparing the above expression to the fourth and fifth line of Eq. (20) we see that it has a form very similar to the A terms present there Making use of the fact that in Eq. (35), inside the trace, we can replace ∂ 3 ∂ξ 3 by either ∂ 3 ∂x 3 or ∂ 3 ∂y 3 , we can cancel parts of the terms in T 7 . Note that due to the factor 1/2 present in Eq. (34) we obtain a rest term when summing T 6 and T 7 .
In order to bring this rest term into the form necessary to cancel the last remaining A terms in Eq. (20) it may be conveniently expressed in the ξ derivatives of Eq. (35) Let us now summarize the two terms as a single expression with derivative in z We may now carry out integration by parts in z (by shifting ∂ ∂z to 1 2 ( ∂ ∂x + ∂ ∂y )), exploiting that δ(x − y) does not depend on z . Using again the fact that inside the trace we can replace ∂ 3 ∂ξ 3 by either ∂ 3 ∂x 3 or ∂ 3 ∂y 3 we arrive at the final expression . which cancels identically with the only remaining A terms in the third line of Eq. (20). This concludes the explicit demonstration that the one-dimensional Lindblad master equation preserves the unit trace property of the reduced density matrix.
Remark. The lesson learned for the discretization of Eq. (20) is that the application of the product rule is not necessary in order to conserve the trace, as long as the reparametrization property of Eq. (28) is fulfilled. This bodes well, as it is well known that discretizations of the difference operator in general violate the product rule, even if they fulfill e.g. the summation by parts property [35]. To summarize, we need difference operators that can integrate by parts in x, y and differentiate in z, z .
However, the standard difference operator is not neutral under reparametrization, as is evidenced by (assuming the same equidistant discretization ∆ in x and y) Our goal thus is to construct a SBP operator that remains neutral under the reparametrization (x, y) → (z, z ).

A reparametrization-neutral summation by parts (RN-SBP) operator
We proceed to construct a novel SBP difference operator, which strictly implements the reparametrization property Eq. (28), restricting ourselves here to the case of an equidistantly discretized function ρ(x, y) with ∆x = ∆y = ∆ and periodic boundary conditions. After choosing an explicit ordering of the discretized density matrix ρ on the now two-dimensional (x,y) grid we introduce the shift operators S + and S − as follows Note that the consecutive application of the two shift operators yields S + S − = S − S + = 1, where 1 is the identity. Our strategy is to combine these shifts together with regular SBP difference operators to achieve the desired reparametrization neutrality. Indeed, I x ⊗ S + shifts rows upward on the discretized grid, while S + ⊗ I y shifts columns to the right, with the inverse operations naturally being I x ⊗ S − and S − ⊗ I y respectively. For periodic boundary conditions the simplest second order periodic SBP operator is constructed using the integration prescription H = ∆1 and where D ≡ H −1 Q. Since Q + Q = 0, we get and the SBP property simplifies to (u, Dv) H = −(Du, v) H .
As indicated by Eq. (40), for the reparametrization property to hold we need to express the x-and y-derivative of the function ρ i,j in terms of its values at the neighboring diagonal corners, i.e. in the variables x + y and x − y. This is possible if we compute the average of the naive finite differences once shifted up and down in rows and columns respectively. With this strategy we arrive at the definition of the following reparametrization neutral summation-by-parts operators (RN-SBP) The explicit expressions applying D x and D y to the function ρ at (x i , y j ) are This construction also preserves the summation by parts property, as we can write with H = H ⊗ H, and It follows that Hence, if u, v ∈ R N 2 , then (u, D x v) H = −(D x u, v) H (using the same argument as in (42)). Similarly we get (u, Let us next define the corresponding RN-SBP operators in z and z as in Eq. (29), which by (44), (45) on the index level becomes Naturally these new operators show the following behavior: if we have two functions u, v ∈ R N 2 and v depends only on z = (x − y), i.e. v is constant along all lines y = x + c then where • denotes the Hadamard (elementwise) product. Similarly if u depends only on z = (x + y) we have Let us now show that the novel RN-SBP operators presented above allow us to implement in a discrete fashion, all manipulations deployed in the proof of the preservation of unit trace in Section 2.
The first set of manipulation we need to realize discretely is related to the terms T 4 , T 5 and T 6 . There, one must first perform summation by parts, which our novel RN-SBP operator implements by construction. Next we need to change x into y derivatives and vice versa via the delta function. Writing explicitly we get Another type of operation is required to treat the T 6 term. In the continuum it amounts to the steps in Eq. (30) which are identical for the RN-SBP operator Note that already the first line of the above operations would not hold if implemented with the naive finite difference operator. The treatment of the T 6 term requires our operator to fulfill an additional property. When the z derivative after summation by parts, acts on δ(x − y)F 3 ( x+y 2 ) we need it to affect only F 3 . Inspecting Eq. (51) we conclude that the RN-SBP operator either the ground or the first excited state, corresponding to the lowest lying φ 0 (x) or next to lowest eigenvector φ 1 (x) of the Hamiltonian where M denotes the heavy quark mass and V (x) a real-valued interaction potential. At each step in the time evolution, we may express the density matrix in the basis of these eigenvectors ρ mn (t) = dx dyφ * m (x)ρ(x, y, t)φ n (y). The survival probability is then read off from the diagonal entries e.g. P 0 (t) = ρ 00 (t).
Eq. (20) has been solved approximately in a previous study [34] using a stochastic unravelling of the master equation in terms of an ensemble of wave functions evolving under a non-linear stochastic Schrödinger equation. As this so called quantum state diffusion approach necessitated additional approximations in the heavy quark velocity, the present study will provide an important crosscheck of the validity of that computation. In [38] the master equation on the other hand has been solved through stochastic unraveling in the recoilless limit, which in the language of Eq. (20) amounts to neglecting all contributions coming from the terms except for the only D(x) in F 1 (x).
Let us remark that standard methods of operator splitting, highly efficient for the solution of the regular Schrödinger equation, fail for Eq. (20), as we are faced with a partial differential equation including variable coefficients. The spatial dependence of the different F terms leads to significant contributions of commutators in the Trotter decomposition, which break the naive counting of the Strang scheme [39]. In order to implement the left hand side of the master equation efficiently we therefore deploy the PETSC [40,41] sparse and distributed matrix library.
Similarly to the values chosen in [34] we discretize the density matrix on a grid of N = 256 points in x and y direction each with periodic boundary conditions. Using the mass M of the heavy quark to express all dimension full quantities, we have for the spatial spacing ∆ = 1/M and a time step of ∆t = 0.1M (∆) 2 . In the simulations we will explicitly set M = 1, from which follows ∆t = 0.1∆. (We have checked that reducing the time step ∆t further does not significantly change the outcome of our simulations.) The interactions among the quark-anti-quark pair are captured in a model potential and dissipation kernel

Re[λ]
exact hermiticity and positivity Initial t=0 tM=100 tM=300 tM=600 tM=900 tM=1200 tM=6000 γ = T /π and as regularization for the Coulomb potential x r = 1/M . The simulation will be performed for three different temperatures, T = 0.05M , T = 0.1M and T = 0.3M . Note that since the the density matrix ρ(x, y) needs to fulfill periodic boundary conditions in each variable x, y independently we are lead to an additional constraint for the functions D(x). I.e. if ρ(x + L, y) = ρ(x, y + L) = ρ(x, y) then D(x + L/2) = D(x), which tells us that the function D must be periodic over half the box size.
In order to carry out the Crank-Nicolson step for Eq. (20) we have to solve a linear system of equations. Exploiting the sparse nature of the update matrix we choose to utilize the distributed sparse matrix format provided by the PETSC library. Contrary to dense matrix algorithms here the solution is found iteratively and thus approximately via the GMRES algorithm. As a compromise between precision and computational speed we select a solution tolerance of ∆ GMRES = 10 −14 and a maximum number of steps in the iterations N GMRES = 100. The error introduced by the approximation to the true solution was the dominant source of error in our simulations.
Let us start by investigating the defining properties of the density matrix for the representative example of T = 0.1M . In Fig. 1 we plot the real-and imaginary part of the eigenvalues of the discretized density matrix at different times during the evolution (colored data points). The values were obtained using the SLEPC distributed eigensolver library [42]. We have chosen the time window such, that at the latest time tM = 6000 the system is very close to a stationary state, which leaves the survival probabilities of the two lowest lying states of interest unchanged (see also Fig. 4). As expected from a density matrix initialized using a single normalized eigenstate of the system Hamiltonian, it contains a single non-vanishing eigenvalue of value unity while the rest of the eigenvalues vanishes. During the time evolution, the real-part of the eigenvalues remains positive, which is a manifestation of the conservation of positivity of ρ. On the other hand we clearly see deviations of the imaginary part of the eigenvalues into the complex plane, which amounts to a violation of the exact hermiticity (denoted by the gray line) of the density matrix. The deviations for our choice of rather large time steps of ∆t = 0.1/M however remain at the permille level at all times. We have checked that replacing the RN-SBP derivative operator with its naive counterpart leaves the positivity and hermiticity properties of the time evolution virtually unchanged. We continue with an inspection of the trace of the density matrix, based on the novel RN-SBP derivative operator and plotted as solid lines in Fig. 2 for T = 0.1M and T = 0.3M . Due to the properties of the RN-SBP operator derived in Section 3, we find that the trace values are excellently preserved with a maximum deviation from unity of 10 −9 . For comparison purposes we also plot the values of the trace as obtained with the naive difference operator as dashed lines. One can clearly see that the violation of the unit trace property grows with time and that the strength of the deviation depends on the particular parameters of the dynamical evolution. We have checked that the minute deviation from unit trace in case of the RN-SBP operator reduces when we lower the tolerance ∆ GMRES for the iterative solution of the the Crank-Nicolson step. Reducing the tolerance further, at some point the errors introduced by the finite ∆ GMRES are no longer the dominant source of error and instead it is the finiteness of the time step ∆t. At this point both ∆t and ∆ GMRES need to be reduced in tandem for the results to further improve.
One may question whether an apparently small deviation from unit trace by less than a percent, as is visible in Fig. 2, has any significant consequences for the physics outcome of our simulation. To this end we compare the We conclude, based on the above comparison, that that the combination of the Crank-Nicolson scheme with our novel RN-SBP operator provides an accurate discrete representation of the dynamics described by Eq. (20).
Having convinced us of the inner workings of the underlying discretization, we proceed to investigate the physics results of our simulation. In Fig. 4 we plot the survival probabilities of the ground and first excited states of the system Hamiltonian P 0 and P 1 for the case of T = 0.05M as gray dashed lines, for T = 0.1M as colored solid lines, and for T = 0.3M as colored dashed lines. As crosscheck of our previous work we also plot the values obtained from an approximate stochastic unravelling of the master equation via the quantum-state-diffusion approach as open circles.
We find that the stochastic unraveling provides a very good description of the dynamics of the ground state up to tM = 3000, which is already longer than what would be needed in the simulation of a realistic heavy-ion collision. At late times deviations from the direct solution of the master equation become visible, which however remain at the 20% level. The first excited state shows similar deviations, which set in at a bit earlier times around tM = 1250.
The simulation of the evolution of the quarkonium system at different temperatures proceeds in a consistent fashion. In a colder environment, such as T = 0.05M (gray dashed lines) the medium is unable to interfere with the quarkonium binding as efficiently as is the case at T = 0.1M . This on the one hand leads to a slower decay of the ground state survival probability and on the other hand produces a less rapid population of the excited states. Conversely in a hotter environment, such as at T = 0.3M (colored dashed lines) the ground state is more efficiently depleted while rapidly populating the excite states.
The relative abundances between the states in thermal equilibrium are expected to be governed by the Boltzmann distribution, which means that the two curves will lie further apart at T = 0.05M and closer together at T = 0.3M than at T = 0.1M . As we see in Fig. 4, at first rising temperatures lead to stronger occupation of the exited states but eventually at high enough temperatures also their contribution will become suppressed.
We see that a steady state is reached at different times for different temperatures. Comparing T = 0.05M and T = 0.1M we see that at higher temperature the steady state emerges already at tM = 3000, while at the lower temperature we need to wait until around tM = 30000. Interestingly for T = 0.3M we find that the relative abundances between the ground and first excited state are established quite quickly, around M t = 1500 but that the overall amplitude of the survival decreases over time, indicating that the excited states are not yet equilibrated 1 .
The properties of the steady state reached at late times may be investigated by inspection of the abundances of the individual states present in the system. In Fig. 5 Figure 5: Comparison of the abundances of the individual states in form of the survival probabilities P i versus the energy of these states (data points). As solid lines we plot exponential fits, motivated by the expected Boltzmann distribution, with the corresponding best fit temperature given in the key.
by the expectation that eventually a Boltzmann distribution will emerge, we also plot exponential fits to the data and provide the best fit value of the "inverse slope parameter", the temperature, in the key. For all systems at T = 0.05M , T = 0.1M and T = 0.3M we find that the fit captures all of the ten lowest lying states very well and a temperature emerges, which, while not exactly at the environment, lies very close to it. Note that as was shown in [34] such a deviation is not unexpected, as thermalization with the same temperature, can only be proven in the classical limit and for small velocities, quantum corrections may lead to small deviations.
In previous studies, such as Refs. [43,44,38], the quarkonium dynamics were investigated in the recoilless limit, which allows Eq. (20) to be unraveled in terms of a stochastic potential. It amounts to neglecting all F 2 and F ij 3 terms and retaining only the D(x) term in F 1 . For illustrative purposes let us compare the full dissipative dynamics to this approximation in the left panel of Fig. 6, similar to a comparison previously carried out in [34]. As dissipative effects are absent, the fluctuations of the environment transfer energy into the quarkonium system which is unable to release it back to the medium. Hence the system heats up unabated and one expects that eventually all states will become destabilized. And indeed the dashed lines clearly show this behavior, as the survival probabilities P 0 and P 1 eventually fall on top of each other. We also emphasize that already at early times a significant deviation from the full dissipative dynamics is observed for the ground state. As a crosscheck we show in the right panel of Fig. 6 the comparison of the dynamics in the recoilless limit obtained from the solution of the approximate master equation (solid line) and via the stochastic potential approach (data points) using here as an exception the parameters of [38], confirming the correctness of the numerics of that study.
The failure of the recoilless limit to thermalize also manifests itself clearly in the total energy E = Tr[Ĥρ] of the system, as shown in Fig. 7. The solid lines represent the fully dissipative dynamics, which asymptote against a constant value at late times. On the other hand the dissipationless dynamics depicted via dashed lines for T = 0.1M and T = 0.05M exhibit an unabated rise that is linear in time at late times. At T = 0.3M the energy in the dissipationless limit eventually runs into a constant too, which is not related to thermalization, but to reaching the "infinite temperature limit" on a lattice with finite extent and lattice spacing.

Conclusion and Outlook
In this study, we have presented an improved numerical open quantum systems treatment of heavy quarkonium at high temperature, via its Lindblad equation. We showed that in order to fulfill the defining properties of the quarkonium density matrix ρ(x, y), guaranteed by the Lindblad equation in the continuum, we need to be able to carry out integration by parts and exploit the reparametrization of the system between the (x, y) and (z = x − y, z = x + y) coordinates. In order for the properties to hold also after discretization we thus developed a novel RN-SBP derivative operator, which not only fulfills the summation-by-parts property, but also remains neutral under the above mentioned reparametrization.
With the novel RN-SBP operator at hand, we presented a numerical implementation of Eq. (20) using the Crank-Nicolson approach. It allowed us to evolve the density matrix in time, while preserving its positivity, hermiticity and trace accurately. In turn we were able to not only crosscheck the validity of previous computations, based on the approximate Quantum State Diffusion approach, but also obtain a more robust result for the steady state occupancies at late times.
There are several directions to explore next: on the one hand it will be interesting to formulate the novel RN-SBP operator for use in higher dimensions, as the ultimate goal is to solve Lindblad equations, such as Eq. (20) for three-dimensional coordinate x and y. Given its simple structure in one dimension, its generalization via shifts along different axes appears straight forward, but has to be verified explicitly. In addition one may want to consider how to formulate compatible second order derivatives, which require less off-diagonal terms, compared to the naive application of the first order derivative operator twice. In order for RN-SBP operators to play out their full strength in precision studies of dissipative dynamics, also higher order incarnations of the first derivative operator need to be formulated.
In summary the novel RN-SBP operator presented here provides an interesting discrete implementation of a continuum derivative property not treated explicitly in the literature so far. In turn, we hope that it will be of benefit in many other numerical settings, be it for the study of dissipative dynamics or beyond.