Variational Matrix Product Operators for the Steady State of Dissipative Quantum Systems

0.55 MB·English

by  Jian Cui

Variationalmatrixproductoperatorsforthesteadystateofdissipativequantumsystems
Jian Cui1, J. Ignacio Cirac2, Mari Carmen Ban˜uls 2
1 QOLS, Blackett Laboratory, Imperial College London, SW7 2BW, United Kingdom and
2 Max-Planck-Institutfu¨rQuantenoptik, Hans-Kopfermann-Str. 1, D-85748Garching, Germany
(Dated:June10,2015)
Wepresentanewvariationalmethod, basedonthematrixproductoperator(MPO)ansatz, forfindingthe
steady state of dissipative quantum chains governed by master equations of the Lindblad form. Instead of
requiringanaccuraterepresentationofthesystemevolutionuntilthestationarystateisattained,thealgorithm
directlytargetsthefinalstate,thusallowingforafasterconvergencewhenthesteadystateisaMPOwithsmall
bonddimension.Ournumericalsimulationsforseveraldissipativespinmodelsoverawiderangeofparameters
illustratetheperformanceofthemethodandshowthatindeedthestationarystateisoftenwelldescribedbya
MPOofverymoderatedimensions.
5
1 Introduction.— The physics of quantum systems out of quantumtrajectories,thelatterhavealsobeenappliedtodis-
0 equilibrium poses unsolved fundamental questions, relating sipative dynamics [30]. The natural extension to operators,
2
to Nature at extreme conditions and to the dynamics after namely matrix product operators (MPO), can be used as an

n longtimeevolution. Progressinthisfieldishoweverhardto ansatz for mixed states [31, 32], which is known to accu-
u
achieve,duetothelackofanalyticaltoolstosolvemanysuch ratelydescribethermalequilibriumstatesforlocalHamiltoni-
J
problems,andthelimitationsofexistingnumericalmethods. ans[33,34]. Suchanextension,incombinationwiththetime
9
In recent times growing attention has been directed to the evolution algorithms, has allowed the numerical exploration

] out-of-equilibriumphysicsofopenquantumsystems,i.e.sys- ofsteadystatesofspinchainsandotheronedimensionalsys-
h tems in interaction with an environment. This interest has temsunderlocaldissipation(seee.g. [3,35–40]).
p
been intensified by the potential applications to the fields of Thismethodisformallysimilartothesearchforaground
-
t condensedmatterphysics,statisticalphysics,andquantumin- stateusingimaginarytimeevolution[26],inthatagivenini-
n
formation processing [1–4]. In particular it has been shown tial state is evolved until reaching a fixed point of the dy-
a
u that dissipation can be used to engineer interesting quantum namics. However, different to the imaginary time evolution
q many-body states and to perform universal quantum compu- method,wherethesequenceofstatesvisitedbythealgorithm
[
tation[1,5],ideaswhichcanbeexploredinthecontextofcur- is not of physical significance, in the simulation of a mas-

2 rentexperimentalsetupsbasedonatomicsystems[6]. Apar- ter equation, the real evolution needs to be followed, so that
v ticularlyinterestingtopicisthatofdissipativequantumphase errorsintheintermediatestatecanseverelyaffecttheconver-
6 transitions(DQPT),namelytransitionsinthenon-equilibrium genceoftheprocedure.
8
steadystateofanopensystem,whichmayarisefromthecom- Abetteralternativecouldbegivenbyavariationalmethod,
7
6 petingeffectsoftheHamiltonianandthedissipativetermsof whichsearchesforthenullvectoroftheLindbladsuperoper-
0 the dynamics. An archetypical example is that of the model ator within the MPO family, in the spirit of the DMRG vari-
. [7,8],butDQPThavealsobeenstudiedinfermionic[9–11], ationalsearch. Suchamethodwouldbepotentiallymoreef-
1
0 bosonic[12]andquantumspinsystems[13–15]. ficient than simulating the full evolution, specially when the
5 Findingthestationarystateofagenericmasterequationis latter traverses intermediate states with a large bond dimen-
1 noteasy,evenfor1Dsystems. Analyticaltreatmentislimited sion, butthetruesteadystateisdescribedbyasmallone, as
:
v toveryspecificproblems,suchasquadraticfermionicmodels is often the case [3, 36, 41]. In this paper we present such
i [16], or systems under special conditions or approximations variationalmethodforthesteadystateofamasterequationin
X
[17,18],andmostoftennumericaltechniquesarenecessary. Lindbladform.Weillustratetheperformanceofthealgorithm
r
a As in the case of pure states an exact numerical treatment with results for several one dimensional models. Notice that
is possible only for small systems, due to the exponentially avariationalmethod,similarinspiritbutrestrictedtodensity
growingcomputationalcost,whichmaybeevenmoresevere matricescontainingonlyfew-bodycorrelations, hasbeenre-
in the case of mixed states. For pure states, parametrizing centlyproposedin[42].
thestateasaTensorNetwork(TN)[19–21]hasprovenanef- Basicconcepts.—Amatrixproductstate(MPS)foraquan-
ficient alternative, that can successfully capture the physical tum system of N d-dimensional components, is a state vec-
properties of quantum many body states in countless situa- torofthe form∣Ψ⟩ = ∑ tr(As1...AsN)∣s ...s ⟩ [43],
tionsofinterest. Thebestexampleisthetremendoussuccess whereeachA isad×D{s×i}Dten1sor,DiNsapa1rametNerofthe
i
ofthedensitymatrixrenormalizationgroup(DMRG)[22,23], representationcalledbonddimension,andthesumrunsover
based on the matrix product state (MPS) ansatz, which pro- all elements of each individual basis s = 1,...d. By suc-
i
vides a quasi-exact solution for one dimensional problems. cessivelyincreasingthebonddimension,D,theMPSfamily
MPS can accurately describe ground states of gapped local definesahierarchyofstatescoveringvectorspacespannedby
Hamiltonians[24,25],andmethodshavebeendefinedtouse the tensor product of the individual bases, {∣s ⟩}. The same
i
themalsoinrealtimeevolution[26–29]. Incombinationwith ansatz can be used to represent operators whose coefficients
2
inatensorproductbasishavethestructureofamatrixprod- ofthisoperator,and,sinceLˆ†Lˆ⪰0,itcorrespondstothelow-
uct, Oˆ = ∑ tr(As1r1...AsNrN)∣s ...s ⟩⟨r ...r ∣ esteigenvalue. IfLˆcanbewrittenasaMPO,alsotheproduct
{si,ri} 1 N 1 N 1 N
Thesearecalledmatrixproductoperators(MPO)[31,32,44]. can, and it is then possible to use standard MPS algorithms
The operators can be vectorized using Choi’s isomorphism, toapproximateitsgroundstate[20,23]. Noticetheparticular
∣s ⟩⟨r ∣ → ∣s r ⟩, which maps any operator Oˆ to a vector case of Hermitian L is especially easy, since the (properly
i i i i α
∣Φ(O)⟩, so that it is possible to work in the vector space of normalized) identity is a steady state, which can be exactly
operatorswiththeusualMPStechniques. writtenasaMPSwithbonddimensionD=1.
Inordertodescribephysicalmixedstates,MPO,orinthis Thefactthatwearetargetingdensitymatricesrequirespar-
case matrix product density operators (MPDO), have to sat- ticularattention,becausenoteveryMPSvectorcanrepresent
isfy additional conditions, namely they have to be normal- a valid physical state. The normalization condition trρ = 1
ized (trρ = 1), Hermitian and positive semidefinite. While translates to ⟨Φ(1)∣Φ(ρ)⟩ = 1, where ∣Φ(1)⟩ is the (unnor-
the first two conditions are easy to impose on the local ten- malized)vectorthatcorrespondstothetracemap,namelythe
sorsoftheansatz,thepositivityinvolvesthefullspectrumof maximallyentangled∣Φ(1)⟩=∑ ∣s ...s ⟩⊗∣s ...s ⟩.
{si} 1 N 1 N
theoperator,andisthusanon-localproperty. Theansatzcan A solution which is not orthogonal to this vector can always
bemodifiedtorepresentpositiveoperators, usingalocalpu- be normalized to ensure the trace condition. In general it is
rification of the state with MPS structure. In this case, each morecomplicatedtodecidewhetheraMPScorrespondstoa
tensorhasastructureAij = ∑ Xik⊗X¯jk, wheretheindex positiveoperator,sincewedonothaveaccesstothefullspec-
n k n n
ksumsovertheancillarydegreeoffreedomandthebarindi- trum. The purification ansatz can guarantee that the search
catescomplexconjugation. Althoughitguaranteespositivity, runsoveronlypositiveoperators, butattheexpenseofmore
working with the purification ansatz is in general computa- costlylocaloptimizations[45]. Henceweusesimplythevec-
tionally more costly [45] and moreover the bond dimension torizedMPOformandrelyonthemathematicalpropertiesof
required to write the purification ansatz may be much larger the problem to provide a physical solution. Since the evolu-
thanthatoftheMPO[46],sothatinpracticeitisnotalways tiongeneratedbyLisaCPmap,itmusthaveapositivefixed
themostconvenientchoice. point, so that if this is non-degenerate, the algorithm should
A variational search for the steady state.— We consider naturally converge towards a MPO approximation of a posi-
a chain of length N, with a quantum system of physical di- tive(andhenceHermitian)operator2 andisthenexpectedto
mension d on each site, and dynamics governed by a master bealmostpositive,withanynon-positivenessbeingcompati-
equation of Lindblad form, dρ = L[ρ], where the rhs is the blewiththetruncationerror.
dt
Lindbladiansuperoperator, Inpractice,wefindthatasuitablewarmupphase(seeSup-
plemental Material [48]) allows us to avoid solutions with
1
L[ρ]=−i[H,ρ]+∑2(2LαρL†α−{L†αLα,ρ}). (1) vanishing trace, and improves the convergence of the algo-
α rithm.3 Although positivity cannot be checked explicitly (it
Theunitarypartoftheevolutionisdeterminedbythesystem is in fact a hard problem [49]), there is a number of neces-
Hamiltonian, H. The effect of the environment is described sary criteria that any physical state needs to satisfy, such as
byasetofLindbladoperators,L . physicallysensiblevaluesofallsinglebodyobservables. Our
α
TheLindbladianactslinearlyonthevectorizedρ,as algorithmperformsasetofsuchtestsandonlyacceptssolu-
tionsthatpassthemall,otherwiserestartingthesearchwitha
Lˆ=−i(H⊗1+1⊗H)
different initial guess. The role of our consistency checks is
+∑1(2L ⊗L¯ −L†L ⊗1−1⊗LTL¯ ). (2) todiscardtheleastsuitableguessesduringthewarmupphase,
2 α α α α α α in order to prepare a suitable initial state for the variational
α
search. Duringthelaterphasesofthealgorithm,thetestsare
The steady state is a fixed point of the evolution, dρs = 0,
dt used as assertions, while we rely on the convergence criteria
andcorrespondstoavector∣Φ(ρs)⟩satisfyingLˆ∣Φ(ρs)⟩=0, (includingthatoftheeffectiveenergy)tostopthecalculation.
i.e. a zero eigenvector of Lˆ. If the Hamiltonian and the in- Afterfindinganacceptablesolutionforagivenbonddimen-
dividualLindbladoperatorshavelocalcharacter,theLindbla- sion,D,wecomputethedesiredexpectationvaluesfromthe
diancanbewrittenasaMPO,1andwecansearchforthebest hermitianpartoftheMPO,(ρ+ρ†)/2,asoftendoneinother
MPSapproximationtoitszeroeigenvector,whichwillgiveus
a vectorized MPO approximation for the steady state. Since
theoperator(2)isnotHermitian,inordertousethestandard
variational search with MPS we consider instead the Hermi-
2Strictlyspeaking,thismayfailifthealgorithmgetsstuckinlocalminima,
tianproductLˆ†Lˆ. Thesteadystateisalsoazeroeigenvector oriftherearedegeneracies,asisknowntohappenalsoforDMRG(see
e.g.[47]).
3NoticethatitisalsopossibletouseaLagrangemultiplierterm,oftheform
∣Φ(1)⟩⟨Φ(1)∣,tofavorsolutionswithnon-vanishingtrace. Thiscanbe
1Strictlyspeaking,itisenoughthatHandeachLαcanthemselvesbewrit- substractedfromLˆ†Lˆ,thusincreasinginoneunitthebonddimensionof
tenasMPO,whichincludesshortrangeinteractionsanddissipation,but theMPO.Inpractice, wefoundthatthewarmupphasewasenoughto
canalsobeappliedtoapproximatepower-lawdecayingterms[44]. obtainphysicalsolutions,withoutincreasingthecomputationalcost.
3
0.5 1 Numerical results.— To illustrate the performance of the
0.4
g/γ=1 g/γ=0.35 algorithmweapplyittoseveralspinchains,wheretheunitary
0.4 0.4 g/γ=0.625 0.8 0.2 g/γ=0.5 and dissipative dynamics show competing effects regarding
2hSiy00..23 0.20 0.0gg5//γγ==02..055.1 2ρtr()00..46 00 0.01 0gg.//0γγ2==01.6 theLcoowhedreimnceensoifotnhaelstDeaicdkyestmatoed.el. A typical example of
1/N 1/N
0.1 0.2 DQPT is exhibited by the Dicke model [7, 8], in which the
collective interaction with a single radiation mode induces
00 2 4 6 00 0.5 1 1.5 2 coherent behavior on a system of N two-level atoms. The
g/γ g/γ
regime of parameters required to observe the DQPT is chal-
FIG.1.Left:Correlation⟨S2⟩forthelowdimensionalDickemodel, lenging, and the experimental observation of the phase tran-
y
asafunctionofg/γforsystemsizesN =10−100(inincreasingly sitionhasonlybeenachievedrecently[50–52]. Itisthusin-
darkershades).Thesolidlineistheresultoffinitesizeextrapolation, terestingtounderstandthebehaviorofsimilarmodelswhich
linearin1/N,asexplicitlyshownintheinsetforseveralvaluesof
may then be easier to realize experimentally. We consider a
g/γ. Right: Purityoftheconvergedsteadystateforthesamesys-
chain of N two-level systems, where each pair of systems
temsizes.Forlargedissipation(shadedregion)weshowonlyresults
couples coherently to a common radiation mode. This can
forthesmallestsystems,forwhichthesimplevariationalsearchcon-
vergesreliably.Theinsetshowstheexplicitdependenceofthepurity be represented by a spin-1/2 chain, governed by a single-
onthesystemsizeforseveralcouplingvalues. particleHamiltonian, H = ∑N gσx, andLindbladoperators
i=1 i
L = γ(σ− +σ− ), for i = 1,...N −1, instead of the sin-
i i i+1
gle collective Lindblad operator of the Dicke model, so that
algorithms to reduce numerical errors.4 The found solution this model can be considered a low-dimensional version of
is normalized and used as initial guess for a larger bond di- the latter. We study the nature of the steady state found by
mension, and finally convergence in D is decided when the thealgorithmatvaryingvaluesofg/γ andincreasingsystem
targetedobservablesareconvergedtothedesiredprecision. sizes, N up to 100, which allows us to perform a finite size
Thealgorithmasdescribedhereisthusformallyequivalent extrapolation and study single site observables and correla-
tothevariationalgroundstatesearchforaMPOHamiltonian tions in the thermodynamic limit. In the Dicke model, the
overtheMPSfamily,andpresentsthesamescaling,onlywith superradiant phase transition (at g/γ2 = N) [8] is visible in
d2 playingtheroleofthephysicaldimension,andwithanef- theseobservables. Inthelowdimensionalversionwedonot
fectiveHamiltonianwhichhasthesquaredbonddimensionof find evidence of such transition, although (short range) cor-
the MPO for Lˆ. The gap of Lˆ†Lˆwill be determinant for the relations appear, as shown in figure 1 for S2 = (∑ σy/N)2.
y i i
convergence of the algorithm. It is interesting to notice that It is remarkable that for all values g/γ ≥ 0.35 and system
thisisnotrelatedtotheeigenvaluesofLˆ,buttoits(squared) sizes N ≤ 100 the steady state is converged with very small
singular values (see [48]). All in all we find that, for typical bond dimension, D < 30, most of them even with D ≤ 20.
cases, thesmallbonddimensionrequiredtoapproximatethe At g =√0 there are however two dark states, namely ∣0⟩⊗N
steadystateasaMPOcompensatesfortheadditionalcompu- and (1/ N)∑N (−1)k∣0⊗(k−1)10⊗(N−k)⟩. Hence, the null
k=1
tationaleffortassociatedtoLˆ†Lˆ,providedtheLindbladianis subspace of Lˆis four-fold degenerate. This hinders the con-
notdegenerate. vergence of the algorithm at very small g/γ, as the steady
Another variational approach has been recently proposed state is no longer the unique and positive zero eigenvector,
[42], which chooses to minimize the trace norm of L[ρ]. and the warmup strategy is not enough to guarantee a phys-
Since this quantity is not efficiently computable, the method ical solution, except for the smallest system sizes (N ≤ 20).
in[42]proceedsbyminimizinganupper-boundtothisnorm For those converged cases, we can detect the peculiarity of
for restricted sets of density matrices. Our eigenvalue mini- thisparameterregionbyanalyzingthepurityofthesolution,
mization is instead equivalent to finding the vector that min- showninfigure1. Indeed,wecanfindpositivesolutionswith
imizes the Euclidean norm ∥Lˆ∣Φ(ρ)⟩∥ with the constraint increasingpurity,whichcanbeupto1forg=0.5 Inprinciple
∥∣Φ(ρ)⟩∥ = 1. Both minimizations have an exact solution onecouldcomplementthemethodwithadditionaltechniques
in the physical steady state, although they are not equivalent to try and select the physical steady states (e.g. finding and
whennotexactlyonthestationarystate. UsingtheEuclidean then processing several orthogonal eigenstates, not necessar-
norm of the vectorized expression is preferable in our case, ilyphysical)evenforlargerchains.Herewehavenevertheless
becausethetracenormrequiresthefulldiagonalizationofthe focusedontheconvergenceinthemostcommonlyoccurring
operators, impossible for the system sizes we are interested situationofauniquesteadystate,wherethemethodcanpro-
in, while the norm of the vectorized operators is efficiently vide the largest gain by directly targeting a MPO with small
calculablefortheMPSansatz. bonddimension.
4Thenormofthenon-Hermitianpartcanalsobeusedasadditionalconsis- 5Noticethatforg=0anymixtureofbothdarkstateswillbeasteadystate,
tencycheck. withpurityin[0.5,1].
4
2Mpz0000....3456 NNNNN=====1234500000 N/2|C(n)|zz1100−05∆5=∆−=n210.∆54=425 2ρtr()000...4681 zσn0.51 2 ∆∆∆=n==6−404 1 0 egitmnhxetaoitessnrrtetmithnheeegfedfistiatceaacintteteisunooastnrtlaeansrxtoyeepstllswuotwtaoriatohretknii,co.whanpTimtpohhrfauooystuahtcnoehteuheseretsdet,neaatecdhehlyidasnrsitgmqtoaeuetpreetrhbecpoocahdninsadedsaleidylrliedomdciweatelsgyncfrsoratiiramobrnae-.
0.2 Ournumericalresultsshowthatforavariatedsetofmodels,
0.2
0.1 withcorrelationscreatedbytheunitaryevolution,thedissipa-
0
−2 0 2 4 −2 0 2 4 tionorboth[48],thesteadystateisindeedwellapproximated
∆ ∆
byaMPOofverysmallbonddimension,D≤30forsizesup
√
FIG. 2. Left: antiferromagnetic order parameter ⟨M2⟩ in the to N = 100. This can be directly compared to the bond di-
z
steady state of the Ising model with local dissipation, for several mensionsrequiredtodescribetheintermediatestatesintime
system sizes, and fixed γ = 1, V = 5, Ω = 1.5, as a func- evolutionmethods.Forinstancein[55]D≈200wasrequired
tion of ∆. The inset shows the exponential decay of correlations, foradissipativeIsingchainoflengthN =40. In[3],theevo-
CzNz/2(n) = ⟨σNz/2σNz/2+n⟩−⟨σNz/2⟩⟨σNz/2+n⟩, for N = 50 in the lution required D of several hundreds [56], for a XXZ chain
variousregionsof∆.Right:Purityoftheconvergedsteadystatefor oflengthN =96,whenthesteadystatehasD=1.
thesamesystemsizes.Theinsetshowsthelocalpolarization,σz,for
n Ourapproachisbasedonthegroundstateoptimizationover
theN =10case.At∆=0theantiferromagneticorderingcanbeap-
preciated,whileforlargervaluesof∣∆∣thesteadystateapproaches MPS for a MPO Hamiltonian, and relies on the guaranteed
totalpolarization. existenceofavalid,positivesteadystate.Thisbasictechnique
is complemented with a warm-up phase or a suitable initial
guess, found to be crucial in practice for convergence to a
DissipativeIsingchain. Acomplementarykindofmodel physicalresultwithsmallbonddimension.
isonewheretheHamiltoniandynamicsinducescorrelations, When the steady state is degenerate, the simplest method
for instance an Ising chain, and the dissipation is purely lo- described in this paper might have problems to find a valid
cal. We consider a nearest neighbor Ising interaction, H = guess for the steady state. Specially in the situation of sev-
V ∑ σzσz +∑ (Ωσx− V−∆σz)+V(σz+σz ),andlocal eraldarkstates,thenullsubspaceoftheLindbladiancontains
4 i<N i i+1 i 2 i√ 2 i 4 1 N
dissipationgivenbyL = γσ+, i = 1,...N. Suchamodel infinitely many vectors which do not correspond to positive
i i
has attracted considerable attention in recent years, as it can operatorsandhencedonotconstitutevalidphysicalstates. 6
beeffectivelyrealizedinatomiclatticesystemsusingRydberg In principle it would be possible to complement the current
states[15, 53,54]. Inorder tocompare toexisting resultsin algorithmwithadditionaltechniques, suchassymmetries, in
the literature [53], we fix γ = 1, V = 5, Ω = 1.5. We com- order to reduce the degeneracy, or to construct a candidate
putethesteadystateforsystemsuptoN = 50andstudythe steadystatefromappropriatecombinationsofseverallinearly
squaredofthestaggeredmagnetization,M =∑ (−1)iσz/N, independentnullvectors,evenifnonpositive.
z i i
equivalent to the antiferromagnetic order parameter defined
in [53], and the purity of the steady state as a function of ∆
(Fig. 2). We find that our convergence criteria are met with SUPPLEMENTARY
smallbonddimensionD≤20. Ourresultsshowtheorderpa-
rametervanishingasN →∞,consistentwithshortrangecor- Basicsof“groundstate”searchingforMPO
relations, which indeed are observed to decay exponentially.
Weobservethatthepurityofthesteadystategrowsforlarge In the variational search for the best MPS approximation
absolutevaluesof∆. Thiscanbeeasilyunderstoodbygoing to the ground state of a Hamiltonian, H, the optimization
to the interaction picture with respect to the single body σz proceedsusinganalternatingleastsquares(ALS)strategy,in
i
terms in the Hamiltonian. This does not change the form of whichtheenergy,⟨Ψ∣H∣Ψ⟩/⟨Ψ∣Ψ⟩,issuccessivelyminimized
thedissipativeterms,butforverylarge∣V −∆∣,theσx terms withrespecttooneofthetensorsintheansatz, whilethere-
i
can be neglected in the rotating wave approximation. In this mainingtensorsarefixed[57]. Inthecaseofthesteadystate
situationthesingledarkstateofthedissipation,thefullypo- optimization,thebasicalgorithmisidentical,wheretheprod-
larized state ∣0⟩⊗N, is also an eigenstate of the Hamiltonian, uct Lˆ†Lˆ plays the role of the Hamiltonian. We use a MPO
andisthusasteadystate. description for the superoperator Lˆ(Eq. 2 in the main text),
The method can also be applied to other models, for in- whichcanbeeasilyconstructedasfortheHamiltoniancase,
stance with coherence induced by both the Hamiltonian and so that also Lˆ†Lˆhas a MPO representation with the squared
theenvironment[48]. bonddimensionoftheformer(seefigure3). Inparticular,in
Conclusion.— We have presented and analyzed a varia- thecasesdescribedinthiswork,theoperatorLˆcontainsonly
tionalalgorithmthatsearchsforaMPOapproximationtothe
steadystateofanopenquantumsystem. Thealgorithmisap-
plicabletoanymodelinwhichtheHamiltonianandtheLind-
bladoperatorscanbeexpressedasMPO.Insteadofsimulat- 6Noticethatthissituationcouldalsobeadversefortimeevolvingnumerical
ing the real time evolution of the system, as done by other methods.
5
localandnearest-neighborterms,buttheconstructioncanbe Physicalsolutions.
extended to more general situations [44]. Running the ALS
optimizationrequiresthecorrespondingeffectiveoperatoron Checking whether a given MPO corresponds to a physi-
eachtensor,i.e. thecontractionof⟨Ψ∣Lˆ†Lˆ∣Ψ⟩exceptforthis cal operator is in general a very hard problem. Neverthe-
tensor,whichcanbecomputedefficientlyusingthesameba- less, we can apply some compatibility checks to discard the
siccomputationalroutinesasthegroundstatealgorithms. mostunphysicalsolutions. Moreconcretely,forthespin-1/2
ThesearchforthegroundstateofLˆ†Lˆcanthenproducea chains studied, we demand that each of the local magnetiza-
MPSapproximationtothevector∣Φ(ρs)⟩,normalizedinEu- tions,σnx,σny,σnz,haveexpectationvalueswithinthephysical
clidean norm, which minimizes ∥Lˆ∣Φ(ρs)⟩∥. But not every range. 8 Wefoundthistobeenoughtoidentifythemostun-
such approximation will represent a possible physical state, physical intermediate solutions, and to obtain solutions with
asthoseneedtobepositiveoperators, withtraceone. Math- non-vanishing trace, so that the steady state could be prop-
ematically, there is a zero-eigenstate of Lˆ corresponding to erlynormalized,butitispossibletodefineadditionaltestsand
a positive operator, so that the algorithm will in general be performmoredemandingphysicalitychecks. Additionally,it
naturally looking for a physical solution (except in the case, mightalsobeinterestingtocomplementtheminimizationof
illustratedinthemaintext,ofseveraldarkstates). Neverthe- Lˆ†LˆwithaLagrangemultipliertermtofavorvectorsthatare
lessweidentifysomeadditionaltechniquesthathelpensuring notorthogonaltotheidentity.
thestabilityandfastconvergenceofthealgorithm.
Decidingconvergence.
Initialwarm-upphase.
Forthegroundstatesearch,convergenceoftheglobalalgo-
rithm can be decided in terms of the relative variation of the
Ingeneralitisconvenienttostartthevariationalsearchwith energyfromonevalueofthebonddimensiontothenext. In
a small bond dimension, D, and to use the result as initial ourcase,theexactsolutionhaszeroeigenvalue. Insteadofits
guess for the search with increased D until the required pre- relativevariation,wechecktheabsolutevalue,sowerequire
cision is attained. In the steady state case, nevertheless, if thatthefoundeigenvalueisbelowsomethreshold,andaddi-
theinitialbonddimensionisverysmall,itmighthappenthat tionally, demand convergence of some physical observables.
the algorithm at this stage converges to some vector not cor- Inparticular,werequirethatthelocalpolarizationsconverge,
respondingtoaphysicalstate,whichconstitutesabadinitial by constructing the N-component vectors, (σα,...σα), for
1 N
guess for the next rounds. This may lead to very slow con- α = x,y,z and requiring that the relative variation (in Eu-
vergence, and to an artificial increase of the required bond clidean norm) when increasing the bond dimension is below
dimensionforthefinalsolution. Hence,forthesmallestbond 10−4. Withthesecriteria,thesteadystateisconvergedforall
dimension,insteadofdirectlystartingthealgorithmonaran- thecasespresentedinthiswork. Additionallywecheckthat
dom vector, we implement a first warmup phase that con- otherobservables(suchastheantiferromagneticorderparam-
structs a more suitable initial state. This is then fed to the eter for the Ising model with local dissipation) and even the
usualDMRG-likesweep,untilconvergence. Wealsoobserve purity, which is far from being a local observable, are con-
thatthealgorithmbehavesbetterwhenthebonddimensionis vergedtoaprecisionbetterthan10−2.
increasedinsmallsteps.
Inparticular,forthe(closeto)reflectionsymmetricmodels
consideredhere,thefirstsweepsovertheMPStensorsarenot Isingchainwithcoherentdissipation.
performedfromlefttorightandback,asintheregularDMRG
algorithms,butsymmetricallyfromtheoutsidein,orfromthe In the main text we have illustrated the performance of
insideout. Wefindthatthisstrategyworksbestwhenstarting the algorithm in two typical examples, where the Hamilto-
fromtheseparableansatz, D = 1, forwhichthealgorithmis nian and the dissipation compete to induce or destroy cor-
extremely fast, and can thus be made to converge to numer- relations. Here we complement those results with a situa-
ical precision, and be repeated if necessary at a low cost. If tion in which both the unitary and the dissipative dynamics
afterthealgorithmhasconvergedforD=1,ifthesolutionis can induce coherence. Namely, we study an Ising model,
notcompatiblewithaphysicalstate,theprocedureisrepeated H =∑ σxσx +g∑ σz,withnearestneighborLindbladop-
i i i+1 i i
from a different initial state. Otherwise, the bond dimension erators,L =µσ++νσ− ,i<N (plusL =µσ+).Forchains
i i i+1 N N
isincreased,andtheobtainedstateisusedasinitialguess. 7 of up to N = 50 sites, we find convergence with remarkably
small bond dimension, D < 20. Figure 4 illustrates the ex-
pectation value of S2 = (∑ σz)2 in the steady state for the
z i i
7Whenexploringtheparameterspaceofacertainmodel,thistechniquecan
beusedforsomeparametervaluesandthentheconvergedresultscanbe 8Atleastwithinsometolerance,whichcanberelaxedforthesmallestbond
usedasguessesforthemodelwithslightlychangedparameter. dimensions.
6
Lˆ†LˆaregivenbythesquaredsingularvaluesofL.
Hereweillustratehowthegapsvaryforthetwomodelspre-
(a)TheMPOforρismappedtoaMPSfor∣Φ(ρ)⟩byChoi’s
sentedinthemaintex,inthecaseofsmallchains,usingexact
isomorphism.
diagonalization. In principle it would be possible to extend
this study and perform a finite size extrapolation, to predict
theperformanceofthealgorithmforincreasinglylargersys-
tems. Suchexhaustivestudyescapesthescopeofthispaper,
butthefinitesystemresultsalreadyallowustoappreciatethe
verydifferentsituationsthealgorithmcanface.
(b)FortheHamiltoniansandLiouvillianoperatorsconsidered
For the low dimensional Dicke model we observe that the
inthispaper,theLindbladmap,L(ρ)(left),is
correspondinglymappedtoaMPOactinglinearlyasLˆ∣Φ(ρ)⟩ gap closes, even for small systems, in the region of small
(right). g/γ. This is consistent with the parameter region on which
wefoundthemostdifficulties. However,noticethatthealgo-
rithmeasilyconvergedforsystemsizesN ≈10,inspiteofthe
closinggap.
For the Ising chain with local dissipation, when we fix
the parameters as in the text and vary only ∆, the gap ex-
hibitsmuchlessdramaticdependence. Nevertheless,wehave
shown that in this case, easily treatable by our method, the
(c)TheMPOsketched
steadystatecanpresentrichfeatures.
abovecanbeusedto
composeLˆ†Lˆ∣Φ(ρ)⟩.
0.2 N=3 0.8 N=3
FIG.3.MPO-MPSrepresentationofthestateandtheLindbladsuper NN==45 0.7 NN==45
operator. 0.15 N=6 0.6 N=6
0.5
0.1
0.4
particularcaseµ = 0.5andchangingvaluesofg andν, after 0.05 0.3
extrapolatingtheresultstoN →∞. Itiseasytocheckthatin 0.2
thethermodynamiclimittheidentityisalwaysasteadystate 00 0.5 1 1.5 2 0.−14 −2 0 2 4 6
g/γ ∆
for the symmetric case, µ = ν. This explains the vanishing
⟨S2⟩observedintheplotwhenν =0.5forallvaluesofg. FIG.5. GapofLˆ†Lˆasafunctionofthemodelparametersforsmall
z
chains, N = 3–6. ForthelowdimensionalDickemodel(left), the
gapclosesatsmallvaluesofg/γ. ThedissipativeIsingchain(right,
1 0.1 0.2 0.3
forγ=1,V =5,Ω=1.5)doesnotexhibitsuchastrongdependence
0.8 gg==01..50 1 ofthegaponthefreeparameter,∆. Inbothmodelsthefinalvalues
2/N0.6 ggg===125...500 0.8 ⟨Φ(ρ)∣Lˆ†Lˆ∣Φ(ρ)⟩thatwehaveacceptedareallsmallerthan10−5.
2iSz0.4 gg==1200..00 µ0.6
h WearethankfultoG.Giedke,D.Porras,H.Weimer,J.von
0.4
0.2 DelftandF.Mintertfordiscussions. J.C.issupportedbythe
0.2
0 “MPG CAS joint doctoral promotion programe(DPP)”, and
0.2 0.4 ν0.6 0.8 1 0.2 0.4 ν0.6 0.8 1 acknowledges Max-Planck-Institute of Quantum Optics, In-
stituteofPhysicsChineseAcademyofSciencesandFreiburg
FIG.4.TotalS2inthesteadystatefortheIsingmodelwithnearest-
z InstituteforAdvancedStudies,wherepartoftheresearchwas
neighbordissipation,asafunctionofν atfixedµ=0.5forvarious
carriedout. ThisworkwaspartiallysupportedbyEUthrough
valuesofg,andafterfinitesizeextrapolationfromsystemsizesup
toN =50.Therightplotshowsthesamequantityforvaryingµand SIQSgrant(FP7600645).
fixedg=1,forN =40.
[1] F.Verstraete,M.M.Wolf, andJ.IgnacioCirac,Nat.Phys.5,
SpectralgapsofLˆ†LˆandL 633(2009).
[2] S.Diehl,E.Rico,M.A.Baranov, andP.Zoller,Nat.Phys.7,
The gap of Lˆ†Lˆis determinant for the convergence of the 971(2011).
[3] Z.CaiandT.Barthel,Phys.Rev.Lett.111,150403(2013).
presented algorithm. It plays the same role as the energy
[4] M.Kastoryano, M.Wolf, andJ.Eisert,Phys.Rev.Lett.110,
gapinconventionalDMRG.Fortimeevolutionalgorithmsthe
110501(2013).
mostrelevantquantityisthespectralgapofL,whichcannot [5] S.Diehl,A.Micheli,A.Kantian,B.Kraus,H.P.Buchler, and
be related directly to the former. Instead, the eigenvalues of P.Zoller,NaturePhysics4,878(2008).
7
[6] M.Mu¨ller,S.Diehl,G.Pupillo, andP.Zoller,inAdvancesin [34] A.Molnar,N.Schuch,F.Verstraete, andJ.I.Cirac,Phys.Rev.
Atomic,Molecular,andOpticalPhysics,AdvancesInAtomic, B91,045138(2015).
Molecular, and Optical Physics, Vol. 61, edited by E. A. [35] T. Prosen and M. Zˇnidaricˇ, Journal of Statistical Mechanics:
PaulBermanandC.Lin(AcademicPress,2012)pp.1–80. TheoryandExperiment2009,P02035(2009).
[7] R.Dicke,Phys.Rev.93,99(1954). [36] L.Bonnes,D.Charrier, andA.M.La¨uchli,Phys.Rev.A90,
[8] H.J.Carmichael,JournalofPhysicsB:AtomicandMolecular 033612(2014).
Physics13,3551(1980). [37] J.Mendoza-Arenas,T.Grujic,D.Jaksch, andS.Clark,Phys.
[9] J. Eisert and T. Prosen, “Noise-driven quantum criticality,” Rev.B87,235130(2013).
arxiv:1012.5013(2010). [38] I.Pizˇorn,Phys.Rev.A88,043635(2013).
[10] B. Horstmann, J. I. Cirac, and G. Giedke, Phys. Rev. A 87, [39] F.W.G.Transchel,A.Milsted, andT.J.Osborne,“Amonte
012108(2013). carlo time-dependent variational principle,” arXiv:1411.5546
[11] M.Ho¨ning,M.Moos, andM.Fleischhauer,Phys.Rev.A86, [quant-ph](2014).
013606(2012). [40] A. H. Werner, D. Jaschke, P. Silvi, T. Calarco, J. Eisert, and
[12] L. M. Sieberer, S. D. Huber, E. Altman, and S. Diehl, Phys. S.Montangero,“Apositivetensornetworkapproachforsim-
Rev.Lett.110,195301(2013). ulating open quantum many-body systems,” arXiv:1412.5746
[13] T.ProsenandI.Pizˇorn,Phys.Rev.Lett.101,105701(2008). [quant-ph](2014).
[14] E. M. Kessler, G. Giedke, A. Imamoglu, S. F. Yelin, M. D. [41] D.Linzner,M.Honing, andM.Fleischhauer,“Bistabilityina
Lukin, andJ.I.Cirac,Phys.Rev.A86,012116(2012). dissipativerydberglatticemodel,”Privatecommunicationswith
[15] C. Ates, B. Olmos, J. P. Garrahan, and I. Lesanovsky, Phys. theauthors.
Rev.A85,043620(2012). [42] H.Weimer,Phys.Rev.Lett.114,040402(2015).
[16] T.Prosen,NewJournalofPhysics10,043026(2008). [43] D. Pe´rez-Garc´ıa, F. Verstraete, M. M. Wolf, and J. I. Cirac,
[17] P.Degenfeld-SchonburgandM.J.Hartmann,Phys.Rev.B89, QuantumInfo.Comput.7,401(2007).
245108(2014). [44] B.Pirvu,V.Murg,J.I.Cirac, andF.Verstraete,NewJournal
[18] A.C.Y.Li,F.Petruccione, andJ.Koch,Sci.Rep.4(2014). ofPhysics12,025012(2010).
[19] J.I.CiracandF.Verstraete,JournalofPhysicsA:Mathematical [45] M. Lubasch, J. I. Cirac, and M.-C. Ban˜uls, Phys. Rev. B 90,
andTheoretical42,504004(2009). 064425(2014).
[20] F.Verstraete,V.Murg, andJ.Cirac,AdvancesinPhysics57, [46] G.D.lasCuevas,N.Schuch,D.Pe´rez-Garc´ıa, andJ.I.Cirac,
143(2008). NewJournalofPhysics15,123021(2013).
[21] R.Oru´s,AnnalsofPhysics349,117 (2014). [47] S. R. White and D. J. Scalapino, Phys. Rev. Lett. 80, 1272
[22] S.White,Phys.Rev.Lett.69,2863(1992). (1998).
[23] U.Schollwo¨ck,AnnalsofPhysics326,96 (2011),january2011 [48] SeeSupplementalMaterial[url],whichincludesRefs.[57].
SpecialIssue. [49] M. Kliesch, D. Gross, and J. Eisert, Phys. Rev. Lett. 113,
[24] F.VerstraeteandJ.I.Cirac,Phys.Rev.B73,094423(2006). 160503(2014).
[25] M.B.Hastings,JournalofStatisticalMechanics: Theoryand [50] K.Baumann,C.Guerlin,F.Brennecke, andT.Esslinger,Na-
Experiment2007,P08024(2007). ture464,1301(2010).
[26] G.Vidal,Phys.Rev.Lett.93,040502(2004). [51] C. Hamner, C. Qu, Y. Zhang, J. Chang, M. Gong, C. Zhang,
[27] A.J.Daley,C.Kollath,U.Schollwo¨ck, andG.Vidal,Journal andP.Engels,NatureCommunications5(2014).
ofStatisticalMechanics:TheoryandExperiment2004,P04005 [52] M. P. Baden, K. J. Arnold, A. L. Grimsmo, S. Parkins, and
(2004). M.D.Barrett,Phys.Rev.Lett.113,020408(2014).
[28] M.Hartmann,J.Prior,S.Clark, andM.Plenio,Phys.Rev.Lett. [53] T. E. Lee, H. Ha¨ffner, and M. C. Cross, Phys. Rev. A 84,
102,057202(2009). 031402(2011).
[29] M. C. Ban˜uls, M. B. Hastings, F. Verstraete, and J. I. Cirac, [54] M.Ho¨ning,D.Muth,D.Petrosyan, andM.Fleischhauer,Phys.
Phys.Rev.Lett.102,240603(2009). Rev.A87,023401(2013).
[30] A.J.Daley,J.M.Taylor,S.Diehl,M.Baranov, andP.Zoller, [55] M. Ho¨ning, Many-body phases of open Rydberg systems and
Phys.Rev.Lett.102,040402(2009). signaturesoftopologyinquantumgases,Ph.D.thesis,Technis-
[31] F.Verstraete,G.-R.JuanJose´, andJ.I.Cirac,Phys.Rev.Lett. cheUniversita¨tKaiserslautern(2014).
93,207204(2004). [56] T.Barthel,privatecommunication.
[32] M.ZwolakandG.Vidal,Phys.Rev.Lett.93,207205(2004). [57] F. Verstraete, D. Porras, and J. Cirac, Phys. Rev. Lett. 93,
[33] M.Hastings,Phys.Rev.B73,085115(2006). 227205(2004).