We provide the explicit form of our construction and explain how exact detailed balance can be achieved through the key algorithmic subroutines and circuits.<\/p>\n
Constructing the full Lindbladian<\/p>\n
Recall that our main result considers the following Lindbladian in the Schr\u00f6dinger picture:<\/p>\n
$$\\begin{array}{c}{\\mathcal{L}}[\\cdot ]\\,:= \\,\\sum _{a\\in A}\\,\\mathop{\\underbrace{-i[{C}^{a},\\cdot ]}}\\limits_{ \\mbox{“} {\\rm{c}}{\\rm{o}}{\\rm{h}}{\\rm{e}}{\\rm{r}}{\\rm{e}}{\\rm{n}}{\\rm{t}}\\mbox{”}}\\,\\,+\\mathop{\\overbrace{{\\int }_{-{\\rm{\\infty }}}^{{\\rm{\\infty }}}\\gamma (\\omega )\\left(\\mathop{\\underbrace{{\\hat{A}}^{a}(\\omega )[\\cdot ]{\\hat{A}}^{a}{(\\omega )}^{\\dagger }}}\\limits_{ \\mbox{“} {\\rm{t}}{\\rm{r}}{\\rm{a}}{\\rm{n}}{\\rm{s}}{\\rm{i}}{\\rm{t}}{\\rm{i}}{\\rm{o}}{\\rm{n}}\\mbox{”}}-\\mathop{\\underbrace{\\frac{1}{2}\\{{\\hat{A}}^{a}{(\\omega )}^{\\dagger }{\\hat{A}}^{a}(\\omega ),\\cdot \\}}}\\limits_{ \\mbox{“} {\\rm{d}}{\\rm{e}}{\\rm{c}}{\\rm{a}}{\\rm{y}}\\mbox{”}}\\right){\\rm{d}}\\omega }}\\limits^{{\\rm{d}}{\\rm{i}}{\\rm{s}}{\\rm{s}}{\\rm{i}}{\\rm{p}}{\\rm{a}}{\\rm{t}}{\\rm{i}}{\\rm{v}}{\\rm{e}}\\,{\\rm{p}}{\\rm{a}}{\\rm{r}}{\\rm{t}}}\\,.\\end{array}$$<\/p>\n
\n (9)\n <\/p>\n
Roughly, it resembles the Davies generator (equation (5<\/a>)) but is carefully modified to maintain quantum detailed balance while ensuring locality and efficiency. We begin by reviewing the detailed balance condition for the Davies generator to illustrate the key ingredients and differences in our construction. We may regroup the above according to the jumps \\({\\mathcal{L}}={\\sum }_{a\\in A}{{\\mathcal{L}}}^{a}\\), and study each term individually, so in what follows we drop the label a for simplicity, that is, substitute Aa\u2009\u2190\u2009A.<\/p>\n Detailed balance of the Davies generator<\/p>\n For a Hermitian jump A, recall that the Davies generator \\({{\\mathcal{L}}}_{Davies}[\\cdot ]\\,:= \\,\\sum _{\\nu }\\gamma (\\nu ){A}_{\\nu }[\\cdot ]{A}_{\\nu }^{\\dagger }-\\) \\(\\frac{\\gamma (\\nu )}{2}\\{{A}_{\\nu }^{\\dagger }{A}_{\\nu },\\cdot \\}\\) satisfies the KMS-detailed balance condition (equation (4<\/a>)), that is, satisfies the superoperator equation<\/p>\n $${\\mathcal{L}}[\\cdot ]={\\rho }^{\\frac{1}{2}}{{\\mathcal{L}}}^{\\dagger }[\\,{\\rho }^{-\\frac{1}{2}}\\cdot {\\rho }^{-\\frac{1}{2}}]{\\rho }^{\\frac{1}{2}},$$<\/p>\n \n (10)\n <\/p>\n where \\({{\\mathcal{L}}}^{\\dagger }\\) is the superoperator adjoint of \\({\\mathcal{L}}\\) defined by \\({\\rm{Tr}}[{({\\mathcal{L}}[X])}^{\\dagger }Y]={\\rm{Tr}}[{X}^{\\dagger }{{\\mathcal{L}}}^{\\dagger }[Y]]\\) for all X,\u00a0Y. Here, detailed balance hinges on the following exact operator-valued symmetries: for all Bohr frequency \u03bd \u2208 B(H),\u00a0<\/p>\n $${\\rho }^{-\\frac{1}{2}}{A}_{\\nu }{\\rho }^{\\frac{1}{2}}={{\\rm{e}}}^{\\frac{\\beta \\nu }{2}}{A}_{\\nu }\\,(\\text{conjugation identity}),$$<\/p>\n $${A}_{-\\nu }={({A}_{\\nu })}^{\\dagger }\\,(\\text{adjoint property}).$$<\/p>\n The conjugation identity is rooted in that A\u03bd are eigenoperators [H,\u2009A\u03bd]\u00a0=\u00a0\u03bdA\u03bd of the commutator [H, \u22c5 ], highlighting a special role played by the Bohr-frequency decomposition A\u03bd. The adjoint property says that for a Hermitian jump A, the transition amplitudes associated with energy difference \u03bd are paired with the reverse difference\u00a0\u2013\u03bd,\u00a0reminiscent of a Fourier transform symmetry of real functions.<\/p>\n As a consequence, the decay part readily satisfies detailed balance by itself because of the adjoint property, because the operator \\({({A}_{\\nu })}^{\\dagger }{A}_{\\nu }={A}_{-\\nu }{A}_{\\nu }\\) in the decay part preserves the energies and commutes with the Hamiltonian (and hence with \u03c1). For the transition part,<\/p>\n $$\\begin{array}{c}\\sum _{\\nu }\\gamma (\\nu ){\\rho }^{\\frac{1}{2}}{({A}_{\\nu })}^{\\dagger }{\\rho }^{-\\frac{1}{2}}\\cdot {\\rho }^{-\\frac{1}{2}}{A}_{\\nu }{\\rho }^{\\frac{1}{2}}\\\\ \\,\\,=\\,\\sum _{\\nu }\\gamma (\\nu ){e}^{\\beta \\nu }{({A}_{\\nu })}^{\\dagger }\\cdot {A}_{\\nu }\\,\\,\\,\\,\\,\\,\\text{(by the conjugation identity)}\\,\\\\ \\,\\,=\\,\\sum _{\\nu }\\gamma (-\\nu ){A}_{-\\nu }\\cdot {({A}_{-\\nu })}^{\\dagger }\\,=\\,\\sum _{\\nu }\\gamma (\\nu ){A}_{\\nu }\\cdot {({A}_{\\nu })}^{\\dagger }.\\end{array}$$<\/p>\n \n (11)\n <\/p>\n The second equality uses the Kubo\u2013Martin\u2013Schwinger condition (in the frequency domain) for the transition weights<\/p>\n $$\\gamma (\\nu ){{\\rm{e}}}^{\\beta \\nu }=\\gamma (-\\nu ).$$<\/p>\n \n (12)\n <\/p>\n The main obstacle towards a Lindbladian with exact detailed balance and efficient implementation is the lack of algorithmic access to A\u03bd in the presence of a dense spectrum, as approximations to A\u03bd can easily break the exact symmetry conditions in detailed balance18<\/a>,20<\/a>,51<\/a>.<\/p>\n Operator Fourier transforms<\/p>\n The key to both efficiency and exact detailed balance is a careful relaxation of A\u03bd that preserves some aspects of the symmetries using the operator Fourier transform damped by a Gaussian filter\u00a0\\(f(t)\\,:= \\,{{\\rm{e}}}^{-{\\sigma }^{2}{t}^{2}}\\sqrt{\\sigma \\sqrt{2\/{\\rm{\\pi }}}}\\) (normalized by \\({\\int }_{-\\infty }^{\\infty }{| f(t)| }^{2}{\\rm{d}}t=1\\)) with a tunable width \u221d 1\/\u03c3 (setting \\(\\sigma =\\frac{1}{\\beta }\\) recovers (equation (6<\/a>)):<\/p>\n $$\\begin{array}{c}\\widehat{A}(\\omega )\\,:= \\,\\frac{1}{\\sqrt{2{\\rm{\\pi }}}}{\\int }_{-{\\rm{\\infty }}}^{{\\rm{\\infty }}}{{\\rm{e}}}^{{\\rm{i}}Ht}A{{\\rm{e}}}^{-{\\rm{i}}Ht}{{\\rm{e}}}^{-{\\rm{i}}\\omega t}f(t){\\rm{d}}t\\\\ \\,\\,\\,=\\,\\frac{1}{\\sqrt{{\\sigma }\\sqrt{2{\\rm{\\pi }}}}}\\sum _{\\nu \\in B(H)}{A}_{\\nu }\\,{{\\rm{e}}}^{-\\frac{{(\\omega -\\nu )}^{2}}{4{{\\sigma }}^{2}}}.\\end{array}$$<\/p>\n \n (13)\n <\/p>\n Like A\u03bd in the Davies generator (equation (5<\/a>)), our operator Fourier transforms \\(\\widehat{A}(\\omega )\\) are labelled by energy differences, but here the parameter \u03c9\u2009\u2208\u2009[\u2212\u221e,\u2009\u221e] takes continuous values, without referring to the true Bohr frequencies. When the uncertainty vanishes \u03c3\u2009\u2192\u20090, we recover \\(\\mathop{\\mathrm{lim}}\\limits_{\\sigma \\to 0}\\sqrt{\\sigma \\sqrt{2{\\rm{\\pi }}}}\\widehat{A}(\\omega )=\\sum _{\\nu \\in B(H)}{\\mathbb{1}}(\\nu =\\omega ){A}_{\\nu }\\); when the energy uncertainty is finite \u03c3\u2009\u2260\u20090, the operator Fourier transforms still maintain exact algebraic properties crucial for detailed balance. First, conjugating with the Gibbs state preserves the form of operator Fourier transforms, albeit causing a shift and rescaling<\/p>\n $$\\begin{array}{c}{\\rho }^{-\\frac{1}{2}}\\widehat{A}(\\omega ){\\rho }^{\\frac{1}{2}}\\,=\\,\\frac{1}{\\sqrt{{\\sigma }\\sqrt{2{\\rm{\\pi }}}}}\\sum _{\\nu \\in B(H)}{A}_{\\nu }{{\\rm{e}}}^{\\frac{\\beta \\nu }{2}}\\exp \\left(-\\frac{{(\\omega -\\nu )}^{2}}{4{{\\sigma }}^{2}}\\right)\\,\\,\\,\\,({\\rm{b}}{\\rm{y}}\\,(13)\\,\\,\\text{and the conjugation identity)}\\,\\\\ \\,\\,\\,\\,=\\,\\frac{1}{\\sqrt{{\\sigma }\\sqrt{2{\\rm{\\pi }}}}}\\sum _{\\nu \\in B(H)}{{\\rm{e}}}^{\\frac{\\beta \\omega }{2}+\\frac{{{\\sigma }}^{2}{\\beta }^{2}}{4}}{A}_{\\nu }\\exp \\left(-\\frac{{(\\omega -\\nu +{{\\sigma }}^{2}\\beta )}^{2}}{4{{\\sigma }}^{2}}\\right)\\,\\,\\,\\,\\,\\text{(by completing the square)}\\,\\\\ \\,\\,\\,\\,=\\,{{\\rm{e}}}^{\\frac{\\beta \\omega }{2}+\\frac{{{\\sigma }}^{2}{\\beta }^{2}}{4}}\\widehat{A}(\\omega +{{\\sigma }}^{2}\\beta ),\\,\\,\\,\\,\\,\\text{(by shift-rescale symmetry)}\\end{array}$$<\/p>\n which reflects the fact that multiplying a Gaussian distribution by an exponential weight must shift the mean \\({{\\rm{e}}}^{{(x-a)}^{2}-2bx}={{\\rm{e}}}^{{(x-a-b)}^{2}-2ab-{b}^{2}}\\). Second, even though the operator Fourier transform is a linear combination of different Bohr frequencies, the adjoint property holds exactly when A\u2009=\u2009A\u2020 as<\/p>\n $$\\begin{array}{l}\\widehat{A}(\\,-\\,\\omega )\\,=\\,\\frac{1}{\\sqrt{2{\\rm{\\pi }}}}{\\int }_{-{\\rm{\\infty }}}^{{\\rm{\\infty }}}{{\\rm{e}}}^{{\\rm{i}}Ht}A{{\\rm{e}}}^{-{\\rm{i}}Ht}{{\\rm{e}}}^{{\\rm{i}}\\omega t}f(t){\\rm{d}}t\\\\ \\,\\,\\,\\,=\\,\\frac{1}{\\sqrt{2{\\rm{\\pi }}}}{\\int }_{-{\\rm{\\infty }}}^{{\\rm{\\infty }}}{{\\rm{e}}}^{{\\rm{i}}Ht}{A}^{\\dagger }{{\\rm{e}}}^{-{\\rm{i}}Ht}{({{\\rm{e}}}^{-{\\rm{i}}\\omega t}f(t))}^{\\ast }{\\rm{d}}t=\\widehat{A}{(\\omega )}^{\\dagger }\\,\\,\\,\\,({\\rm{s}}{\\rm{i}}{\\rm{n}}{\\rm{c}}{\\rm{e}}\\,\\,{f}^{\\ast }(t)=\\,f(t)).\\end{array}$$<\/p>\n The above two exact symmetries appear to be absent in previous approaches that attempted to directly measure energy differences. Now we show that quantum detailed balance is a consequence of these exact symmetries and related algebraic properties of the Gaussian uncertainty in operator Fourier transforms, which hold exactly despite not measuring energies to high precision.<\/p>\n Exact detailed balance with finite uncertainty<\/p>\n Although the operator Fourier transform does not yield an exact representation of A\u03bd, the uncertainty has a specific structure with distinctive symmetries\u2014arising from the interplay between the Gaussian weighing and the exponential form of the Boltzmann factors\u2014and consequently can be exactly compensated by an appropriate shift in the transition weights. In the following, we prove that the transition part (equation (9<\/a>)) satisfies detailed balance (equation (10<\/a>)): \\({\\mathcal{T}}\\,[\\cdot ]={\\rho }^{\\frac{1}{2}}{{\\mathcal{T}}}^{\\dagger }[\\,{\\rho }^{-\\frac{1}{2}}\\cdot {\\rho }^{-\\frac{1}{2}}]{\\rho }^{\\frac{1}{2}}\\). The above shift-rescale and adjoint symmetries yield<\/p>\n $$\\begin{array}{c}{\\rho }^{\\frac{1}{2}}{{\\mathcal{T}}}^{\\dagger }[{\\rho }^{-\\frac{1}{2}}\\cdot {\\rho }^{-\\frac{1}{2}}]{\\rho }^{\\frac{1}{2}}\\\\ \\,\\,=\\,{\\int }_{-{\\rm{\\infty }}}^{{\\rm{\\infty }}}\\gamma (\\omega ){\\left({\\rho }^{-\\frac{1}{2}}\\widehat{A}(\\omega ){\\rho }^{\\frac{1}{2}}\\right)}^{\\dagger }[\\cdot ]{\\rho }^{-\\frac{1}{2}}\\widehat{A}(\\omega ){\\rho }^{\\frac{1}{2}}{\\rm{d}}\\omega \\,\\,\\,\\,\\text{(by definition)}\\\\ \\,\\,=\\,{\\int }_{-{\\rm{\\infty }}}^{{\\rm{\\infty }}}\\gamma (\\omega ){{\\rm{e}}}^{\\beta \\omega +\\frac{{{\\sigma }}^{2}{\\beta }^{2}}{2}}\\widehat{A}(-\\omega -{{\\sigma }}^{2}\\beta )[\\cdot ]\\widehat{A}{(-\\omega -{{\\sigma }}^{2}\\beta )}^{\\dagger }{\\rm{d}}\\omega ,\\end{array}$$<\/p>\n \n (14)\n <\/p>\n showing that we may compensate for the uncertainty by imposing a shifted KMS condition (equation (12<\/a>))<\/p>\n $$\\begin{array}{c}\\gamma (\\omega ){{\\rm{e}}}^{\\beta \\omega +\\frac{{{\\sigma }}^{2}{\\beta }^{2}}{2}}\\,=\\,\\gamma (\\,-\\,\\omega -{{\\sigma }}^{2}{\\beta }^{2}),\\,\\,\\,\\text{such that}\\\\ \\,\\,(14)\\,=\\,{\\int }_{-{\\rm{\\infty }}}^{{\\rm{\\infty }}}\\gamma (\\,-\\,\\omega -{{\\sigma }}^{2}\\beta )\\widehat{A}(\\,-\\,\\omega -{{\\sigma }}^{2}\\beta )[\\cdot ]\\widehat{A}{(-\\omega -{{\\sigma }}^{2}\\beta )}^{\\dagger }{\\rm{d}}\\omega ={\\mathcal{T}}\\,[\\cdot ].\\end{array}$$<\/p>\n \n (15)\n <\/p>\n Therefore, we can take any transition weight function satisfying the KMS condition \u03b30(\u03bd)e\u03b2\u03bd\u2009=\u2009\u03b30(\u2212\u03bd) and pretend we underestimated the Bohr frequency by \\(\\frac{{\\sigma }^{2}\\beta }{2}\\), that is,\u00a0substitute \\(\\nu \\leftarrow {\\omega }_{+}:= \\omega +\\frac{{\\sigma }^{2}\\beta }{2}\\) to obtain our shifted \u03b3(\u03c9). The canonical examples include the Metropolis and Glauber weights<\/p>\n $${\\gamma }^{M}(\\omega )\\,:= \\,\\exp (-\\beta \\,\\max ({\\omega }_{+},0)).\\,(\\text{shifted Metropolis})$$<\/p>\n \n (16)\n <\/p>\n $${\\gamma }^{G}(\\omega )\\,:= \\,\\frac{1}{1+{{\\rm{e}}}^{\\beta {\\omega }_{+}}}.\\,(\\text{shifted Glauber})$$<\/p>\n \n (17)\n <\/p>\n See ref.\u200933<\/a> for the original step-by-step detailed derivation; the above streamlined derivation draws partly from subsequent simplification and development (Lemma 7.1 of ref.\u200952<\/a> and Lemma IX.2 of ref. 49<\/a>).<\/p>\n Under the natural normalization \u03b3(\u03c9)\u2009\u2208\u2009[0,\u20091], the shifted symmetry (equation (15<\/a>)) implies a low transition rate \\(\\gamma (0)\\le \\exp (-{\\sigma }^{2}{\\beta }^{2}\/2)\\) around \u03c9\u2009=\u20090. Therefore, to avoid unnecessary idling of the process, it is imperative to choose an uncertainty \u03c3\u2009\u2264\u20091\/\u03b2 not exceeding the temperature; in the main text, we have simply set \\(\\sigma =\\frac{1}{\\beta }\\). This is why implementing a single step uses\u00a0~\u03b2 Hamiltonian simulation time in our construction. By contrast, for classical or commuting systems with gapped periodic spectrum, the uncertainty can be discretized, and there is no need to scale the Hamiltonian simulation time with \u03b2.<\/p>\n Achieving full detailed balance by tuning the coherent term<\/p>\n Even if the transition part (equation (9<\/a>)) satisfies detailed balance exactly, the decay part D of the Lindbladian may still break detailed balance when it does not commute with the Hamiltonian H.<\/p>\n A second insight of our construction is to prescribe and efficiently implement a dedicated coherent term C that perfectly cancels out the deviation from quantum detailed balance in a uniquely quantum way, as shown by the following lemma.<\/p>\n Lemma 1: prescribing the coherent term<\/p>\n (Lemma II.1 of ref.\u200933<\/a>) For any full-rank density operator \\(0\\,\\prec \\,\\rho \\in {{\\mathbb{C}}}^{d\\times d}\\) and Hermitian operator \\(D\\in {{\\mathbb{C}}}^{d\\times d}\\), there is a unique Hermitian operator \\(C\\in {{\\mathbb{C}}}^{d\\times d}\\) (up to adding any scalar multiples of the identity I) such that the superoperator<\/p>\n $$-\\frac{1}{2}\\{D,\\cdot \\}-{\\rm{i}}[C,\\cdot ]$$<\/p>\n \n (18)\n <\/p>\n satisfies \u03c1-detailed balance. For a Gibbs state \\(\\rho \\propto \\exp (-\\beta H)\\), we can express C as<\/p>\n $$\\begin{array}{c}\\,C=\\sum _{\\nu \\in B(H)}\\frac{{\\rm{i}}}{2}\\,\\tanh \\left(\\frac{\\beta \\nu }{4}\\right){D}_{\\nu }\\\\ \\text{where}\\,\\,{D}_{\\nu }:= \\sum _{{E}_{i}-{E}_{j}=\\nu }|{\\psi }_{i}\\rangle \\langle {\\psi }_{i}|D|{\\psi }_{j}\\rangle \\langle {\\psi }_{j}|.\\end{array}$$<\/p>\n \n (19)\n <\/p>\n Proof<\/p>\n The detailed balance condition is equivalent to the following matrix being self-adjoint:<\/p>\n $$\\begin{array}{c}{\\rho }^{-\\frac{1}{4}}\\left(-\\frac{D}{2}-{\\rm{i}}C\\right){\\rho }^{\\frac{1}{4}}\\,=\\,\\frac{1}{2}\\,{\\rho }^{-\\frac{1}{4}}\\left(\\sum _{\\nu \\in B(H)}{D}_{\\nu }\\left(\\tanh \\left(\\frac{\\beta \\nu }{4}\\right)-1\\right)\\right){\\rho }^{\\frac{1}{4}}\\\\ \\,\\,\\,\\,\\,\\,\\,=\\,-\\frac{1}{2}\\sum _{\\nu \\in B(H)}\\frac{{D}_{\\nu }}{\\cosh \\left(\\frac{\\beta \\nu }{4}\\right)}.\\end{array}$$<\/p>\n This is self-adjoint because \\({({D}_{\\nu })}^{\\dagger }={D}_{-\\nu }\\), \\(\\cosh (x)=\\cosh (-x)\\), and B(H)\u2009=\u2009\u2212B(H). \u25a1<\/p>\n See also ref.\u200953<\/a> on prescribing fixed points in general, without necessarily assuming detailed balance. The coherent term C is a Hermitian matrix obtained by reweighing the given D operator with the profile \\({\\rm{i}}\\,\\tanh (\\beta \\nu \/4)\/2\\) on each Bohr frequency component D\u03bd. The coherent term is completely determined by \u03c1 and D (ref.\u200954<\/a>), and we found a particularly simple and useful closed-form representation in the time domain (recall \\({\\rm{sinhc}}(x)\\,:= \\,\\sinh (x)\/x\\) for real x\u2009\u2260\u20090)<\/p>\n $$C={\\int }_{-\\infty }^{\\infty }\\frac{{\\rm{i}}}{\\sinh (2{\\rm{\\pi }}t)}({{\\rm{e}}}^{{\\rm{i}}\\beta Ht}D{{\\rm{e}}}^{-{\\rm{i}}\\beta Ht}-D){\\rm{d}}t=\\frac{\\beta }{2{\\rm{\\pi }}}{\\int }_{-\\infty }^{\\infty }\\frac{-1}{{\\rm{sinhc}}(2{\\rm{\\pi }}t)}{\\int }_{0}^{1}{{\\rm{e}}}^{{\\rm{i}}s\\beta Ht}[H,D]{{\\rm{e}}}^{-{\\rm{i}}s\\beta Ht}{\\rm{d}}s{\\rm{d}}t,$$<\/p>\n \n (20)\n <\/p>\n which in turn implies the algorithmic efficiency and quasi-locality of the coherent term. The above integral form and the exponential decay in t allow using the linear combination of unitaries technique to implement a block encoding of C (refs.\u200933<\/a>,55<\/a>) by truncating and discretizing the time integral. For the Davies generator (equation (5<\/a>)), we have that [H,\u2009D]\u2009=\u20090. Therefore, the dissipative part readily satisfies detailed balance, and D\u03bd\u2009=\u20090 for all \u03bd\u2009\u2260\u20090, implying that the above coherent term prescription simply vanishes.<\/p>\n Quasi-locality<\/p>\n A salient feature of our construction is the quasi-locality of the Lindbladian terms, inherited from the operator Fourier transform in systems that feature a Lieb\u2013Robinson bound. For all \\(\\omega \\in {\\mathbb{R}}\\) and geometrically local jump A, truncating the Hamiltonian to a local Hamiltonian patch H\u2113 within distance \u2113 from the jump yields an error<\/p>\n $$\\left\\Vert {\\widehat{A}}_{(H)}(\\omega )-{\\widehat{A}}_{({H}_{{\\ell }})}(\\omega )\\right\\Vert \\le \\frac{1}{\\sqrt{2{\\rm{\\pi }}}}{\\int }_{-\\infty }^{\\infty }\\left\\Vert {{\\rm{e}}}^{{\\rm{i}}Ht}A{{\\rm{e}}}^{-{\\rm{i}}Ht}-{{\\rm{e}}}^{{\\rm{i}}{H}_{{\\ell }}t}A{{\\rm{e}}}^{-{\\rm{i}}{H}_{{\\ell }}t}\\right\\Vert |\\,f(t)|{\\rm{d}}t.$$<\/p>\n \n (21)\n <\/p>\n For a wide variety of local Hamiltonian systems, off-the-shelf Lieb\u2013Robinson bounds30<\/a> for the Heisenberg dynamics eiHtAe\u2212iHt state that the integrand in equation (21<\/a>) exponentially reduces with the distance \u2113 but degrades with the evolution time t. As the Gaussian weight function f(t) of equation (13<\/a>) effectively dampens the integral to small values of t\u2009~\u20091\/\u03c3, the operator \\({\\widehat{A}}_{(H)}(\\omega )\\) is well-approximated by a Hamiltonian patch with radius scaling with the energy uncertainty \u03c3, which can be independent of the system size. See, for example, Appendix A of ref.\u200949<\/a> for a quantitative estimate for both the transition and coherent parts. Of course, after truncation, the resulting strictly local Lindbladian may no longer satisfy exact quantum detailed balance.<\/p>\n Fixed point and its uniqueness<\/p>\n The detailed balance condition (equation (10<\/a>)) directly implies that the Gibbs state \u03c1\u00a0is a fixed point\u00a0for the Lindbladian \\({\\mathcal{L}}\\). <\/p>\n $${\\mathcal{L}}[\\,\\rho ]={\\rho }^{\\frac{1}{2}}{{\\mathcal{L}}}^{\\dagger }\\,[{\\rho }^{-\\frac{1}{2}}({\\rho }^{\\frac{1}{2}}){\\rho }^{-\\frac{1}{2}}]\\,{\\rho }^{\\frac{1}{2}}={\\rho }^{\\frac{1}{2}}{{\\mathcal{L}}}^{\\dagger }[I]{\\rho }^{\\frac{1}{2}}=0,$$<\/p>\n where the last equality used the trace-preservation property \\({{\\mathcal{L}}}^{\\dagger }[I]=0\\) of Lindbladians.<\/p>\n We now explain why the Gibbs state is the unique fixed point whenever the set of jump operators has no (nontrivial) invariant subspaces, which holds, for example, when the jumps include all single-site Pauli X and Z operators. Decomposing each jump operator by the operator Fourier transform \\({A}^{a}\\propto {\\int }_{-\\infty }^{\\infty }{\\widehat{A}}^{a}(\\omega ){\\rm{d}}\\omega \\) cannot create new invariant subspaces; likewise, multiplying by strictly positive transition weights \u03b3(\u03c9) cannot. Thus, the resulting set of Lindblad operators \\(\\sqrt{\\gamma (\\omega )}{\\widehat{A}}^{a}(\\omega )\\) has no invariant subspaces, which is known56<\/a> to imply the uniqueness of the fixed point. But this uniqueness argument says little about the quantitative convergence rate, and the mixing times can depend on the particular Hamiltonian.<\/p>\n A single-qubit example<\/p>\n Let us illustrate the details of our Lindbladian construction through a pedagogical example. Consider a single-qubit Hamiltonian with a single jump:<\/p>\n $$\\begin{array}{l}H=Z=\\left[\\begin{array}{ll}1 & 0\\\\ 0 & -1\\end{array}\\right]=| 0\\rangle \\langle 0| -| 1\\rangle \\langle 1| \\quad \\text{and}\\\\ A=X=\\left[\\begin{array}{ll}0 & 1\\\\ 1 & 0\\end{array}\\right]=| 1\\rangle \\langle 0| +| 0\\rangle \\langle 1| .\\end{array}$$<\/p>\n \n (22)\n <\/p>\n The eigenvalues of the Hamiltonian are\u00a0\u00b11, and the Bohr frequencies, their differences, are B(H)\u2009=\u2009{2,\u20090,\u2009\u22122}. We can decompose the jump by the Bohr frequencies, into the \u03bd\u2009=\u20092 and \u03bd\u2009=\u2009\u22122 components as follows:<\/p>\n $${A}_{2}=| 0\\rangle \\langle 1| \\,\\,\\text{and}\\,\\,{A}_{-2}=| 1\\rangle \\langle 0| .$$<\/p>\n \n (23)\n <\/p>\n We can then directly obtain the Davies generator for any transition weights satisfying \u03b3(2)e2\u03b2\u2009=\u2009\u03b3(\u22122). By contrast, our operator Fourier transform yields<\/p>\n $$\\widehat{A}(\\omega )=\\frac{1}{\\sqrt{{\\sigma }\\sqrt{2{\\rm{\\pi }}}}}\\left({{\\rm{e}}}^{-\\frac{{(\\omega -2)}^{2}}{4{{\\sigma }}^{2}}}{A}_{2}+{{\\rm{e}}}^{-\\frac{{(\\omega +2)}^{2}}{4{{\\sigma }}^{2}}}{A}_{-2}\\right).$$<\/p>\n \n (24)\n <\/p>\n Thus, for all \\(\\omega \\in {\\mathbb{R}}\\) the resulting \\(\\widehat{A}(\\omega )\\) has contributions coming from both the A2 and A\u22122 components; these components only get separated in the \u03c3\u2009\u2192\u20090 limit. Nevertheless, detailed balance still holds at every finite uncertainty \u03c3 with suitable choices of \u03b3(\u03c9) (equation (15<\/a>)) and the additional coherent term.<\/p>\n Recovering the Davies generator<\/p>\n As we show here, our Lindbladian exactly recovers the Davies generator in the \u03c3\u2009\u2192\u20090 limit. Moreover, the Davies generator reduces to Glauber dynamics when the Hamiltonian H is a classical function of bitstrings and the jumps A map a bitstring to another: in this case, inputs that are a probabilistic mixture of bitstrings undergo a classical Glauber dynamics (see section \u2018Quantum MCMC by master equations for thermalization<\/a>\u2019).<\/p>\n For any bounded, continuous function \u03b30 satisfying the KMS condition \u03b30(\u03bd)e\u03b2\u03bd\u2009=\u2009\u03b30(\u2212\u03bd), consider the shift \\({\\gamma }_{\\beta \\sigma }(\\omega )\\,:= \\,{\\gamma }_{0}\\left(\\omega +\\frac{{\\beta }^{2}{\\sigma }^{2}}{2}\\right),\\) then the transition part (similarly for the decay term) converges to<\/p>\n $$\\mathop{\\mathrm{lim}}\\limits_{\\sigma \\to 0}{\\int }_{-\\infty }^{\\infty }{\\gamma }_{\\beta \\sigma }(\\omega )\\widehat{A}(\\omega )[\\cdot ]\\widehat{A}{(\\omega )}^{\\dagger }{\\rm{d}}\\omega =\\sum _{\\nu \\in B(H)}{\\gamma }_{0}(\\nu ){A}_{\\nu }[\\cdot ]{({A}_{\\nu })}^{\\dagger }.$$<\/p>\n When we expand the left-hand side as a sum over \\({A}_{{\\nu }_{1}}[\\cdot ]{({A}_{{\\nu }_{2}})}^{\\dagger }\\), the coefficients for each \u03bd1,\u00a0\u03bd2<\/p>\n $${\\mathrm{lim}}_{\\sigma \\to 0}{\\int }_{-\\infty }^{\\infty }\\frac{{\\gamma }_{\\beta \\sigma }(\\omega )}{\\sigma \\sqrt{2{\\rm{\\pi }}}}{{\\rm{e}}}^{-\\frac{{(\\omega -{\\nu }_{1})}^{2}}{4{\\sigma }^{2}}}{{\\rm{e}}}^{-\\frac{{(\\omega -{\\nu }_{2})}^{2}}{4{\\sigma }^{2}}}{\\rm{d}}\\omega ={\\mathbb{1}}({\\nu }_{1}={\\nu }_{2}){\\gamma }_{0}({\\nu }_{1}),$$<\/p>\n approach the corresponding value in the Davies generator. By contrast, the coherent term vanishes, because the decay term commutes with the Hamiltonian in the \u03c3\u2009\u2192\u20090 limit and tanh(0)\u2009=\u20090 (see Lemma 1).<\/p>\n Implementation<\/p>\n In this section, we give low-level quantum algorithmic implementations of our Lindbladian dynamics.<\/p>\n A phase estimation for operators<\/p>\n The most general form of the discrete operator Fourier transform is a subroutine similar to quantum phase estimation: it combines controlled Hamiltonian simulation with the quantum Fourier transform, as shown in Fig. 3a<\/a>, acting on the time\u2013frequency and system registers. To approximately implement the (continuous) operator Fourier transform, we need to discretize the time integral by introducing a time mesh \\(\\bar{t}\\in {S}_{{t}_{0}}\\) and a corresponding frequency mesh \\(\\bar{\\omega }\\in {S}_{{\\omega }_{0}}\\) for the QFT, each having M points, resulting in the discretized operation<\/p>\n $$\\begin{array}{l}\\mathop{\\underbrace{{\\sum }_{\\bar{t}\\in {S}_{{t}_{0}}}f(\\bar{t})| \\bar{t}\\rangle \\otimes A\\to {\\sum }_{\\bar{\\omega }\\in {S}_{{\\omega }_{0}}}| \\bar{\\omega }\\rangle \\otimes \\widehat{A}(\\bar{\\omega })}}\\limits_{{\\rm{discrete}}\\,{\\rm{operator}}\\,{\\rm{Fourier}}\\,{\\rm{transform}}}\\\\ \\text{where}\\,\\widehat{A}(\\bar{\\omega })\\,:= \\,\\frac{1}{\\sqrt{M}}\\sum _{\\bar{t}\\in {S}_{{t}_{0}}}{e}^{-i\\bar{\\omega }\\bar{t}}f(\\bar{t}){e}^{iHt}A{e}^{-iHt}.\\end{array}$$<\/p>\n \n (25)\n <\/p>\n As given in Corollary C.2 of ref.\u200922<\/a>, choosing<\/p>\n $$\\begin{array}{l}\\,{\\omega }_{0}=2\\sigma \\sqrt{\\frac{2{\\rm{\\pi }}}{M}},\\,{t}_{0}=\\frac{1}{2\\sigma }\\sqrt{\\frac{2{\\rm{\\pi }}}{M}},\\\\ {S}^{\\lceil M\\rfloor }:= \\,\\{-\\lceil (M-1)\/2\\rceil ,\\ldots ,-1,0,1,\\ldots ,\\lfloor (M-1)\/2\\rfloor \\},\\end{array}$$<\/p>\n \n (26)\n <\/p>\n $${\\rm{and}}\\quad {S}_{{\\omega }_{0}}:= {\\omega }_{0}\\cdot {S}^{\\lceil M\\rfloor },\\quad {S}_{{t}_{0}}:= {t}_{0}\\cdot {S}^{\\lceil M\\rfloor },$$<\/p>\n \n (27)\n <\/p>\n the above equation (25<\/a>) recovers the advertised continuous operator Fourier transform (equation (13<\/a>)) in the M\u2009\u2192\u2009\u221e limit. As f is a smooth Gaussian function, we can achieve any finite precision \u03f5-approximation of the dissipative part of \\({\\mathcal{L}}\\) with a moderately scaling dimension \\(M \\sim {\\rm{Poly}}(\\parallel H\\parallel ,\\beta ,1\/{\\epsilon },\\sigma +1\/\\sigma ,| A| )\\), which requires only \\(\\log (M)=\\widetilde{{\\mathcal{O}}}(1)\\)-many ancilla qubits, and no more than \\({\\mathcal{O}}(\\sqrt{\\log (1\/{\\epsilon })}\\,\/\\sigma )\\) (controlled) Hamiltonian simulation time, because of the truncation of the Gaussian tail.<\/p>\n A general-purpose Lindbladian simulation algorithm<\/p>\n Another key algorithmic component is an improved black-box Lindbladian simulation subroutine. It achieves the sought nearly linear dependence on t even for Lindbladians with high Kraus rank, as long as the Lindblad operators are given in a block-encoded form, which represents an improvement on previous Lindbladian simulation algorithms51<\/a>,57<\/a>.<\/p>\n Definition 1: block encoding for Lindblad operators<\/p>\n (Definition I.2 of ref.\u200922<\/a>) We say that a unitary U is a block encoding for Lindblad operators {Lj:\u00a0j\u2009\u2208\u2009J}, if<\/p>\n $$(\\langle {0}^{b}|\\otimes I)\\cdot U\\cdot (|{0}^{q}\\rangle \\otimes I)=\\sum _{j\\in J}|\\,j\\rangle \\,\\otimes {L}_{j}\\,\\text{for some}\\,b\\le q\\in {\\mathbb{N}}.$$<\/p>\n In particular, discretized operator Fourier transforms (equation (25<\/a>)) naturally give block encodings of this form, where J refers to the discretized set of frequency labels \\({S}_{{\\omega }_{0}}\\). Given a block encoding U of the Lindblad operators as above, we can directly obtain a block encoding of the decay term D by a single use of U and U\u2020 by standard multiplication of block-encoded matrices11<\/a>. As a consequence of the exponentially decaying tails in the integral representation of equation (20<\/a>), by using an additional\u00a0~\u00a0\u03b2-time controlled Hamiltonian simulation and the linear combination of unitaries (LCU) technique55<\/a>, we can also obtain a sufficiently good approximate block encoding of the coherent term C.<\/p>\n We begin with a brief analysis of the first-order simulation circuit shown in Fig. 3b<\/a>. Assuming the system register is in the pure state |\u03c8\u27e9, the first three gates act as follows (for simplicity, we drop the superscripts Aa\u2009\u2190\u2009A):<\/p>\n $$\\begin{array}{cc}|0\\rangle |{0}^{q}\\rangle |\\psi \\rangle & \\mathop{\\to }\\limits^{(1)}|0\\rangle \\otimes \\sum _{\\bar{\\omega }\\in {S}_{{\\omega }_{0}}}|\\bar{\\omega }\\rangle \\otimes \\hat{A}(\\bar{\\omega })|\\psi \\rangle \\\\ & \\mathop{\\to }\\limits^{(2)}\\sum _{\\bar{\\omega }\\in {S}_{{\\omega }_{0}}}(\\mathop{\\underbrace{\\sqrt{1-\\delta \\gamma (\\bar{\\omega })}}}\\limits_{1-\\frac{\\delta \\gamma (\\bar{\\omega })}{2}+{\\mathcal{O}}({\\delta }^{2})}|0\\rangle +\\sqrt{\\delta \\gamma (\\bar{\\omega })}|1\\rangle )|\\bar{\\omega }\\rangle \\hat{A}(\\bar{\\omega })|\\psi \\rangle \\\\ & \\mathop{\\to }\\limits^{(3)}|0\\rangle |{0}^{q}\\rangle \\left({I}-\\frac{\\delta }{2}D\\right)|\\psi \\rangle +|1\\rangle \\sum _{\\bar{\\omega }\\in {S}_{{\\omega }_{0}}}\\sqrt{\\delta \\gamma (\\bar{\\omega })}|\\bar{\\omega }\\rangle \\hat{A}(\\bar{\\omega })|\\psi \\rangle -\\frac{\\delta }{2}|0\\rangle \\otimes |{0}^{q}\\perp \\rangle +{\\mathcal{O}}({\\delta }^{2}),\\end{array}$$<\/p>\n \n (28)\n <\/p>\n where \\(| {0}^{q}\\perp \\rangle \\) is a quantum state such that \\(\\parallel | {0}^{q}\\perp \\rangle \\parallel \\le 1\\) and \\((\\langle {0}^{q}| \\otimes I)\\cdot | {0}^{q}\\perp \\rangle =0\\), see Theorem III.1 of ref.\u200922<\/a>, for details.<\/p>\n Let \\(| {\\psi }^{{\\prime} }\\rangle \\) be the resulting state in equation (28<\/a>), and let us denote the dissipative part\u00a0\\({{\\mathcal{L}}}^{{\\prime} }:= {{\\mathcal{L}}}^{a}+{\\rm{i}}[{C}^{a},\\cdot ]\\). Tracing out the first q\u2009+\u20091 qubits, we get that \\(| {\\psi }^{{\\prime} }\\rangle \\) is \\({\\mathcal{O}}({\\delta }^{2})\\)-close to the desired state, ignoring the coherent term. We now show that<\/p>\n $${\\parallel ({\\mathcal{I}}+\\delta {{\\mathcal{L}}}^{{\\prime} })[| \\psi \\rangle \\langle \\psi | ]-{{\\rm{Tr}}}_{q+1}[| {\\psi }^{{\\prime} }\\rangle \\langle {\\psi }^{{\\prime} }| ]\\parallel }_{1}={\\mathcal{O}}({\\delta }^{2})$$<\/p>\n \n (29)\n <\/p>\n by observing that<\/p>\n $$\\begin{array}{c}{{\\rm{T}}{\\rm{r}}}_{q+1}[|{\\psi }^{{\\prime} }\\rangle \\langle {\\psi }^{{\\prime} }|]\\,=\\,{{\\rm{T}}{\\rm{r}}}_{q}[(\\langle 0|\\otimes I)\\cdot |{\\psi }^{{\\prime} }\\rangle \\langle {\\psi }^{{\\prime} }|\\cdot (|0\\rangle \\otimes I)]+{{\\rm{T}}{\\rm{r}}}_{q}[(\\langle 1|\\otimes I)\\cdot |{\\psi }^{{\\prime} }\\rangle \\langle {\\psi }^{{\\prime} }|\\cdot (|1\\rangle \\otimes I)]\\\\ \\,\\,\\,\\,\\,\\,=\\,\\left(I-\\frac{\\delta }{2}D\\right)|\\psi \\rangle \\langle \\psi |\\left(I-\\frac{\\delta }{2}D\\right)\\\\ \\,\\,\\,\\,\\,\\,+\\,\\delta \\sum _{\\overline{\\omega }\\in {S}_{{\\omega }_{0}}}\\gamma (\\overline{\\omega })\\hat{A}(\\overline{\\omega })|\\psi \\rangle \\langle \\psi |\\hat{A}{(\\overline{\\omega })}^{\\dagger }+{\\mathcal{O}}({\\delta }^{2})\\\\ \\,\\,\\,\\,\\,\\,=\\,({\\mathcal{I}}+\\delta {{\\mathcal{L}}}^{{\\prime} })[|\\psi \\rangle \\langle \\psi |]+{\\mathcal{O}}({\\delta }^{2}).\\end{array}$$<\/p>\n Convexity and the triangle inequality imply that equation (29<\/a>) also holds for mixed input states. By including the \u03b4-time Hamiltonian simulation for Ca, we get an \\({\\mathcal{O}}({\\delta }^{2})\\) approximation of \u03b4-time evolution by \\({{\\mathcal{L}}}^{a}\\). Once again, owing to linearity and the triangle inequality, this also implies that performing the circuit Fig. 3b<\/a> for a uniformly random jump Aa we get a \u03b4-time evolution by \\(\\frac{1}{|A|}\\sum _{a\\in A}{{\\mathcal{L}}}_{a}\\) up to \\({\\mathcal{O}}({\\delta }^{2})\\) error in trace distance, see Corollary III.1 of ref.\u200922<\/a>. Repeating the entire argument after replacing the jump operators Aa by Aa \u2297 I, we can see that actually the above results in a \\({\\mathcal{O}}({\\delta }^{2})\\)-precise implementation of the quantum channel \\(\\exp \\left(\\frac{\\delta }{|A|}\\sum _{a\\in A}{{\\mathcal{L}}}_{a}\\right)\\) in the completely bounded 1-1 superoperator norm, that is, the diamond-norm.<\/p>\n Choosing \\(\\delta =\\varTheta \\left(\\frac{{\\epsilon }}{t}\\right)\\) ensures that the error in a single time step is bounded by \\({\\mathcal{O}}\\left(\\frac{{{\\epsilon }}^{2}}{{t}^{2}}\\right)\\), and repeating the process \\(\\varTheta \\left(\\frac{{t}^{2}}{{\\epsilon }}\\right)\\)-times induces an error that is bounded by \u03f5 for the entire time-t evolution. The complexity is then \\(\\varTheta \\left(\\frac{{t}^{2}}{{\\epsilon }}\\right)\\) times the cost of implementing the circuit in Fig. 3b<\/a>.<\/p>\n Building on the compression techniques in ref.\u200957<\/a>, we can bootstrap the first-order weak-measurement circuit of Fig. 3b<\/a> by observing that for very small time steps, the circuit is very close to identity. Exploiting this by effectively only using the circuit Fig. 3b<\/a> in parts of the trajectories that are nontrivial, we can achieve almost linear scaling in t and polylogarithmic scaling in the desired diamond-norm accuracy of the simulation. However, we need to apply some modifications, as it turns out the original compressed measurement scheme described in ref.\u200957<\/a> does not work as intended. Thus, we use a variant of the analogous measurement scheme described in ref.\u200958<\/a> (see Appendix F of ref.\u200922<\/a>, for more details). This amounts to our ultimate near-linear-time simulation result.<\/p>\n Theorem 2: almost linear-time Lindbladian simulation<\/p>\n (Theorem III.2 of ref.\u200922<\/a>) Suppose U is a unitary block encoding of the Lindblad operators of \\({\\mathcal{L}}\\) as in Definition 1, and V is a block encoding of the coherent term C. Let t\u2009\u2265\u20091 and \u03f5\u2009\u2264\u20091\/2, then we can simulate the map \\({e}^{{\\mathcal{L}}t}\\) to error \u03f5 in diamond norm using<\/p>\n $${\\mathcal{O}}((q+\\log (t\/{\\epsilon }))\\log (t\/{\\epsilon }))\\,\\,({\\rm{resettable}})\\,{\\rm{ancilla}}\\,{\\rm{qubits}},$$<\/p>\n \n (30)\n <\/p>\n $${\\mathcal{O}}\\left(t\\frac{\\log (t\/{\\epsilon })}{\\log \\,\\log (t\/{\\epsilon })}\\right)\\,\\,({\\rm{controlled}})\\,{\\rm{uses}}\\,{\\rm{of}}\\,U,V,{U}^{\\dagger }\\,{\\rm{and}}\\,{V}^{\\dagger },$$<\/p>\n \n (31)\n <\/p>\n $$\\,\\text{and}\\,\\,{\\mathcal{O}}(t(q+1){\\rm{polylog}}(t\/{\\epsilon }))\\,\\,\\,\\text{other two-qubit gates},$$<\/p>\n \n (32)\n <\/p>\n where q is the number of ancilla qubits used for the block encodings.<\/p>\n To place the above general result in context, we give some simple bounds on the resources required for implementing our Lindbladian. For instance, if the jump operators are K different Pauli strings, then \\(q={\\log }_{2}(K)+1\\) ancilla qubits suffice for block encoding them. The overall gate complexity should be dominated by the controlled Hamiltonian simulation subroutine. Thus, we focus on estimating the required Hamiltonian simulation time per call to U and V (assuming an operator Fourier transform width \\(\\sigma =\\frac{1}{\\beta }\\)). Under the normalization \\(\\parallel \\sum _{a\\in A}{A}^{a\\dagger }{A}^{a}\\parallel \\le 1\\), the Lindblad operators from our construction (equation (7<\/a>)) can be \u03f5-accurately block-encoded by discretized operator Fourier transform using \\({\\mathcal{O}}(\\beta \\sqrt{\\log (1\/{\\epsilon })})\\) (controlled) Hamiltonian simulation time by truncating the Gaussian integrand in equation (6<\/a>). Meanwhile, as implied by Corollary III.2 of ref.\u200933<\/a>, a slightly subnormalized coherent term C\/\u03b1 for \\(\\alpha ={\\mathcal{O}}\\left(\\log \\left(\\frac{\\beta \\,\\parallel \\,H\\,\\parallel }{{\\epsilon }}\\right)\\right)\\) can be \u03f5-accurately block-encoded by LCU by using a mere \\({\\mathcal{O}}(\\beta \\,\\log (1\/{\\epsilon }))\\) (controlled) Hamiltonian simulation time. The extra subnormalization factor \u03b1 can be absorbed into the number of uses of the block encoding V using state-of-the-art block-encoded Hamiltonian simulation techniques4<\/a>,11<\/a>.<\/p>\n To approximate \\({e}^{{\\mathcal{L}}t}\\) to \u03f5-diamond distance, we can control the accumulation of errors by setting an increased accuracy goal of \\({\\rm{poly}}({\\epsilon }\/(t\\beta \\parallel H\\parallel ))\\) for the approximate block encodings U, V. This results in an overall \\({\\mathcal{O}}(\\beta \\sqrt{\\log (t\\beta \\parallel H\\parallel \/{\\epsilon })})\\) and \\({\\mathcal{O}}(\\beta {\\log }^{2}(t\\beta \\parallel H\\parallel \/{\\epsilon }))\\) (controlled) Hamiltonian simulation time overhead for implementing sufficiently accurate block encodings U and V, respectively.<\/p>\n Comparison with physically derived master equations<\/p>\n Our synthetic Lindbladian is mainly presented as an algorithmic construction. In this section, we compare it with physically motivated master equations derived from first-principles open-system calculations. Overall, we believe that our Lindbladian can serve as a self-contained toy model of thermalization. Nevertheless, for quantitative modelling of particular system\u2013bath interactions, the role of strict detailed balance is less clear, and it may be preferable to use a physically derived master equation.<\/p>\n The textbook setup in open-system thermalization24<\/a>,59<\/a> considers a system governed by a Hamiltonian HS that couples weakly to a large thermal bath with Hamiltonian HB, which are together governed by the total Hamiltonian \\({H}_{{\\rm{tot}}}={H}_{{\\rm{S}}}\\otimes {I}_{{\\rm{B}}}+{I}_{{\\rm{S}}}\\otimes {H}_{{\\rm{B}}}+\\lambda \\sum _{a\\in A}{A}^{a}\\otimes {B}^{a}\\). The Hermitian operators {Aa:\u00a0a\u2009\u2208\u2009A} act on the system (mirroring our jump operators, hence the same notation) and {Ba:\u00a0a\u2009\u2208\u2009A} acts on the bath, whereas \u03bb represents the coupling strength. Tracing out the bath, we can obtain an effective master equation governing the system dynamics under the assumptions that the thermal bath is Markovian and the coupling \u03bb is sufficiently weak.<\/p>\n The aforementioned Davies generator was originally derived in the weak-coupling limit \u03bb\u2009\u2192\u20090 (relative to all other energy scales). In this limit, the rotating-wave approximation (or the secular approximation) removes the cross terms \\({A}_{{\\nu }_{1}}[\\cdot ]{A}_{{\\nu }_{2}}^{\\dagger }\\); this perfect isolation of Bohr frequencies causes the Davies generator to satisfy detailed balance, but at the same time, makes it implausible for many-body systems with exponentially small level spacings. To focus on issues related to detailed balance, here we studied only the simplest form of the Davies generator. In principle, different jumps Aa may have different transition weights, as they correspond to the Fourier transform of the bath correlation function \\(\\langle {{\\rm{e}}}^{{\\rm{i}}{H}_{B}t}{B}^{a}{{\\rm{e}}}^{-{\\rm{i}}{H}_{B}t}{B}^{a}\\rangle \\), but for algorithmic purposes, it is natural to use universal transition weights. Also, we have studied only the dissipative part of the Davies generator (equation (5<\/a>)). However, the full Davies generator also includes a Lamb-shift term that somewhat resembles our coherent term, but depends on additional details of the bath correlation functions. As the Lamb-shift term commutes with the Hamiltonian, adding this term keeps the Gibbs state stationary.<\/p>\n Recently, there have been several attempts to derive more realistic master equations for many-body systems by avoiding the \u03bb\u2009\u2192\u20090 limit needed for the rotating-wave approximation. This essentially translates to introducing Lindblad operators with a finite energy uncertainty, or, equivalently, in the time domain, integrals over Heisenberg-evolved jump operators weighted by a finite-width window function. In particular, the coarse-grained master equation59<\/a> takes a very similar form as ours (equation (7<\/a>)), but its operator Fourier transform features a uniform integral \\(\\frac{1}{\\sqrt{2{\\rm{\\pi }}T}}{\\int }_{-\\frac{T}{2}}^{\\frac{T}{2}}{{\\rm{e}}}^{{\\rm{i}}Ht}{A}^{a}{{\\rm{e}}}^{-{\\rm{i}}Ht}{{\\rm{e}}}^{-{\\rm{i}}\\omega t}{\\rm{d}}t\\), and its coherent term also resembles ours (equation (19<\/a>)) but does not seem to strictly enforce detailed balance. The universal Lindblad equation60<\/a> combines (square root of) the transition weights into the operator Fourier transforms and is closer to subsequent constructions61<\/a>,62<\/a>. Very recently63<\/a> derived a master equation with exact KMS-detailed balance using slightly different operator Fourier transform weights compared with ours. All the above recent master equations feature some finite energy uncertainty, derived from various system\u2013bath parameters, that parallels our tunable Gaussian width \u03c3.<\/p>\n","protected":false},"excerpt":{"rendered":"We provide the explicit form of our construction and explain how exact detailed balance can be achieved through…\n","protected":false},"author":3,"featured_media":306723,"comment_status":"","ping_status":"","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[25],"tags":[26264,26265,10046,50568,10047,492,18678,159,63805,67,132,68],"class_list":{"0":"post-306722","1":"post","2":"type-post","3":"status-publish","4":"format-standard","5":"has-post-thumbnail","7":"category-physics","8":"tag-computational-science","9":"tag-computer-science","10":"tag-humanities-and-social-sciences","11":"tag-information-theory-and-computation","12":"tag-multidisciplinary","13":"tag-physics","14":"tag-quantum-simulation","15":"tag-science","16":"tag-thermodynamics","17":"tag-united-states","18":"tag-unitedstates","19":"tag-us"},"share_on_mastodon":{"url":"https:\/\/pubeurope.com\/@us\/115381091920798541","error":""},"_links":{"self":[{"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/posts\/306722","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/users\/3"}],"replies":[{"embeddable":true,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/comments?post=306722"}],"version-history":[{"count":0,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/posts\/306722\/revisions"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/media\/306723"}],"wp:attachment":[{"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/media?parent=306722"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/categories?post=306722"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/tags?post=306722"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}