familyi={P\bibinithyphendelim L\bibinitperiod},
given={Anatoly\bibnamedelima Yu.},
giveni={A\bibinitperiod\bibinitdelim Y\bibinitperiod}}}%
given={Andrew\bibnamedelima J.},
giveni={A\bibinitperiod\bibinitdelim J\bibinitperiod}}}%
given={Paul\bibnamedelima I.},
giveni={P\bibinitperiod\bibinitdelim I\bibinitperiod}}}%
given={P.\bibnamedelimi Aaron},
giveni={P\bibinitperiod\bibinitdelim A\bibinitperiod}}}%
given={Allison\bibnamedelimb J.\bibnamedelimi R.},
giveni={A\bibinitperiod\bibinitdelim J\bibinitperiod\bibinitdelim R\bibinitperiod}}}%
given={Loren\bibnamedelima J.},
giveni={L\bibinitperiod\bibinitdelim J\bibinitperiod}}}%
given={Mark\bibnamedelima H.},
giveni={M\bibinitperiod\bibinitdelim H\bibinitperiod}}}%
given={Jed\bibnamedelima D.},
giveni={J\bibinitperiod\bibinitdelim D\bibinitperiod}}}%
given={Jeremy\bibnamedelima P.},
giveni={J\bibinitperiod\bibinitdelim P\bibinitperiod}}}%
given={Mark\bibnamedelima W.},
giveni={M\bibinitperiod\bibinitdelim W\bibinitperiod}}}%
given={Mohammad\bibnamedelima H.},
giveni={M\bibinitperiod\bibinitdelim H\bibinitperiod}}}%
\strng{namehash}{9fed10ecd9131a424b47594c63ce52f8}
\strng{fullhash}{70f539fc440fb2366f65bbb70557a443}
\strng{bibnamehash}{9fed10ecd9131a424b47594c63ce52f8}
\strng{authorbibnamehash}{9fed10ecd9131a424b47594c63ce52f8}
\strng{authornamehash}{9fed10ecd9131a424b47594c63ce52f8}
\strng{authorfullhash}{70f539fc440fb2366f65bbb70557a443}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{journaltitle}{arXiv:1911.03446}
\field{title}{Scaling advantage in quantum simulation of geometrically frustrated magnets}
\field{year}{2019}
\verb{urlraw}
\verb http://arXiv.org/abs/1911.03446
\endverb
\verb{url}
\verb http://arXiv.org/abs/1911.03446
\endverb
\endentry
\entry{Gardas2018}{article}{}
\name{author}{4}{}{%
given={Wojciech\bibnamedelima H.},
giveni={W\bibinitperiod\bibinitdelim H\bibinitperiod}}}%
\strng{namehash}{e0001ae412e438bd813cd87a4c1f4c6b}
\strng{fullhash}{e0001ae412e438bd813cd87a4c1f4c6b}
\strng{bibnamehash}{e0001ae412e438bd813cd87a4c1f4c6b}
\strng{authorbibnamehash}{e0001ae412e438bd813cd87a4c1f4c6b}
\strng{authornamehash}{e0001ae412e438bd813cd87a4c1f4c6b}
\strng{authorfullhash}{e0001ae412e438bd813cd87a4c1f4c6b}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{issn}{20452322}
\field{journaltitle}{Sci. Rep.}
\field{number}{1}
\field{title}{{Defects in Quantum Computers}}
\field{volume}{8}
\field{year}{2018}
\field{pages}{4539}
\range{pages}{1}
\verb{doi}
\verb 10.1038/s41598-018-22763-2
\endverb
\endentry
\entry{Nishimura2020}{article}{}
\name{author}{3}{}{%
given={Helmut\bibnamedelima G.},
giveni={H\bibinitperiod\bibinitdelim G\bibinitperiod}}}%
\strng{namehash}{f25a4ed184db84dd3861b14173952264}
\strng{fullhash}{f25a4ed184db84dd3861b14173952264}
\strng{bibnamehash}{f25a4ed184db84dd3861b14173952264}
\strng{authorbibnamehash}{f25a4ed184db84dd3861b14173952264}
\strng{authornamehash}{f25a4ed184db84dd3861b14173952264}
\strng{authorfullhash}{f25a4ed184db84dd3861b14173952264}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{issn}{2469-9926}
\field{journaltitle}{Phys. Rev. A}
\field{number}{4}
\field{title}{{Griffiths-McCoy singularity on the diluted Chimera graph: Monte Carlo simulations and experiments on quantum hardware}}
\field{volume}{102}
\field{year}{2020}
\verb{doi}
\verb 10.1103/physreva.102.042403
\endverb
\endentry
\entry{King2020}{article}{}
\name{author}{5}{}{%
given={Andrew\bibnamedelima D.},
giveni={A\bibinitperiod\bibinitdelim D\bibinitperiod}}}%
given={Edward\bibnamedelima D.},
giveni={E\bibinitperiod\bibinitdelim D\bibinitperiod}}}%
familyi={P\bibinithyphendelim L\bibinitperiod},
familyi={L\bibinithyphendelim B\bibinitperiod},
\strng{namehash}{ddeaac7442877612eef52d7970778ad7}
\strng{fullhash}{ddeaac7442877612eef52d7970778ad7}
\strng{bibnamehash}{ddeaac7442877612eef52d7970778ad7}
\strng{authorbibnamehash}{ddeaac7442877612eef52d7970778ad7}
\strng{authornamehash}{ddeaac7442877612eef52d7970778ad7}
\strng{authorfullhash}{ddeaac7442877612eef52d7970778ad7}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{eprinttype}{arXiv}
\field{journaltitle}{arXiv:2007.10555}
\field{month}{7}
\field{title}{{Quantum Artificial Spin Ice}}
\field{year}{2020}
\verb{eprint}
\verb 2007.10555
\endverb
\verb{urlraw}
\verb https://arxiv.org/abs/2007.10555 http://arxiv.org/abs/2007.10555
\endverb
\verb{url}
\verb https://arxiv.org/abs/2007.10555%20http://arxiv.org/abs/2007.10555
\endverb
\endentry
\entry{Kairys2020}{article}{}
\name{author}{7}{}{%
given={Andrew\bibnamedelima D},
giveni={A\bibinitperiod\bibinitdelim D\bibinitperiod}}}%
given={Travis\bibnamedelima S},
giveni={T\bibinitperiod\bibinitdelim S\bibinitperiod}}}%
\strng{namehash}{3fc44ef289d292232e2c5bd2c3f18bb6}
\strng{fullhash}{3d42d4cf3ddae0c3a2ff5373fc05577e}
\strng{bibnamehash}{3fc44ef289d292232e2c5bd2c3f18bb6}
\strng{authorbibnamehash}{3fc44ef289d292232e2c5bd2c3f18bb6}
\strng{authornamehash}{3fc44ef289d292232e2c5bd2c3f18bb6}
\strng{authorfullhash}{3d42d4cf3ddae0c3a2ff5373fc05577e}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{eprinttype}{arXiv}
\field{journaltitle}{arXiv:2003.01019}
\field{title}{{Simulating the Shastry-Sutherland Ising Model using Quantum Annealing}}
\field{year}{2020}
\verb{eprint}
\verb 2003.01019
\endverb
\verb{urlraw}
\verb http://arxiv.org/abs/2003.01019
\endverb
\verb{url}
\verb http://arxiv.org/abs/2003.01019
\endverb
\endentry
\entry{kirkpatrick:sa}{article}{}
\name{author}{3}{}{%
given={C\bibnamedelima Daniel},
giveni={C\bibinitperiod\bibinitdelim D\bibinitperiod}}}%
given={Mario\bibnamedelima P},
giveni={M\bibinitperiod\bibinitdelim P\bibinitperiod}}}%
\list{publisher}{1}{%
{American association for the advancement of science}%
\strng{namehash}{c5fe22fb669c5a6c1b2d8ac99dca41f5}
\strng{fullhash}{c5fe22fb669c5a6c1b2d8ac99dca41f5}
\strng{bibnamehash}{c5fe22fb669c5a6c1b2d8ac99dca41f5}
\strng{authorbibnamehash}{c5fe22fb669c5a6c1b2d8ac99dca41f5}
\strng{authornamehash}{c5fe22fb669c5a6c1b2d8ac99dca41f5}
\strng{authorfullhash}{c5fe22fb669c5a6c1b2d8ac99dca41f5}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{journaltitle}{science}
\field{number}{4598}
\field{title}{Optimization by simulated annealing}
\field{volume}{220}
\field{year}{1983}
\field{pages}{671\bibrangedash 680}
\range{pages}{10}
\endentry
\entry{Brady2020}{article}{}
\name{author}{5}{}{%
given={Lucas\bibnamedelima T},
giveni={L\bibinitperiod\bibinitdelim T\bibinitperiod}}}%
given={Christopher\bibnamedelima L},
giveni={C\bibinitperiod\bibinitdelim L\bibinitperiod}}}%
given={Alexey\bibnamedelima V},
giveni={A\bibinitperiod\bibinitdelim V\bibinitperiod}}}%
\strng{namehash}{02241706ee2b45dacd2ccf5257af1e8a}
\strng{fullhash}{02241706ee2b45dacd2ccf5257af1e8a}
\strng{bibnamehash}{02241706ee2b45dacd2ccf5257af1e8a}
\strng{authorbibnamehash}{02241706ee2b45dacd2ccf5257af1e8a}
\strng{authornamehash}{02241706ee2b45dacd2ccf5257af1e8a}
\strng{authorfullhash}{02241706ee2b45dacd2ccf5257af1e8a}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{journaltitle}{arXiv preprint arXiv:2003.08952}
\field{title}{Optimal protocols in quantum annealing and qaoa problems}
\field{year}{2020}
\endentry
\entry{Mbeng2019}{article}{}
\name{author}{3}{}{%
given={Glen\bibnamedelima Bigan},
giveni={G\bibinitperiod\bibinitdelim B\bibinitperiod}}}%
\strng{namehash}{02b7947b7bfed8a691651797e07fe448}
\strng{fullhash}{02b7947b7bfed8a691651797e07fe448}
\strng{bibnamehash}{02b7947b7bfed8a691651797e07fe448}
\strng{authorbibnamehash}{02b7947b7bfed8a691651797e07fe448}
\strng{authornamehash}{02b7947b7bfed8a691651797e07fe448}
\strng{authorfullhash}{02b7947b7bfed8a691651797e07fe448}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{journaltitle}{arXiv:1906.08948}
\field{month}{6}
\field{title}{{Quantum Annealing: a journey through Digitalization, Control, and hybrid Quantum Variational schemes}}
\field{year}{2019}
\verb{urlraw}
\verb http://arxiv.org/abs/1906.08948
\endverb
\verb{url}
\verb http://arxiv.org/abs/1906.08948
\endverb
\endentry
\entry{Gross1984}{article}{}
\name{author}{2}{}{%
given={D.\bibnamedelimi J.},
giveni={D\bibinitperiod\bibinitdelim J\bibinitperiod}}}%
\strng{namehash}{16f9f974660fa0272d6885b714f5c569}
\strng{fullhash}{16f9f974660fa0272d6885b714f5c569}
\strng{bibnamehash}{16f9f974660fa0272d6885b714f5c569}
\strng{authorbibnamehash}{16f9f974660fa0272d6885b714f5c569}
\strng{authornamehash}{16f9f974660fa0272d6885b714f5c569}
\strng{authorfullhash}{16f9f974660fa0272d6885b714f5c569}
\field{extraname}{2}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{issn}{05503213}
\field{journaltitle}{Nucl. Phys. B}
\field{number}{4}
\field{title}{{The simplest spin glass}}
\field{volume}{240}
\field{year}{1984}
\field{pages}{431\bibrangedash 452}
\range{pages}{22}
\verb{doi}
\verb 10.1016/0550-3213(84)90237-2
\endverb
\endentry
\entry{Bapst2012}{article}{}
\name{author}{2}{}{%
\strng{namehash}{f7d684a3fe16c7da69bb1882d1f6cdb7}
\strng{fullhash}{f7d684a3fe16c7da69bb1882d1f6cdb7}
\strng{bibnamehash}{f7d684a3fe16c7da69bb1882d1f6cdb7}
\strng{authorbibnamehash}{f7d684a3fe16c7da69bb1882d1f6cdb7}
\strng{authornamehash}{f7d684a3fe16c7da69bb1882d1f6cdb7}
\strng{authorfullhash}{f7d684a3fe16c7da69bb1882d1f6cdb7}
\field{extraname}{2}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{issn}{17425468}
\field{journaltitle}{J. Stat. Mech.: Theory Exp.}
\field{number}{6}
\field{title}{{On quantum mean-field models and their quantum annealing}}
\field{volume}{2012}
\field{year}{2012}
\field{pages}{P06007}
\range{pages}{-1}
\verb{doi}
\verb 10.1088/1742-5468/2012/06/p06007
\endverb
\verb{urlraw}
\verb https://doi.org/10.1088/1742-5468/2012/06/p06007
\endverb
\verb{url}
\verb https://doi.org/10.1088/1742-5468/2012/06/p06007
\endverb
\endentry
\entry{shin2014quantum}{article}{}
\name{author}{4}{}{%
given={Seung\bibnamedelima Woo},
giveni={S\bibinitperiod\bibinitdelim W\bibinitperiod}}}%
given={John\bibnamedelima A},
giveni={J\bibinitperiod\bibinitdelim A\bibinitperiod}}}%
\strng{namehash}{5c3bfc4d726caf513f9cf032d754e0a5}
\strng{fullhash}{5c3bfc4d726caf513f9cf032d754e0a5}
\strng{bibnamehash}{5c3bfc4d726caf513f9cf032d754e0a5}
\strng{authorbibnamehash}{5c3bfc4d726caf513f9cf032d754e0a5}
\strng{authornamehash}{5c3bfc4d726caf513f9cf032d754e0a5}
\strng{authorfullhash}{5c3bfc4d726caf513f9cf032d754e0a5}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{journaltitle}{arXiv preprint arXiv:1401.7087}
\field{title}{How" quantum" is the D-Wave machine?}
\field{year}{2014}
\endentry
\entry{Albash2020ComparingRM}{article}{}
\name{author}{2}{}{%
\list{publisher}{1}{%
\strng{namehash}{90cc5ae30836566a8e362c5f25fe5c3a}
\strng{fullhash}{90cc5ae30836566a8e362c5f25fe5c3a}
\strng{bibnamehash}{90cc5ae30836566a8e362c5f25fe5c3a}
\strng{authorbibnamehash}{90cc5ae30836566a8e362c5f25fe5c3a}
\strng{authornamehash}{90cc5ae30836566a8e362c5f25fe5c3a}
\strng{authorfullhash}{90cc5ae30836566a8e362c5f25fe5c3a}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{journaltitle}{Physical Review Applied}
\field{number}{1}
\field{title}{Comparing relaxation mechanisms in quantum and classical transverse-field annealing}
\field{volume}{15}
\field{year}{2021}
\field{pages}{014029}
\range{pages}{1}
\endentry
\entry{amin:freeze-out}{article}{}
\name{author}{1}{}{%
given={Mohammad\bibnamedelima H},
giveni={M\bibinitperiod\bibinitdelim H\bibinitperiod}}}%
\list{publisher}{1}{%
\strng{namehash}{66344a2770a9da8b5f33733719d53940}
\strng{fullhash}{66344a2770a9da8b5f33733719d53940}
\strng{bibnamehash}{66344a2770a9da8b5f33733719d53940}
\strng{authorbibnamehash}{66344a2770a9da8b5f33733719d53940}
\strng{authornamehash}{66344a2770a9da8b5f33733719d53940}
\strng{authorfullhash}{66344a2770a9da8b5f33733719d53940}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{journaltitle}{Physical Review A}
\field{number}{5}
\field{title}{Searching for quantum speedup in quasistatic quantum annealers}
\field{volume}{92}
\field{year}{2015}
\field{pages}{052323}
\range{pages}{1}
\endentry
\entry{spin-bath-polarization}{misc}{}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labeltitlesource}{title}
\field{title}{{D-Wave System Documentation -- Spin-Bath Polarization Effect}}
\field{year}{2021}
\verb{urlraw}
\verb https://docs.dwavesys.com/docs/latest/c_qpu_2.html#spin-bath-polarization-effect
\endverb
\verb{url}
\verb https://docs.dwavesys.com/docs/latest/c_qpu_2.html#spin-bath-polarization-effect
\endverb
\endentry
\entry{chen2020pausing}{article}{}
\name{author}{2}{}{%
given={Daniel\bibnamedelima A.},
giveni={D\bibinitperiod\bibinitdelim A\bibinitperiod}}}%
\list{publisher}{1}{%
{American Physical Society}%
\strng{namehash}{aa56b96f8a58f6592a2e2a4385d9623b}
\strng{fullhash}{aa56b96f8a58f6592a2e2a4385d9623b}
\strng{bibnamehash}{aa56b96f8a58f6592a2e2a4385d9623b}
\strng{authorbibnamehash}{aa56b96f8a58f6592a2e2a4385d9623b}
\strng{authornamehash}{aa56b96f8a58f6592a2e2a4385d9623b}
\strng{authorfullhash}{aa56b96f8a58f6592a2e2a4385d9623b}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{journaltitle}{Physical Review Applied}
\field{month}{07}
\field{number}{1}
\field{title}{Why and When Pausing is Beneficial in Quantum Annealing}
\field{volume}{14}
\field{year}{2020}
\field{pages}{014100\bibrangedash}
\range{pages}{-1}
\verb{doi}
\verb 10.1103/PhysRevApplied.14.014100
\endverb
\verb{urlraw}
\verb https://link.aps.org/doi/10.1103/PhysRevApplied.14.014100
\endverb
\verb{url}
\verb https://link.aps.org/doi/10.1103/PhysRevApplied.14.014100
\endverb
\endentry
\entry{LidarWhaley:03}{inbook}{}
\name{author}{2}{}{%
given={D.\bibnamedelimi A.},
giveni={D\bibinitperiod\bibinitdelim A\bibinitperiod}}}%
given={K.\bibnamedelimi B.},
giveni={K\bibinitperiod\bibinitdelim B\bibinitperiod}}}%
\list{location}{1}{%
\list{publisher}{1}{%
\strng{namehash}{abd5ee69f1e09450c7c7d84213628026}
\strng{fullhash}{abd5ee69f1e09450c7c7d84213628026}
\strng{bibnamehash}{abd5ee69f1e09450c7c7d84213628026}
\strng{authorbibnamehash}{abd5ee69f1e09450c7c7d84213628026}
\strng{authornamehash}{abd5ee69f1e09450c7c7d84213628026}
\strng{authorfullhash}{abd5ee69f1e09450c7c7d84213628026}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{booktitle}{Irreversible Quantum Dynamics}
\field{series}{Lecture Notes in Physics}
\field{volume}{622}
\field{year}{2003}
\field{pages}{83}
\range{pages}{1}
\verb{urlraw}
\verb https://www.springer.com/gp/book/9783540402237
\endverb
\verb{url}
\verb https://www.springer.com/gp/book/9783540402237
\endverb
\endentry
\entry{kubo1957statistical}{article}{}
\name{author}{1}{}{%
\list{publisher}{1}{%
{The Physical Society of Japan}%
\strng{namehash}{3db80a54dbc22f6e0f8ef62d81d5e4c4}
\strng{fullhash}{3db80a54dbc22f6e0f8ef62d81d5e4c4}
\strng{bibnamehash}{3db80a54dbc22f6e0f8ef62d81d5e4c4}
\strng{authorbibnamehash}{3db80a54dbc22f6e0f8ef62d81d5e4c4}
\strng{authornamehash}{3db80a54dbc22f6e0f8ef62d81d5e4c4}
\strng{authorfullhash}{3db80a54dbc22f6e0f8ef62d81d5e4c4}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{journaltitle}{J. Phys. Soc. Jpn.}
\field{number}{6}
\field{title}{Statistical-mechanical theory of irreversible processes. I. General theory and simple applications to magnetic and conduction problems}
\field{volume}{12}
\field{year}{1957}
\field{pages}{570\bibrangedash 586}
\range{pages}{17}
\verb{doi}
\verb 10.1143/JPSJ.12.570
\endverb
\endentry
\entry{martin1959theory}{article}{}
\name{author}{2}{}{%
given={Paul\bibnamedelima C},
giveni={P\bibinitperiod\bibinitdelim C\bibinitperiod}}}%
\list{publisher}{1}{%
\strng{namehash}{a4333f1686e2d7af55bd6904099924f3}
\strng{fullhash}{a4333f1686e2d7af55bd6904099924f3}
\strng{bibnamehash}{a4333f1686e2d7af55bd6904099924f3}
\strng{authorbibnamehash}{a4333f1686e2d7af55bd6904099924f3}
\strng{authornamehash}{a4333f1686e2d7af55bd6904099924f3}
\strng{authorfullhash}{a4333f1686e2d7af55bd6904099924f3}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{journaltitle}{Phys. Rev.}
\field{number}{6}
\field{title}{Theory of many-particle systems. I}
\field{volume}{115}
\field{year}{1959}
\field{pages}{1342}
\range{pages}{1}
\verb{doi}
\verb 10.1103/PhysRev.115.1342
\endverb
\endentry
\entry{albash2015consistency}{article}{}
\name{author}{5}{}{%
given={Paul\bibnamedelima A.},
giveni={P\bibinitperiod\bibinitdelim A\bibinitperiod}}}%
given={Daniel\bibnamedelima A.},
giveni={D\bibinitperiod\bibinitdelim A\bibinitperiod}}}%
\list{publisher}{1}{%
{American Physical Society}%
\strng{namehash}{43b9b3f41dd1b283dcb4acfbf0102764}
\strng{fullhash}{43b9b3f41dd1b283dcb4acfbf0102764}
\strng{bibnamehash}{43b9b3f41dd1b283dcb4acfbf0102764}
\strng{authorbibnamehash}{43b9b3f41dd1b283dcb4acfbf0102764}
\strng{authornamehash}{43b9b3f41dd1b283dcb4acfbf0102764}
\strng{authorfullhash}{43b9b3f41dd1b283dcb4acfbf0102764}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{issue}{4}
\field{journaltitle}{Phys. Rev. A}
\field{month}{4}
\field{title}{Consistency tests of classical and quantum models for a quantum annealer}
\field{volume}{91}
\field{year}{2015}
\field{pages}{042314}
\range{pages}{1}
\verb{doi}
\verb 10.1103/PhysRevA.91.042314
\endverb
\verb{urlraw}
\verb https://link.aps.org/doi/10.1103/PhysRevA.91.042314
\endverb
\verb{url}
\verb https://link.aps.org/doi/10.1103/PhysRevA.91.042314
\endverb
\endentry
\entry{ronnow:speedup}{article}{}
\name{author}{9}{}{%
given={Troels\bibnamedelima F.},
giveni={T\bibinitperiod\bibinitdelim F\bibinitperiod}}}%
given={Sergei\bibnamedelima V.},
giveni={S\bibinitperiod\bibinitdelim V\bibinitperiod}}}%
given={John\bibnamedelima M.},
giveni={J\bibinitperiod\bibinitdelim M\bibinitperiod}}}%
given={Daniel\bibnamedelima A.},
giveni={D\bibinitperiod\bibinitdelim A\bibinitperiod}}}%
\list{publisher}{1}{%
{American Association for the Advancement of Science}%
\strng{namehash}{0acc9ae807bdef2ecf429a5a8ab9abea}
\strng{fullhash}{85e4f918b798b4fb4d14d5047b8f759c}
\strng{bibnamehash}{0acc9ae807bdef2ecf429a5a8ab9abea}
\strng{authorbibnamehash}{0acc9ae807bdef2ecf429a5a8ab9abea}
\strng{authornamehash}{0acc9ae807bdef2ecf429a5a8ab9abea}
\strng{authorfullhash}{85e4f918b798b4fb4d14d5047b8f759c}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{abstract}{Quantum machines offer the possibility of performing certain computations much faster than their classical counterparts. However, how to define and measure quantum speedup is a topic of debate. R{ø}nnow et al. describe methods for fairly evaluating the difference in computational power between classical and quantum processors. They define various types of quantum speedup and consider quantum processors that are designed to solve a specific class of problems.Science, this issue p. 420The development of small-scale quantum devices raises the question of how to fairly assess and detect quantum speedup. Here, we show how to define and measure quantum speedup and how to avoid pitfalls that might mask or fake such a speedup. We illustrate our discussion with data from tests run on a D-Wave Two device with up to 503 qubits. By using random spin glass instances as a benchmark, we found no evidence of quantum speedup when the entire data set is considered and obtained inconclusive results when comparing subsets of instances on an instance-by-instance basis. Our results do not rule out the possibility of speedup for other classes of problems and illustrate the subtle nature of the quantum speedup question.}
\field{issn}{0036-8075}
\field{journaltitle}{Science}
\field{number}{6195}
\field{title}{Defining and detecting quantum speedup}
\field{volume}{345}
\field{year}{2014}
\field{pages}{420\bibrangedash 424}
\range{pages}{5}
\verb{doi}
\verb 10.1126/science.1252319
\endverb
\verb{eprint}
\verb http://science.sciencemag.org/content/345/6195/420.full.pdf
\endverb
\verb{urlraw}
\verb http://science.sciencemag.org/content/345/6195/420
\endverb
\verb{url}
\verb http://science.sciencemag.org/content/345/6195/420
\endverb
\endentry
\entry{vestgaarden2008nonlinearly}{article}{}
\name{author}{3}{}{%
\list{publisher}{1}{%
\strng{namehash}{b057237d55421bf4de7630fc8db3d7a6}
\strng{fullhash}{b057237d55421bf4de7630fc8db3d7a6}
\strng{bibnamehash}{b057237d55421bf4de7630fc8db3d7a6}
\strng{authorbibnamehash}{b057237d55421bf4de7630fc8db3d7a6}
\strng{authornamehash}{b057237d55421bf4de7630fc8db3d7a6}
\strng{authorfullhash}{b057237d55421bf4de7630fc8db3d7a6}
\field{sortinit}{1}
\field{sortinithash}{50c6687d7fc80f50136d75228e3c59ba}
\field{labelnamesource}{author}
\field{labeltitlesource}{title}
\field{journaltitle}{Physical Review B}
\field{number}{1}
\field{title}{Nonlinearly driven Landau-Zener transition in a qubit with telegraph noise}
\field{volume}{77}
\field{year}{2008}
\field{pages}{014514}
\range{pages}{1}
\endentry
\enddatalist
\endrefsection
\endinput

\end{singlespace}

% Appendices
\phantomsection
\addcontentsline{toc}{chapter}{Appendices}%
\markboth{Appendices}{Appendices}%
\chapter*{Appendices}
\renewcommand\thesection{\Alph{section}}
\renewcommand*{\thesubsection}{\Alph{section}.\arabic{subsection}}
\begingroup
\numberwithin{equation}{section}
% Appendix source files
% \section{A Long Proof}
% \label{app:long_proof}

% \subsection{Part one}
% A long proof that nobody is gonna read goes here.

% \subsection{Part two}
% Part two of the long proof that nobody is gonna read goes here.
\section{Quantum trajectories for time-dependent adiabatic master equations}
\subsection{Error estimates}
\label{app:A}
Let us assume that the system is in a pure state at time $t$, i.e., $ρ_S(t) = |ψ(t)⟩⟨$|, and let us consider a single time-step. In a single trajectory, the evolution of $\ket{\psi(t)}$ involves two possibilities: no jump or a jump. The ensemble average of trajectories after one finite time-step mainly involves two kinds of errors: the error associated with the norm of the state vector (Sec. <ref>), and the error associated with the probability elements in a finite time-step (Sec. <ref>).

§.§.§ Error associated with the norm squared of a no-jump trajectory

The Schrödinger equation of the effective Hamiltonian in the case of a no-jump trajectory is given by

\begin{equation}
\label{eqt:Schrapp2}
\frac{d\ket{\psi(t)}}{dt} = -i H_{\text{eff}}(t) \ket{\psi(t)} \ .
\end{equation}

The resulting state after one time-step $\Delta t$ is:

\begin{equation}
\label{eqt:euler}
\ket{\tilde{\psi}(t+\Delta t)} =V_{\text{eff}}(t+\Delta t, t)\ket{\psi(t)}  \ ,
\end{equation}

V_eff(t+Δt, t) = Texp[-i∫_t^t+ΔtH_eff(t')dt']
is a non-unitary contractive evolution operator ($\vertiii{ V_{\text{eff}}(t+\Delta t, t)} \leq 1$ for every induced norm) and $\ket{\tilde{\psi}(t+\Delta t)}$ is unnormalized since $H_{\text{eff}}(t)$ is not hermitian.

We can expand the time-ordered exponential for $V_{\text{eff}}(t+\Delta t, t)$ as

\begin{align}
\label{Eqt: TDtaylor}
&{}{\cal T}\exp\left[-i\int_t^{t+\Delta t}H_{\text{eff}}(t')dt'\right]    \\
&=\mathds{1} - i \int_{t}^{t+\Delta t} dt' H_{\text{eff}}(t')   + \left(-i\right)^2 \frac{1}{2!} \mathcal{T}\left(\int_{t}^{t+\Delta t} dt' H_{\text{eff}}(t')\right)^2 +
\left(-i\right)^3 \frac{1}{3!}\mathcal{T}\left(\left(\int_{t}^{t+\Delta t} dt' H_{\text{eff}}(t')\right)^3\right) \cdots \nonumber
\end{align}

For any $H_{\text{eff}}(t)$ that is $C^{K}$ on an open interval containing $[0, t_f]$ ($H^{(k)}_{\text{eff}}(t)$ is continuous and bounded for $1\leq k \leq K$):

\begin{equation}
\int_{t}^{t+\Delta t} dt' H_{\text{eff}}(t') =
\sum\limits_{k = 0}^{K-1}\frac{H^{(k)}_{\text{eff}}(t)}{(k+1)!}(\Delta t)^{k+1} + O((\Delta t)^{K+1})\, ,
\end{equation}

as $\Delta t \rightarrow 0$. Assume $H_{\text{eff}}(t)$ is $C^{K}$ and $K \geq 2$. Let $\dot{H}_{\eff}(t)\equiv H^{(1)}_{\text{eff}}(t)$ denote the time derivative of $H_\eff(t)$.
We have:

\begin{align}
\label{Eqt: integral}
\int_{t}^{t+\Delta t} dt' H_{\text{eff}}(t')&= H_{\text{eff}}(t)\Delta t + \frac{1}{2!}\dot{H}_{\text{eff}}(t) (\Delta t)^2 + O((\Delta t)^3) \,,\\
\mathcal{T}\left(\int_{t}^{t+\Delta t} dt' H_{\text{eff}}(t')\right)^2 &= 2\int_{t}^{t+\Delta t} dt' H_{\text{eff}}(t')  \int_{t}^{t'} dt'' H_{\text{eff}}(t'') \nonumber\\
&= H_{\text{eff}}(t)^2(\Delta t)^2 + \frac{1}{2}H_{\text{eff}}(t)\dot{H}_{\text{eff}}(t)(\Delta t)^3 + \frac{1}{2}\dot{H}_{\text{eff}}(t)H_{\text{eff}}(t)(\Delta t)^3 + O((\Delta t)^4)\ \,, \\
\mathcal{T}\left(\int_{t}^{t+\Delta t} dt' H_{\text{eff}}(t')\right)^3 &= H_{\text{eff}}(t)^3(\Delta t)^3 + O((\Delta t)^4) \,.
\end{align}

Therefore, we have, to second order in $\Delta t$:

V_eff(t+Δt, t) = 1 -i ( H_eff(t)Δt + 1/2!Ḣ_eff(t) (Δt)^2 ) + (-i)^2 1/2! H_eff(t)^2(Δt)^2  + O((Δt)^3) ,

so that:

\begin{align}
&{} V^{\dagger}_{\text{eff}}(t+\Delta t, t) V_{\text{eff}}(t+\Delta t, t) \nonumber \\
=&{}  \mathds{1} + i\Delta t \left(H^{\dagger}_{\text{eff}}(t) - H_{\text{eff}}(t)\right)
+ \frac{1}{2!}(\Delta t)^2\left[ i\dot{H}^{\dagger}_{\text{eff}}(t) - i\dot{H}_{\text{eff}}(t)
- \left(H^{\dagger 2}_{\text{eff}}(t) + H^2_{\text{eff}}(t)\right)\right] + O((\Delta t)^3) \ .
\end{align}

We approximate $N(t,\Delta t)\equiv \|\ket{\tilde{\psi}(t+\Delta t)} \|^2$ by:

N(t,Δt) = ⟨ψ̃(t+Δt)|ψ̃(t+Δt)|≈⟩⟨ψ(t)|(1 + i Δt(H^†_eff(t) - H_eff(t) ))|ψ(t)⟩ ≡E(t,Δt)  .

The approximation error is given by:

\begin{align}
\delta (t,\Delta t) &\equiv N(t,\Delta t)-E(t,\Delta t)\\
&= \frac{1}{2}(\Delta t)^2\bra{\psi(t)}\left[ i\dot{H}^{\dagger}_{\text{eff}}(t) - i\dot{H}_{\text{eff}}(t)
- \left(H^{\dagger 2}_{\text{eff}}(t) + H^2_{\text{eff}}(t)\right)\right] \ket{\psi(t)}+ O((\Delta t)^3) \\
&\leq \frac{(\Delta t)^2}{2}\vertiii{i(\dot{H}^{\dagger}_{\text{eff}}(t) - \dot{H}_{\text{eff}}(t)) - ( H^{\dagger 2}_{\text{eff}}(t) + H^2_{\text{eff}}(t))} + \vertiii{O((\Delta t)^3)}\ ,
\label{Eqt: errornorm}
\end{align}

where $\|\cdot\|$ is the operator norm (largest singular value).

Eq. (<ref>) gives the relation between the error of the norm square approximation and the time-step. It is also helpful to consider the sources of error here as they will be used later. This error mainly arises from two sources during the truncations of Taylor expansion:

* The truncation of the Taylor expansion of the integral [Eq. (<ref>)] to keep only $H_{\text{eff}}(t)\Delta t$. This turns Eq. (<ref>) into $\exp\left(-iH_{\text{eff}}(t)\Delta t\right)$.  For this to hold, we require

\begin{align}
\label{Eqt: case1}
\left\lVert \frac{1}{2}\dot{H}_{\text{eff}}(t) (\Delta t)^2
+ O((\Delta t)^3) \right\rVert \ll  \lVert H_{\text{eff}}(t)\Delta t \rVert \, ,
\end{align}

\begin{align}
\Delta t \ll 2\frac{\vertiii{H_{\text{eff}}(t)}}{\vertiii{\dot{H}_{\text{eff}}(t)}}  \, ,
\label{Eqt: 1stcon}
\end{align}

assuming $\dot{H}_{\text{eff}}(t)\neq 0$; if it is then the condition becomes $\Delta t \ll \left(\frac{K!\vertiii{H_{\text{eff}}(t)}}{\vertiii{{H^{(K)}_{\text{eff}}(t)}}}\right)^{1/K}$ for the lowest value of $K$ such that $H^{(K)}_{\text{eff}}(t)\neq 0$, where the superscript denotes the $K$-th derivative.

* Keeping only the first order term in Eq. (<ref>) afterwards, i.e., $\mathds{1} - iH_{\text{eff}}(t)$.
This requires, in addition to Eq. (<ref>):

\begin{align}
\label{Eqt: case2}
\left\lVert\frac{1}{2}H^2_{\text{eff}}(t)(\Delta t)^2 + O((\Delta t)^3)
\right\rVert \ll  \vertiii{H_{\text{eff}}(t)\Delta t} \, ,
\end{align}

implying, for all $t$ such that the denominators do not vanish:

Δt ≪min{ H_eff(t)/Ḣ_eff(t), 1/H_eff(t) }  ,
where we ignored the factor of $2$.

In conclusion, the norm after one time-step is related to the jump probabilities as:

\begin{align}
N(t,\Delta t) &=
\bra{\psi(t)}\left(\mathds{1} + i (\Delta t)\left(H^{\dagger}_{\text{eff}}(t) - H_{\text{eff}} \right)\right)\ket{\psi(t)} \notag \\
&\qquad + \delta(t,\Delta t) \\
&= 1 - \Delta t \sum\limits_{i}\braket{{A}_{i}^{\dagger}(t) {A}_{i} (t)} + \delta(t,\Delta t) \\
&= 1 - \Delta p(t) + \delta(t,\Delta t) \,,\label{Eqt: squareofnorm}
\end{align}

where we used Eq. (<ref>) and defined the (approximate) jump probabilities as:

\begin{align}
&{}\Delta p(t)  = \sum_i \Delta p_i(t) \\
&{}\Delta p_i(t) = \Delta t \braket{ {A}_{i}^{\dagger}(t) {A}_{i} (t)}  \, .
\label{eqt:typeofprob}
\end{align}

This explains Eqs. (<ref>) and (<ref>).

§.§.§ Error associated with probability elements

We have defined the $\Delta p_i(t)$ and $\Delta p(t)$ in a fixed time-step in Eq. (<ref>).
Note that even in the time-independent case, where both the Hamiltonian and Lindblad operators are time-independent, $1- \Delta p(t)$ is never exactly equal to the norm squared of the state vector after one time-step, and $\Delta p_i(t)$ is not exactly the jump probability inside the time-step. They are only approximations.

For a finite time-step $\Delta t$, with $p_0$ the probability of having no jump inside the interval $[t, t+\Delta t]$ and $p_1$ the probability of having one jump inside the interval $[t, t+\Delta t]$, we have:
p_0 = e^-∫_t^t+Δt λ(t') dt' = 1-p_1   .

Note that as $\Delta t \rightarrow 0$, $p_{0} = 1 - \lambda(t)\Delta t + o(\Delta t)$, $p_{1} \simeq \lambda(t)\Delta t$.
We shall focus on the case where the time-step is sufficiently small such that the probability of two or more jumps occurring within a single time-step is negligible.[E.g., the quantum jump can be described by a Poisson process with a state-dependent inhomogeneous jump rate, with two or more jumps as successive one and no-jump processes, and $\sum_{n\geq 2}p_{n} = o(\Delta t)$ as $\Delta t\rightarrow 0$.]

First, we can expand the exponential as

\begin{align}
&=1 - \int_t^{t+\Delta t}
\lambda(t')dt' + \frac{1}{2}\left(\int_t^{t+\Delta t} \lambda(t') dt'
\right)^2 \cdots \nonumber
\end{align}

During the time window from $t$ to $t+\Delta t$, $\lvert\psi(t')\rangle$ [required to calculate $\lambda(t')$] is obtained by solving the Schrödinger equation with the effective Hamiltonians [Eq. (<ref>)] and renormalizing the solution during the integration. As we show in Sec. <ref>, in any finite interval the norm squared of the unnormalized state vector [Eq. (<ref>)] is equal to $p_0$. This is the reason why we can use the waiting time distribution (or Gillespie algorithm <cit.>) as our simulation method in the main text.

Second, since

\begin{align}
\int_{t}^{t+\Delta t} dt' \lambda(t') = \sum\limits_{k = 0}^{\infty}\frac{\lambda^{(k)}(t)}{(k+1)!}(\Delta t)^{k+1} \ ,
\end{align}

we have

\begin{align}
\label{eqt:p0allterms}
= 1 - \lambda(t)\Delta t + e_p \ ,
\end{align}

where the error $e_p$ associated with the probability elements in a fixed time-step is:

\begin{equation}
e_p = - \frac{1}{2}\frac{d}{dt}
\lambda(t)(\Delta t)^2 + \frac{1}{2}
\lambda^2(t)(\Delta t)^2 + O((\Delta t)^3) \,.
\end{equation}

This should be much smaller than the first order term $\lambda(t) \Delta t$. Therefore, we need

\begin{align}
\Delta t \ll \left|\frac{ \lambda(t) }{  \lambda^2(t)-\frac{d}{dt} \lambda(t)}\right|  \, .
\label{Eqt: Poissoncon}
\end{align}

In the time-independent case, this reduces to
Δt ≪1/ λ  .

§.§ Proof of equivalence between the master equation and trajectories formulations

Our goal in this section is to show how the master equation, Eq. (<ref>), can be recovered from the quantum trajectories formulation, and to find a bound on the time-step $\Delta t$. This generalizes the proof for the time-independent case found in <cit.>.

§.§.§ To jump or not to jump

The probability elements $\Delta p(t)$ and $\Delta p_i (t)$ are important for determining whether a jump occurs and if a jump does occur, which jump type occurs.
In order to determine if a jump occurs or not, we draw a random number $\epsilon$, uniformly distributed between 0 and 1. If $\Delta p(t) < \epsilon$, which is almost always the case since $\Delta p(t)$ is very small, no jump occurs. In the case of no jump, $\ket{\psi(t)}$  evolves according to the effective Schrödinger equation, Eq. (<ref>). At time $t+\Delta t$ we simply renormalize the solution of Eq. (<ref>):

\begin{align}
\left|\psi(t + \Delta t)\right> &= \frac{1}{\sqrt{\braket{\tilde{\psi}(t+\Delta t)| \tilde{\psi}(t+\Delta t)}}}\ket{\tilde{\psi}(t+\Delta t)}\\&\stackrel{\eqref{Eqt: squareofnorm}}{=} \frac{1}{\sqrt{1 - \Delta p(t) + \delta}}\ket{\tilde{\psi}(t+\Delta t)} \ .
\end{align}

If $\Delta p(t) > \epsilon$, the state undergoes an abrupt jump and we choose the new wavefunction among the different states $A_{i}\ket{\psi(t)}$ and renormalize:

\begin{align}
\left|\psi(t + \Delta t)\right> &= \frac{{A}_{i}(t)\ket{\psi(t)}}{\sqrt{\bra{\psi(t)} {A}_{i}^{\dagger}(t) {A}_{i} (t)\ket{\psi(t)}}} \\
\sqrt{\frac{\Delta t}{\Delta p_i(t)}} {A}_{i}(t)\ket{\psi(t)}\ .
\end{align}

Which type of jumps occurs is determined according to the probability

\begin{align}
\Pi_i(t) &= \frac{\Delta p_i(t)}{\Delta p(t)} = \frac{\braket{\psi(t)|{A}_{i}^{\dagger}(t){A}_{i}(t)|\psi(t)}\Delta t}{\sum\limits_{i}\braket{\psi(t)|{A}_{i}^{\dagger}(t){A}_{i}(t)|\psi(t)}\Delta t} \\
&= \frac{\braket{{A}_{i}^{\dagger}(t){A}_{i}(t)}}{\lambda(t)}\ ,
\label{eq:Pi(t)}
\end{align}

λ(t) = ∑_i⟨A_i^†(t)A_i(t)⟩
is the time-dependent jump rate.

§.§.§ Averaging over trajectories

Let $H_{\text{eff}}(t)$ be $C^{K}$ with $K \geq 2$.
We first express the mean value $\bar{\sigma}_{S}(t)$ as a sum over the non-Hermitian evolution [with probability $1 - \Delta p(t)$] and the jump trajectories [with probability $\Delta p(t)$], so that as $\Delta t \rightarrow 0$ we have:

\begin{align}
\label{Eqt: averaging}
\bar{\sigma}_{S}(t + \Delta t)
&= \left[1 - \Delta p(t)\right]\frac{\ket{\tilde{\psi}(t + \Delta t)}}{\sqrt{1 - \Delta p(t) + \delta}}\frac{\bra{\tilde{\psi}(t + \Delta t)}}{\sqrt{1 - \Delta p(t) + \delta}}  \nonumber \\
& + \Delta p(t) \sum\limits_{i} \Pi_{i}(t) \sqrt{\frac{\Delta t}{\Delta p_i(t)}}{A}_{i}(t)  \ket{\psi(t)}\sqrt{\frac{\Delta t}{\Delta p_i(t)}}\bra{\psi(t)}{A}_{i}^{\dagger}(t) \\
&\stackrel{\eqref{eq:Pi(t)}}{=}   \frac{1 - \Delta p(t)}{1 - \Delta p(t) + \delta}\ket{\tilde{\psi}(t + \Delta t)}\bra{\tilde{\psi}(t + \Delta t)}
+ \Delta t \sum\limits_{i}{A}_{i}(t)\sigma_{S}(t){A}_{i}^{\dagger}(t)
\label{Eqt: noninnerproductform}
\end{align}

Combining Eq. (<ref>) with Eq. (<ref>) we have:
|ψ̃(t+Δt)⟩ =
[ 1 -i ( H_eff(t)Δt + 1/2!Ḣ_eff(t) (Δt)^2 ) + (-i)^2 1/2! H_eff(t)^2(Δt)^2  + O((Δt)^3) ] |ψ(t)⟩   .
Recall that in Eq. (<ref>) we gave conditions allowing us to neglect the $O((\Delta t)^2)$ terms.

\begin{align}
\label{eq:B3}
\bar{\sigma}_{S}(t + \Delta t)
&=  \frac{1 - \Delta p(t)}{1 - \Delta p(t) + \delta}\left( \mathds{1} - i H_{\text{eff}}(t) \Delta t
+ O((\Delta t)^2)\right)\bar{\sigma}_{S}(t) \left(\mathds{1} + i H^{\dagger}_{\text{eff}}(t) \Delta t
+O((\Delta t)^2)\right)\\
&\phantom{==}+ \Delta t \sum\limits_{i}{A}_{i}(t)\bar{\sigma}_{S}(t){A}_{i}^{\dagger}(t)\nonumber  \, .
\end{align}

where we have replaced $\sigma_{S} (t)$ by $\bar{\sigma}_{S} (t)$ after averaging over many trajectories.
Rearranging this expression into a form that exposes the terms that will become the master equation,
the expression for the averaged state at $t + \Delta t$ becomes:

\begin{align}
\bar{\sigma}_{S}(t + \Delta t) &= \bar{\sigma}_{S}(t) + i\Delta t\left(\bar{\sigma}_{S}(t)H_{\text{eff}}^{\dagger}(t) - H_{\text{eff}}(t)\bar{\sigma}_{S}(t) \right)+ \Delta t \sum_{i} {A}_{i}(t) \bar{\sigma}_{S}(t){A}_{i}^{\dagger}(t) \\
- \frac{\delta}{1 - \Delta p(t) + \delta} \left[ \bar{\sigma}_{S}(t) +  i\bar{\sigma}_{S}(t)\left(H_{\text{eff}}^{\dagger}(t) -H_{\text{eff}}(t)\right) \bar{\sigma}_{S}(t) \Delta t \right] + O((\Delta t)^2)\ .
\label{eq:B5b}
\end{align}

Note that $\delta$ is $O((\Delta t)^2)$ as $\Delta t \rightarrow 0$ [Eq. (<ref>)], and $\Delta p(t) = \Delta t \sum_{i}\braket{\psi(t)|{A}_{i}^{\dagger}(t){A}_{i}(t)|\psi(t)}$ is $O(Δt)$ as $Δt →0$, so that:
\bes
\begin{align}
\frac{\delta}{1 - \Delta p(t) + \delta} &=
\frac{O((\Delta t)^2)}{1 - O(\Delta t) + O((\Delta t)^2)} \\
&=O((\Delta t)^2)\, .
\end{align}
\ees
Therefore line~\eqref{eq:B5b} can be absorbed into $O((Δt)^2)$, and we are left with:
\begin{align}
\frac{\bar{\sigma}_{S}(t + \Delta t) - \bar{\sigma}_{S}(t)}{\Delta t} =&
-i\left(H_{\text{eff}}(t)\bar{\sigma}_{S}(t) - \bar{\sigma}_{S}(t)H_{\text{eff}}^{\dagger}(t) \right)\notag \\
&+  \sum_{i} {A}_{i}(t) \bar{\sigma}_{S}(t){A}_{i}^{\dagger}(t)
+ O(\Delta t) \ ,
\end{align}
which becomes the master equation, Eq.~\eqref{eqt:effdiagonalform}, in the $Δt →0$ limit.

\subsubsection{Upper bound on $\Delta t$}
\label{ssec:upperbound}
The above proof takes $Δt →0$. We would like to know how small the time-step $Δt$ should be in order for the approximations made to be valid. In Eq.~\eqref{eq:B3}, we expanded the time-ordered exponential, and kept only the first order terms. This is equivalent to the criteria in Eqs.~\eqref{Eqt: case1} and ~\eqref{Eqt: case2}, summarized as a single condition in Eq.~\eqref{eqt:Deltat2}.
As shown in Sec.~\ref{sssec:en}, this also automatically makes the error in the norm squared approximation $δ$ small. We also need to satisfy Eq.~\eqref{Eqt: Poissoncon}, in order to accurately approximate the probability elements. Taken together, therefore:
\beq
\Delta t \ll  \min \left\{ \frac{\vertiii{H_{\text{eff}}(t)}}{\vertiii{\dot{H}_{\text{eff}}(t)}}, \frac{1}{\vertiii{H_{\text{eff}}(t)}}, \left|\frac{ \lambda(t) }{  \lambda^2(t)-\dot \lambda(t)}\right| \right\} \ .
\label{eq:B7}
\eeq
In practice, choosing a constant time-step that satisfies Eq.~\eqref{eq:B7} in the whole timespan $[0, t_f]$ is sufficient, though one might prefer to implement an adaptive time-step tailored to the instantaneous value of the R.H.S.
\subsection{On the validity of waiting times (quantum time-dependent operators)}
\label{sec:C}
Here we show the validity of using the waiting time distribution in the case of time-dependent operators. The argument presented here is based on Ref.~\cite{Breuer:2002} and we extend it to the time-dependent case.

Let us denote by $|ψ(t)⟩$ and $|ψ̃(t)⟩$ the normalized and unnormalized state vectors respectively, and let us assume they are equal at time $t$. This can happen when $t = 0$ or any time immediately after each jump. Let $t^+ ≡t+τ$, where $τ$ can be as large as is possible until the next jump occurs,  and
\begin{equation}
V_{\text{eff}}(t^+, t) = {\cal T}\exp\left[-i\int_t^{t^+}H_{\text{eff}}(t')dt'\right] \ ,
\end{equation}
\bes
\begin{align}
\ket{\tilde{\psi}(t + \tau)} &=
V_{\text{eff}}(t^+, t)\ket{\psi(t)}  \ ,\\
\ket{\psi(t^+)} &=
\frac{V_{\text{eff}}(t^+, t)\ket{\psi(t)}}{\left\lVert V_{\text{eff}}(t^+, t)\ket{\psi(t)} \right\rVert} \,.
\label{eq:C2b}
\end{align}
\ees
Then, starting from $t$, for any future $t^+ > t$, we have
\bes
\begin{align}
\frac{d\phantom{,}}{dt^{+}}\left\lVert \ket{\tilde{\psi}(t^{+})} \right\rVert^2 &= \phantom{}\frac{d\phantom{,}}{dt^{+}}\left\lVert{\cal T}\exp\left[-i\int_t^{t^+}H_{\text{eff}}(t')dt'\right]\ket{\psi(t)}\right\rVert^2 \\
&= \frac{d\phantom{,}}{dt^{+}}\bra{\psi(t)}V_\eff^{\dagger}(t^+,t) V_\eff(t^+,t)\ket{\psi(t)} \\
&=\bra{\psi(t)}V_\eff^{\dagger}(t^+,t)(+i)H^{\dagger}_{\text{eff}}(t^{+})V_\eff^{}(t^+,t)\ket{\psi(t)} +
\bra{\psi(t)}V_\eff^{\dagger}(t^+,t)(-i)H_{\text{eff}}(t^{+})V_\eff^{}(t^+,t)\ket{\psi(t)} \\
&=-\bra{\psi(t)}V_\eff^{\dagger}(t^+,t)\left(\sum\limits_{i}{A}^{\dagger}_{i}(t^{+}){A}_{i}(t^{+})\right)V_\eff^{}(t^+,t)\ket{\psi(t)} \\
&=  -\left\lVert V_{\text{eff}}(t^{+}, t)\ket{\psi(t)} \right\rVert^2 \sum\limits_{i} \bra{\psi(t^{+})} {A}^{\dagger}_{i}(t^{+}){A}_{i}(t^{+}) \ket{\psi(t^{+})}  \, ,
\end{align}
\ees
where in the last equality we used Eq.~\eqref{eq:C2b}.

Let $N(t^+) ≡ |ψ̃(t^+)⟩ ^2$, as in Eq.~\eqref{Eqt: Step2}.  We have
\begin{equation}
\frac{d\phantom{,}}{dt^{+}} N(t^{+}) = -N(t^{+})  \lambda(t^+)   \, ,
\end{equation}
where $λ$ is the time-dependent jump rate [Eq.~\eqref{eq:lambda}].
The solution to this differential equation with the initial condition $N(t) = 1$ is
\begin{equation}
N(t^{+}) = \exp\left(-\int_t^{t^{+}} \lambda(t') dt'\right) = p_0(t^{+})  \, ,
\end{equation}
where $p_0$ is the probability of not having any jump inside the interval [Eq.~\eqref{eq:p_0}], which we have now shown to be equal to the the norm squared of the unnormalized state vector for any finite interval $[t, t+τ]$.

No commutators of operators at different times appears in the derivation. The use of the waiting time distribution is therefore valid for time-dependent operators as long as the correlation matrix is positive.

\subsection{Effect of the annealing schedule on the $8$-qubit chain example}
\label{app:Chain}

We provide the functional form for the (D-Wave 1) annealing schedule functions $A(t)$ and $B(t)$ used in Sec.~\ref{sec:8qubitchain} in Fig.~\ref{fig:5a}.  We compare this non-linear schedule to a linear annealing schedule with an overall energy scale chosen to closely match the energy scale at which the first annealing schedule curves cross.  The dynamics associated with the linear annealing schedule [shown in Fig.~\ref{fig:5b}] is somewhat different than the non-linear one. This can be attributed to the differences between the two schedules. In order to ensure that the two sets of schedules intersect at the same energy scale, this requires the linear schedule to start from a significantly smaller energy scale.  This results in some important effects on the dynamics.  First, this energy scale is not large enough to ensure that the thermal state at $s = 0$ has negligible weight on excited states.  Since we use the ground state as the initial state of the simulation, which is now sufficiently different from the thermal state at $s=0$, the dissipative dynamics causes visible changes in the state immediately, which can be see as both the dip in Fig~\ref{fig:5a} near $s=0$ and the large number of excitations at the beginning of the anneal in Fig.~\ref{fig:5c}.  Second, even after this initial dip, the lower energy scale associated with the linear schedule means that thermal depopulation of the ground state occurs generally sooner relative to the non-linear schedule studied in the main text, although at a lower rate because of the linear form of the schedule.  Furthermore, because the transverse field remains significant compared to the Ising term for longer in the anneal than in the linear case, repopulation of the ground state due to thermal relaxation occurs for a longer period of time.  However, ultimately, the ground state probability at the end of the anneal is not significantly different than that shown in Fig.~\ref{fig:3qubits8qubits}. The jump statistics are shown in Fig.~\ref{fig:5c}, and closely resemble the non-linear case shown in the inset of Fig.~\ref{fig:8qubitssingle}.

\begin{figure*}[htb!] %  figure placement: here, top, bottom, or page
\centering
\subfigure[\ ]{\includegraphics[width=0.32\textwidth]{qt/AnnealingSchedules}\label{fig:5a}}
\subfigure[\ ]{\includegraphics[width=0.32\textwidth]{qt/EnergyPop_8Qubit_tf=10us_T=20mK_LinearSchedule_8GHz_v1}\label{fig:5b}}
\subfigure[\ ]{\includegraphics[width=0.32\textwidth]{qt/BarPlot_8Qubit_tf=10us_T=20mK_LinearSchedule_8GHz}\label{fig:5c}}
\caption{(a): Dotted lines: Linear schedule; Solid lines: D-Wave schedule used in Fig.~\ref{fig:3qubits8qubits}; (b): Same as in Fig.~\ref{fig:3qubits8qubits} (the evolution of the population in the instantaneous ground state for the $8$-qubit problem), but with a linear schedule. Inset: the convergence of the ground state population towards the AME results. (c): A histogram of the net number of jumps to the instantaneous ground state (GS). It shares the same pattern as in the non-linear schedule case (the inset of Fig.~\ref{fig:8qubitssingle}), but the net number of jumps out of the ground state is smaller due to the smaller energy scale in the linear schedule.}
\label{fig:8qubitslinear}
\end{figure*}

\subsection{Proof that eigenstates of $H_S$ are modified only under jumps}
\label{app:eigenstate-proof}

Recall that the Lindblad operators are defined in Eq.~\eqref{eq:Lindblad2} of the main text as:
\begin{align}
L_{\alpha, \omega}(t) &= \sum_{a,b} \delta_{\omega, \eps_b(t) - \eps_a(t)} \bra{\eps_a(t)} A_\alpha \ket{\eps_b(t)}
| \eps_a(t) \rangle  \langle \eps_b(t)| \,.
\end{align}
After inverting Eq.~\eqref{eq:L} we have
\begin{equation}
A_{i,\omega}(t) = \sum_{\alpha}u^{*}_{i,\alpha}(\omega)L_{\alpha, \omega}(t)
\end{equation}
for the new Lindblad operators corresponding to a diagonalized $γ$ matrix of decay rates.

Assume for simplicity that the $γ$ matrix is already diagonal, so that $u$ is just an identity transformation and $i$ is a relabeling of $α$; then
\begin{align}
A_{i,\omega}(t) &= \sum_{a,b} \delta_{\omega, \eps_b(t) - \eps_a(t)} \bra{\eps_a(t)} A_\alpha \ket{\eps_b(t)}
| \eps_a(t) \rangle  \langle \eps_b(t)| \nonumber \\
A^{\dagger}_{i,\omega}(t) &= \sum_{a,b} \delta_{\omega, \eps_b(t) - \eps_a(t)} \bra{\eps_b(t)} A_\alpha \ket{\eps_a(t)}
| \eps_b(t) \rangle  \langle \eps_a(t)| \,.
\end{align}

Consider the drift term
\[
-\frac{1}{2}\sum_{i} \left[ A_{i}^{\dagger}(t)A_{i}(t) - \braket{A_{i}^{\dagger}(t)A_{i}(t)}\right]\ket{\psi(t)}dt\ .
\]
Since $A_i(t)$ comes from the redefinition where the index $i$ includes the Bohr frequencies [Eq.~\eqref{eqt:Aired}]: $√(γ'_i(ω))A_i,ω(t)→A_i(t)$, the drift term becomes the following after we reintroduce the Bohr frequencies:
\begin{equation}
- \frac{1}{2}\sum_{i}\sum_{\omega} \gamma'_i(\omega)\left[ A_{i,\omega}^{\dagger}(t)A_{i,\omega}(t) - \bra{\psi(t)} A_{i,\omega}^{\dagger}(t)A_{i,\omega}(t) \ket{\psi(t)} \right]\ket{\psi(t)}dt \,.
\label{eq:explicit}
\end{equation}

\subsubsection{When $\ket{\psi(t)}$ is an eigenstate of $H_{S}(t)$}
If the $|ψ(t)⟩$ is an eigenstate of $H_S(t)$ (denoted as $|_b(t)⟩$ here),
\begin{align}
A_{i,\omega}(t)\ket{\psi(t)} = A_{i,\omega}(t)\ket{\eps_b(t)} &=  \sum_{a} \delta_{\omega, \eps_b(t) - \eps_a(t)} \bra{\eps_a(t)} A_\alpha \ket{\eps_b(t)}
| \eps_a(t) \rangle \nonumber \\
\bra{\psi(t)}A^{\dagger}_{i,\omega}(t) = \bra{\eps_b(t)}A^{\dagger}_{i,\omega}(t) &= \sum_{a} \delta_{\omega, \eps_b(t) - \eps_a(t)} \bra{\eps_b(t)} A_\alpha \ket{\eps_a(t)} \langle \eps_a(t)| \,.
\end{align}

Concentrate on the term of the parenthesis inside each summation of Eq.~\eqref{eq:explicit}, i.e.
\begin{equation}
\left[ A_{i,\omega}^{\dagger}(t)A_{i,\omega}(t) - \bra{\eps_b(t)} A_{i,\omega}^{\dagger}(t)A_{i,\omega}(t) \ket{\eps_b(t)}  \right]\ket{\eps_b(t)} \,.
\label{eq:paraterm}
\end{equation}
The first term is:
\begin{align}
&\phantom{{}={}} A_{i,\omega}^{\dagger}(t) A_{i,\omega}(t)\ket{\eps_b(t)}  \nonumber\\
&= \left(\sum_{a^{\prime},b^{\prime}} \delta_{\omega, \eps_{b^{\prime}}(t) - \eps_{a^{\prime}}(t)} \bra{\eps_{b^{\prime}}(t)} A_\alpha \ket{\eps_{a^{\prime}}(t)}
| \eps_{b^{\prime}}(t) \rangle  \langle \eps_{a^{\prime}}(t)| \right) \sum_{a} \delta_{\omega, \eps_b(t) - \eps_a(t)} \bra{\eps_a(t)} A_\alpha \ket{\eps_b(t)}
| \eps_a(t) \rangle \nonumber\\
&= \sum_{a,b^{\prime}} \delta_{\omega, \eps_{b^{\prime}}(t) - \eps_{a}(t)} \delta_{\omega, \eps_b(t) - \eps_a(t)}\bra{\eps_{b^{\prime}}(t)} A_\alpha \ket{\eps_{a}(t)}  \bra{\eps_a(t)} A_\alpha \ket{\eps_b(t)}
| \eps_{b^{\prime}}(t) \rangle  \nonumber\\
&= \sum_{a,b^{\prime}} \delta_{0, \eps_{b^{\prime}}(t) - \eps_{b}(t)} \delta_{\omega, \eps_b(t) - \eps_a(t)}\bra{\eps_{b^{\prime}}(t)} A_\alpha \ket{\eps_{a}(t)}  \bra{\eps_a(t)} A_\alpha \ket{\eps_b(t)}
| \eps_{b^{\prime}}(t) \rangle   \, ,
\label{eq:firstterm}
\end{align}
where the sum over $b^'$ denotes the sum over $|_b^'⟩$ sharing the same energy as $|_b⟩$.
The second term is:
\begin{align}
&\phantom{{}={}} \bra{\eps_b(t)} A_{i,\omega}^{\dagger}(t)A_{i,\omega}(t) \ket{\eps_b(t)} \ket{\eps_b(t)}  \nonumber\\
&= \left(\sum_{a^{\prime}} \delta_{\omega, \eps_b(t) - \eps_{a^{\prime}}(t)} \bra{\eps_b(t)} A_\alpha \ket{\eps_{a^{\prime}}(t)} \langle \eps_{a^{\prime}}(t)| \right) \left(\sum_{a} \delta_{\omega, \eps_b(t) - \eps_a(t)} \bra{\eps_a(t)} A_\alpha \ket{\eps_b(t)}
| \eps_a(t) \rangle \right)\ket{\eps_b(t)} \nonumber\\
&= \left(\sum_{a} \delta_{\omega, \eps_b(t) - \eps_{a}(t)} \bra{\eps_b(t)} A_\alpha \ket{\eps_{a}(t)} \bra{\eps_a(t)} A_\alpha \ket{\eps_b(t)} \right)\ket{\eps_b(t)} \,.\nonumber\\
\label{eq:secondterm}
\end{align}

Subtracting Eq.~\eqref{eq:secondterm} from Eq.~\eqref{eq:firstterm} yields the drift term [Eq.~\eqref{eq:paraterm}], which is not zero, but a linear combination of degenerate eigenstates with the same energy $_b(t)$. Before the jump happens the environment leads to the redistribution of $|_b(t)⟩$ to other states in the same energy manifold. (The Lamb shift $H_LS(t) = ∑_i,ωS_i(ω) A^†_i,ω(t) A_i,ω(t) $ also yields the same effect.) Since they all share the same energy, this does not affect the overlap with the ground state. If the evolution by $H_S(t)$ is adiabatic, such a linear combination will stay in the same energy manifold and this explains the square-pulse like behavior in the overlapping with the ground state in Fig.~\ref{fig:8qubitssingle} of the main text.

\subsubsection{No degeneracy in $\eps_b(t)$}
If there are no degenerate states with energy $_b(t)$,
\begin{align}
&\phantom{{}={}} A_{i,\omega}^{\dagger}(t) A_{i,\omega}(t)\ket{\eps_b(t)}  = \sum_{a} \delta_{\omega, \eps_b(t) - \eps_a(t)}\bra{\eps_{b}(t)} A_\alpha \ket{\eps_{a}(t)}  \bra{\eps_a(t)} A_\alpha \ket{\eps_b(t)}
| \eps_{b}(t) \rangle
\end{align}
This cancels with $⟨_b(t)| A_i,ω^†(t)A_i,ω(t) |_b(t)⟩ |_b(t)⟩$ [Eq.~\eqref{eq:secondterm}] and the drift term [Eq.~\eqref{eq:paraterm}] becomes zero.
\newline

\subsection{Derivation of Eq.~\eqref{eq:inctunnrate}}
\label{app:lindbladfermi}

When the state is $|ψ_1(t)⟩$, its jump rate $λ_1 →0(t) $ to $|ψ_0(t)⟩$ comprises the summation of Lindblad terms responsible for the $1 →0$ transition:
\begin{equation}
\lambda_{1 \rightarrow 0}(t) = \sum_{\alpha \in \{1\rightarrow 0\}}\braket{{A}_\alpha^{\dagger}(t){A}_{\alpha}(t)} \,.
\end{equation}
The summation is over the number of qubits $n$. Since each qubit is coupled to its own environment with an independent noise source, the $γ$ matrix in Eq.~\eqref{eq:unitarytransform} is already diagonal. From Eq.~\eqref{eq:Lindblad2}
we know that
\begin{equation}
L_{\alpha, \omega_{10}}(t) = \sum_{a,b} \delta_{\omega_{10}, \eps_b(t) - \eps_a(t)} \bra{\eps_a(t)} \sigma^{z}_{\alpha} \ket{\eps_b(t)} | \eps_a(t) \rangle  \langle \eps_b(t)|    \,.
\end{equation}

Assume that Bohr frequency $ω_10(t)$ is due only to the $1→0$ transition (even if it is not, the other terms would be annihilated by the matrix element $⟨ψ_1(t) |…|ψ_1(t) ⟩$). The Lindblad operators $A_α(t)$ have the form:
\begin{equation}
{A}_{\alpha}(t) = \sqrt{\gamma_{\alpha}(\omega_{10})}\langle\psi_0(t)|\sigma^{z}_{\alpha} |\psi_1(t)\rangle| \psi_0(t) \rangle  \langle \psi_1(t)|  \,.
\end{equation}
\begin{align}
\lambda_{1 \rightarrow 0}(t) &= \sum_{\alpha}  \langle \psi_1(t) | {A}^{\dagger}_{\alpha}(t)  {A}_{\alpha}(t)   |\psi_1(t) \rangle    \\
&=  \sum_{\alpha=1}^{n} \gamma_\alpha(\omega_{10})|\langle\psi_0(t)|\sigma^{z}_{\alpha} |\psi_1(t)\rangle|^2 \nonumber  \,.
\end{align}
Here $γ_α(ω_10)$ is evaluated with respect to the Ohmic spectral density.

\section{Experimental iterative reverse annealing}
\subsection{\label{append:spectrum}Spectrum of the $p$-spin problem with $p=2$ and scaling of the minimum gap with $N$}
The spectrum of $p$-spin problem and the ordering of energies of different spin sectors can be found in~\cite{Bapst2012}. For the particular case of $n=4, p=2$ and the D-Wave 2000Q annealing schedule, we plot the spectrum in Fig.~\ref{fig:myspectrum}. From Fig.~\ref{fig:myspectrum} we can see that the instantaneous first excited state $ϵ_1(s)$ and the instantaneous ground state $ϵ_0(s)$ converge at $s=1$.

The energy gap between instantaneous eigenenergies $ϵ_i(s)$ and $ϵ_j(s)$ is $Δ_ij(s) = ϵ_i(s) - ϵ_j(s)$. The corresponding minimum energy gap is $Δ_ij = min_s Δ_ij(s)$. For $p=2$, we are interested instead in the minimum energy gap between the ground state $ϵ_0(s)$ and the second excited state $ϵ_2(s)$, since the $1$st excited state $ϵ_1(s)$ and the ground state converge at $s=1$. It is denoted by $Δ$.
\begin{equation}
\Delta = \Delta_{20} = \min_{s} \Delta_{20}(s) = \min_{s} \epsilon_{2}(s) - \epsilon_{0}(s) \,.
\end{equation}
We calculate and plot in Figure~\ref{fig:RA_20_size}(c) the value of $Δ$, for $N = {4, ⋯, 22}$. We also show the $s$ where the minimum gap is found, i.e.
\begin{equation}
s_{\mathrm{min}} = \mathrm{argmin} \Delta_{20}(s) \,.
\end{equation}

\begin{figure}[h!]
\centering
\includegraphics[width=0.6\linewidth]{ra2/pspin_p2_spectrum_4.pdf}
\caption{Full spectrum of 4 qubits.}
\label{fig:myspectrum}
\end{figure}

% \begin{figure}[h!]
%     \centering
%     \includegraphics[width=0.6\linewidth]{ra2/gaprevealplot_3.pdf}
%     \caption{$\Delta$ for $N = \{4, \cdots, 22\}$. Inset: the value of $s_{\mathrm{min}}$ for each $N$.}
%     \label{fig:gapvsN}
% \end{figure}

\subsection{\label{append:successbound}Bound on success probabilities, for closed system and open-system dynamics with collective system-bath interaction}
If the dynamics preserve spin symmetry (for example in closed system Schr\"{o}dinger equation and open-system simulation with the collective system-bath coupling assumption), there exists a natural upper bound on the maximum success probabilities achievable in reverse annealing. The upper bound is the population of the initial state in the maximum spin sector.

% We see that the collective system-bath coupling model gives success probabilities lower than that of independent system-bath coupling model. For the computational initial state of $\ket{0}^{\otimes N-1}\ket{1}$, only $\frac{1}{N}$ of its population is in the maximum spin subspace. Recall also that the (instantaneous) ground state and the first excited state belong to maximum spin subspace. Since the collective system-bath coupling model preserves spin symmetry and the dynamics is restricted to each spin subspace, at most $\frac{1}{N}$ of the initial population can be relaxed to the instantaneous ground state and the instantaneous first excited state, and get to the correct solution at the end of the reverse anneal.

% For the example of $n=4$, there are two degenerate ground states and they belong to the maximum spin subspace, $S=\frac{n}{2}=2$. For $H_0$ all the bases of the maximum spin subspace are
% ${\scriptstyle\left\{\ket{\uparrow\uparrow\uparrow\uparrow}, \ket{\downarrow\downarrow\downarrow\downarrow}, \frac{\ket{\uparrow\uparrow\uparrow\downarrow}+\ket{\uparrow\uparrow\downarrow\uparrow}+\ket{\uparrow\downarrow\uparrow\uparrow}+\ket{\downarrow\uparrow\uparrow\uparrow}}{2},\right. \\ \left. \frac{\ket{\uparrow\downarrow\downarrow\downarrow}+\ket{\downarrow\uparrow\downarrow\downarrow}+\ket{\downarrow\downarrow\uparrow\downarrow}+\ket{\downarrow\downarrow\downarrow\uparrow}}{2}, \frac{\ket{\uparrow\uparrow\downarrow\downarrow}+\ket{\uparrow\downarrow\uparrow\downarrow}+\ket{\uparrow\downarrow\downarrow\uparrow}+\ket{\downarrow\uparrow\uparrow\downarrow}+\ket{\downarrow\uparrow\downarrow\uparrow}+\ket{\downarrow\downarrow\uparrow\uparrow}}{\sqrt{6}}\right\}}$,
% $\{\ket{0000}, \ket{1111}, \frac{\ket{0001}+\ket{0010}+\ket{0100}+\ket{1000}}{2}, \frac{\ket{0111}+\ket{1011}+\ket{1101}+\ket{1110}}{2},\\ \frac{\ket{0011}+\ket{0101}+\ket{0110}+\ket{1001}+\ket{1010}+\ket{1100}}{\sqrt{6}}\}$.

For the example of $N=4$, there are two degenerate ground states ($|0000⟩$ ($|0⟩ ≡|↑⟩$, $|1⟩ ≡|↓⟩$) and $|1111⟩$) of the problem Hamiltonian ($H_0$) and they both belong to the maximum spin subspace of $S=S_max=N/2=2$. While the computational basis state with one spin down, for example, $|0001⟩$ is a first excited state of $H_0$, it is not a basis of the maximum spin subspace $S=2$. A superposition of the computational basis states with one spin down, i.e. $|0001⟩+|0010⟩+|0100⟩+|1000⟩/√(4)$, however belong to the maximum spin subspace $S=2$, and this is the case for~\cite{Passarelli2019}. Unfortunately, the D-Wave machine does not allow such an initialization.

In general, suppose that the initial computational basis state is $|ψ(t=0)⟩_comp$, with a particular (normalized) magnetization $m_0$ (Eq.~\ref{eq:mo}). It is represented as a linear combination of states, each of which has a fixed value of total spin momentum $S$ and (normalized) magnetization $m_0$,
\begin{equation}
\label{eq:compinitialstate}
\ket{\psi(t=0)}_{\text{comp}} = \sum_{S=S_{\text{min}}}^{S_{\text{max}}}a_S\ket{S, m_{S}=m_0}.
\end{equation}
% \begin{equation}
%     \ket{\psi(t=0)}_{\text{comp}} = a \ket{S=S_{\text{max}},m=m_0} + b\ket{S=S_{\text{max}}-1,m=m_0} + \dots.
% \end{equation}
From the theory of addition of angular momentum, we know that $S_min = N/2-⌊N/2 ⌋$ and $S_max = N/2$. The total spin momentum (integer or half-integer) $S ∈{S_min, S_min + 1, …, S_max}$. $|S, m_S=m_0⟩$ is a simultaneous eigenstate of $S^2 = ∑_α∈{x,y,z}S^α$ and $S^z=1/2∑_i=1^Nσ_i^z$ (or $H_p$), with eigenvalues:
\begin{align}
\bold{S}^2\ket{S, m_{S}=m_0} &= S(S+1)\ket{S, m_{S}=m_0} \,,\\
S^z\ket{S, m_{S}=m_0} &= \Big(\frac{N}{2}m_0\Big)\ket{S, m_{S}=m_0}
\end{align}
In Eq.~(\ref{eq:compinitialstate}), $|a_S|^2$ is the initial state's population in that spin subspace.
% $\ket{S=S_{\text{max}}-k}$ is a basis of spin subspace $S=S_{\text{max}}-k$ and $|a|^2$ is the initial state's population in the maximum spin subspace $S_{\text{max}}$.
% Projected to $\bra{S=S_{\text{max}}}$,
% $\braket{S=S_{\text{max}}|0001} = a$.

% Since we know the rhs to be $\frac{1}{2}$ for $N=4$, we find $a=\frac{1}{2}$.

For the $p$-spin Hamiltonian, in a closed system the state in each spin subspace develops independently due to the preservation of spin symmetry of the Hamiltonian~\cite{nishimori:reverse-pspin-2}. At the end of a single cycle of reverse annealing ($r=1$), we have from Eq.~(\ref{eq:compinitialstate}) that
\begin{equation}
U(2t_\text{inv},0)\ket{S, m_{S}=m_0} = \sum_{M=-S}^{S} c_{M} \ket{S, m_{S}=2M/N},
\end{equation}
where $c_M$ is the relative phase developed in the other bases of the spin $S$ subspace at the end of the time evolution. The (unnormalized) magnetization $M = (N/2m_S) ∈{-S, -S+1, …, S}$.

Therefore, the final state of $1$-cycle reverse annealing (denoted by $|ψ⟩_fin$) can be expressed as:
\begin{align}
\ket{\psi}_{\text{fin}} &= U(2t_\text{inv},0)\ket{\psi(t=0)}_{\text{comp}} \nonumber\\
&= \sum_{S=S_{\text{min}}}^{S_{\text{max}}}a_S\left(\sum_{M=-S}^{S} c_M \ket{S, m_{S}=2M/N}\right)\,.
\end{align}
We know that the all-up state ($|0⟩^⊗N = |↑>^⊗N = |up⟩$) and all-down state ($|1⟩^⊗N = |↓>^⊗N = |down⟩$) are both solution states of $H_0$, and moreover they lie in the maximum spin subspace $S=S_max$. Particularly, $|up⟩ = |S_max, 1⟩$ and  $|down⟩ = |S_max, -1⟩$. Therefore, projecting $|ψ⟩_fin$ onto the $|up⟩$ gives:
\begin{align}
&=  \sum_{S=S_{\text{min}}}^{S_{\text{max}}}a_S\left(\sum_{M=-S}^{S} c_{M} \braket{\text{up}|S, m_{S}=2M/N}\right)\\
&= a_{S_\text{max}}c_{S_\text{max}}\,.
\end{align}
Similarly, $⟨down|ψ|_⟩fin = a_S_maxc_-S_max$.

The total success probability is thus bounded by:
\begin{align}
&\phantom{==}|\braket{\text{up}|\psi}_{\text{fin}}|^2 + |\braket{\text{down}|\psi}_{\text{fin}}|^2 \nonumber\\
&= |a_{S_\text{max}}|^2(|c_{S_\text{max}}|^2+|c_{-S_\text{max}}|^2) \nonumber\\
&\leq |a_{S_\text{max}}|^2 \,,
\end{align}
where equality is met with %$\ket{S=S_{\text{max}}}_{\text{fin}} = \ket{\text{up}}$.
\newline \rightline{$\square$}

The upper bound $|a_S_max|^2$ is the population of the initial state in the maximum spin subspace. For the initial state of $|0001⟩$, we have $a_S_max = 1/√(4)$ since $|S=S_max=2, m_S=0.5⟩ = |0001⟩+|0010⟩+|0100⟩+|1000⟩/√(4)$. Therefore, the total success probability is bounded by $|a_S_max|^2 = 1/4$. For the other examples in the main text, the initial state of $|0011⟩$ has $a_S_max = 1/√(42) = 1/√(6)$; while the initial state of $|00000001⟩$ has $a_S_max = 1/√(8)$.

By expressing the above equations in terms of density matrix, we can also generalize the same conclusion to open-system simulation of collective system-bath coupling, where the population of each spin-sector is preserved during dynamics.

\subsection{\label{append:svmc}SVMC and SVMC-TF for reverse annealing}
We outline the algorithm of SVMC and SVMC-TF with application to reverse annealing.
\begin{spacing}{0.8}
\begin{algorithm}[H]
\scriptsize
\label{alg:1}
\caption{SVMC and SVMC-TF ($p$-spin reverse annealing)}
\hspace*{\algorithmicindent}
% \textbf{Input}: $\vec{v}$, $\vec{w}$ and $W_\mathrm{max}$.
\begin{algorithmic}[1]
% \State $\ket{\psi(0)}$ denotes the  uniform superposition of all the allowed states satisfying
% the weight constraint: $\sum_{i=1}^n w_i x_i  = 	W_x \leq W_\mathrm{max}$.
% \State Choose a number of blocks $p$ for the QAOA.
% \State Choose a number of samples/shots $S$.
\State $\mathcal{H}(s) = -\frac{A(s)}{2} (\sum_{i}^{N} \sin \theta_i) - \frac{B(s) N}{2}\left(\frac{1}{N}\sum_{i}^{N} \cos \theta_i\right)^{p}$, with a specified reverse annealing $s=t$ dependence $s(t)$.
% \For{every $\vec{\beta} \in [0, \pi]^{p}$ and $\vec{\gamma} \in [0, 2\pi]^{p}$}
\Procedure{SVMC}{}
\For{$k=1$ to $K$ (number of samples)}
\State Initial computational basis state: $\ket{0}_i\rightarrow \theta^{k}_{i,t}: 0$, $\ket{1}_i\rightarrow \theta^{k}_{i,t}: \pi$.
\For{$t=1$ to $T$ (total number of sweeps)}
\For{$i=1$ to $n$ (number of qubits)}
%\State Randomly choose $i \in [1, n]$.
\State Randomly choose a new angle $\theta_{i,t}^{'k} \in [0, \pi]$.
\State Calculate $\Delta E$.
\If{$\Delta E \leq 0$}
\State $\theta_{i,t}^{k} \rightarrow \theta_{i,t}^{'k}$
\ElsIf {$p < \exp (-\beta\Delta E)$ where $p\in [0,1]$ is drawn with uniform probability.}
\State $\theta_{i,t}^{k} \rightarrow \theta_{i,t}^{'k}$
\EndIf
\EndFor
\EndFor
\State Take the mean of $K$ samples: $\bar{\theta}_{i,t} = (\sum_{k=1}^{K}\theta_{i,t}^{'k})/K$.
\EndFor
\EndProcedure
\State \Return $\bar{\theta}_{i,t}$.
\Procedure{SVMC-TF}{}
\For{$k=1$ to $K$ (number of samples)}
\State Initial computational basis state: $\ket{0}_i\rightarrow \theta^{k}_{i,t}: 0$, $\ket{1}_i\rightarrow \theta^{k}_{i,t}: \pi$.
\For{$t=1$ to $T$ (total number of sweeps)}
\For{$i=1$ to $n$ (number of qubits)}
%\State Randomly choose $i \in [1, n]$.
\State $\theta_{i,t}^{'k} = \theta_{i,t}^{k} + \epsilon_i(s(t/T))$,
\State a random $\epsilon_i(s(t/T) \in [-\min \left(1,\frac{A(s(t/T)}{B(s(t/T)}\right) \pi, \min \left(1,\frac{A(s(t/T)}{B(s(t/T)}\right) \pi]$
\State Calculate $\Delta E$.
\If{$\Delta E \leq 0$}
\State $\theta_{i,t}^{k} \rightarrow \theta_{i,t}^{'k}$
\ElsIf {$p < \exp (-\beta\Delta E)$ where $p\in [0,1]$ is drawn with uniform probability.}
\State $\theta_{i,t}^{k} \rightarrow \theta_{i,t}^{'k}$
\EndIf
\EndFor
\EndFor
\State Take the mean of $K$ samples: $\bar{\theta}_{i,t} = (\sum_{k=1}^{K}\theta_{i,t}^{'k})/K$.
\EndFor
\EndProcedure
\State \Return $\bar{\theta}_{i,t}$.
% \EndFor
% \State Angle searching: Search for the angles %$({\vec{\beta}_{\text{opt}}, \vec{\gamma}_{\text{opt}}})$,
% $({\vec{\beta}^{*}, \vec{\gamma}^{*}})$,
% that maximize $\overline{V(x)}$.
% \State Calculation step: Calculate the \textit{solution} based on the bit string distribution produced by QAOA with $({\vec{\beta}_{\text{opt}}, \vec{\gamma}_{\text{opt}}})$.
% \State
% \begin{align*}
% \Delta E &=  -A(t/T) \sum_i \sin \theta^{'}_i  + B(t/T) \left( - \sum_i h_i \cos \theta^{'}_i + \sum_{i,j} J_{ij} \cos \theta^{'}_i \cos \theta^{'}_j \right) \\
% &- -A(t/T) \sum_i \sin \theta_i  + B(t/T) \left( - \sum_i h_i \cos \theta_i + \sum_{i,j} J_{ij} \cos \theta_i \cos \theta_j \right)
% \end{align*}

\Procedure{Projection onto computational basis}{}
\If{$0\leq\theta_{i,t}^{k}\leq \pi/2$}
\State $\ket{\psi^k(t)}_i = \ket{0}$
\ElsIf {$\pi/2\leq\theta_{i,t}^{k}\leq \pi$}
\State $\ket{\psi^k(t)}_i = \ket{1}$
\EndIf
\EndProcedure
% \Procedure{Projection2 (Computational basis)}{}
% \If{$0\leq\bar{\theta}_{i,t}\leq \pi/2$}
%     \State $\ket{\psi(t)}_i = \ket{0}$
% \ElsIf {$\pi/2\leq\bar{\theta}_{i,t}\leq \pi$}
%     \State $\ket{\psi(t)}_i = \ket{1}$
% \EndIf
% \EndProcedure
% \Procedure{Projection3 (Computational basis)}{}
% \If{$\lfloor\cos(\bar{\theta}_{i,t})\rfloor = 0$}
%     \State $\ket{\psi(t)}_i = \ket{0}$
% \ElsIf {$\lfloor\cos(\bar{\theta}_{i,t})\rfloor = -1$}
%     \State $\ket{\psi(t)}_i = \ket{1}$
% \EndIf
% \EndProcedure
\State Remark: $\Delta E$, for example, due to the update of qubit $2$, can be expressed in the following form:
\begin{align*}
\phantom{=}&\Delta E_2 \\
&\left(-\frac{A(s)}{2} \left(\sum_{\substack{i=1\\i \ne 2}}^{N}\sin \theta_i + \sin \theta^{'}_2\right) - \frac{B(s)N}{2}\left(\frac{1}{N}\left(\sum_{\substack{i=1\\i \ne 2}}^{N} \cos \theta_i + \cos \theta^{'}_2\right)\right)^{2}  \right) \\
- &\left(-\frac{A(s)}{2} \left(\sum_{i=1}^{N}\sin \theta_i\right) - \frac{B(s)N}{2}\left(\frac{1}{N}\left(\sum_{i=1}^{N} \cos \theta_i\right)\right)^{2} \right) \,.
\end{align*}
\end{algorithmic}
\end{algorithm}
\end{spacing}

% For example, the update of qubit $2$ to a new angle results in a change in energy, denoted by $\Delta E_2$, as:
%\phantom{=}&\Delta E_2 \\
%&\left(-\frac{A(s)}{2} \left(\sum_{\substack{i=1\\i \ne 2}}^{4}\sin \theta_i + \sin \theta^{'}_2\right) - \frac{B(s)}{8}\left(\sum_{\substack{i=1\\i \ne 2}}^{4} \cos \theta_i + \cos \theta^{'}_2\right)^{2} \right) \\
%- &\left(-\frac{A(s)}{2} \left(\sum_{i=1}^{4}\sin \theta_i\right) - \frac{B(s)}{8}\left(\sum_{i=1}^{4} \cos \theta_i\right)^{2} \right) \,.

% \subsection{\label{append:svmctflarge}SVMC-TF results for $N=16,32$}

% In Fig.~\ref{fig:partial16} we report the simulation results of partial success probability of $N=16$ reverse annealing simulated by SVMC-TF, with $\tau = 10^3$ (sweeps) and $10^4$ (sweeps), respectively. The number of samples is $10^4$.

% \begin{figure}[h!]
% % \subfigure[]{\includegraphics[width = 0.4939\columnwidth]{svmcplot_100sweeps.pdf}\label{fig:4qubitsunitarya}}
% \subfigure[]{\includegraphics[width = 0.4939\columnwidth]{ra2/svmcplot_1000sweeps_n16_ge.pdf}}
% \subfigure[]{\includegraphics[width = 0.4939\columnwidth]{ra2/svmcplot_10000sweeps_n16_ge.pdf}}
% \caption{Partial success probability for SVMC-TF. $N=16$. (a) $\tau = 10^3$ sweeps, (b) $\tau = 10^4$ sweeps. Error bars are 2$\sigma$ over $10^4$ samples.}
% \label{fig:partial16}
% \end{figure}

% In Fig.~\ref{fig:partial32} we report the simulation results of partial success probability of $N=32$ reverse annealing simulated by SVMC-TF, with $\tau = 10^3$ (sweeps) and $10^4$ (sweeps), respectively. The number of samples is $10^4$.

% \begin{figure}[h!]
% % \subfigure[]{\includegraphics[width = 0.4939\columnwidth]{svmcplot_100sweeps.pdf}\label{fig:4qubitsunitarya}}
% \subfigure[]{\includegraphics[width = 0.4939\columnwidth]{ra2/svmcplot_1000sweeps_n32_ge.pdf}}
% \subfigure[]{\includegraphics[width = 0.4939\columnwidth]{ra2/svmcplot_10000sweeps_n32_ge.pdf}}
% \caption{Partial success probability for SVMC-TF. $N=32$. (a) $\tau = 10^3$ sweeps, (b) $\tau = 10^4$ sweeps. Error bars are 2$\sigma$ over $10^4$ samples.}
% \label{fig:partial32}
% \end{figure}

\endgroup

% In case your dissertation has multiple volumes.
% \addvolumecontents{thesis_part2}
% \addvolumecontents{thesis_part3}
% \addvolumecontents[lof]{thesis_part2}

\end{document}