Exact results on the quench dynamics of the entanglement entropy in the toric code
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] O c t Exact results on the quench dynamics of the entanglement entropy in the toric code
Armin Rahmani and Claudio Chamon
Physics Department, Boston University, Boston, MA 02215, USA (Dated: September 18, 2018)We study quantum quenches in the two-dimensional Kitaev toric code model and compute exactly the time-dependent entanglement entropy of the non-equilibrium wave-function evolving from a paramagnetic initialstate with the toric code Hamiltonian. We find that the area law survives at all times. Adding disorder to thetoric code couplings makes the entanglement entropy per unit boundary length saturate to disorder-independentvalues at long times and in the thermodynamic limit. There are order-one corrections to the area law fromthe corners in the subsystem boundary but the topological entropy remains zero at all times. We argue thatbreaking the integrability with a small magnetic field could change the area law to a volume scaling as expectedof thermalized states but is not sufficient for forming topological entanglement due to the presence of an excessenergy and a finite density of defects.
I. INTRODUCTION
Cold atomic systems have provided an ideal playground forthe experimental studies of the unitary evolution of thermallyisolated quantum systems.
Apart from the experimental in-terest, studies of quantum evolution have shed light on somefundamental questions such as thermalization.
A widely-studied scenario is the unitary evolution of aquantum state when the Hamiltonian is suddenly changed –a “quantum quench”, which has been the subject of manyrecent studies (see Refs. [12–15] for a few examples). A chal-lenging question, which has been addressed theoretically inone dimensional systems, is the evolution of the entanglemententropy following the quench.
More generally, such stud-ies are concerned with the fate of an out-of-equilibrium quan-tum state with a Hamiltonian that changes with some time-dependent protocol.The case where the coupling constants are changed acrossa quantum phase transition is of special interest. In recentexperiments with spinor condensates, for example, the forma-tion of ferromagnetic domains was observed in a quench to aHamiltonian with a symmetry-breaking ferromagnetic groundstate. Even though the interactions after the quench favor acertain type of order, it is not at all obvious that the ordershould emerge out of equilibrium. Since symmetry-breakingorder has a local order parameter, the non-equilibrium statecan support the formation of ordered domains. But what aboutquenches across a topological phase transition? Would any ofthe topological characteristics of the equilibrium phase showup in the non-equilibrium wave-function?Another interesting aspect of the non-equilibrium dynam-ics is the issue of thermalization and whether a large systemcan act as an effective heat bath for its subsystems. If thermal-ization occurs, due to the presence of some excess energy, thenon-equilibrium state at long times should be in some sense close to a thermal state at a finite temperature, which impliesvolume scaling for the entanglement entropy.The qualitative discussions above motivate a detailed anal-ysis of a quantum quench across a topological phase transitioninto the topological phase. We study the specific case of theKitaev toric code model with spin-polarized initial states. Arelated quantum quench problem with topologically orderedinitial states was recently studied numerically in Ref. [19]. We focus on the time evolution of the entanglement entropy be-cause it is intimately related to topological order through thetopological entanglement entropy. Entanglement entropyis also related to thermalization since, as mentioned above,volume scaling of the entropy as opposed to the area law isgenerically expected of thermal states.We compute exactly the bipartite entanglement entropy asa function of time. This is an example where exact calculationcan be done in 2+1D (thus far, analytical results for the time-evolution of entanglement had focused on 1+1D systems). Weshow that the entanglement entropy respects the area law atall times and if one adds disorder to the coupling constants,the recurrences or revivals are suppressed at long times andthe entanglement entropy per unit “area” (boundary length) P equals S/P = ln(4 /e ) in the thermodynamic limit when the initial polarization is inthe x or z direction.We find the dependence of this saturation value on the ini-tial state for an arbitrary direction of the initial magnetiza-tion. A new ingredient that comes in when considering theseinitial states is order-one corrections to the entanglement en-tropy from the convex and concave corners in the subsystemboundary. We find that the general form of the generated en-tanglement entropy is given by S ( t ) = f ( t ) P + f > ( t ) C > + f < ( t ) C < where P is the perimeter of the subsystem and C > ( C < ) thenumber of convex (concave) corners in its boundary. Thefunctions of time f , f > and f < are independent of the sub-system geometry. From the form of the entanglement entropyabove, we find that the topological entropy vanishes at alltimes.We study the effects of breaking the integrability in a per-turbative approach. We identify the mechanism by which anintegrability-breaking magnetic field can change the area lawto the volume law as expected of thermalized states. The phys-ical picture is that the area law survives at low orders in anexpansion in the integrability-breaking field. Going to higherorders however increases the thickness of a ring-like regionaround the subsystem boundary that contributes to the entan-glement entropy. Eventually at high enough orders, the ringwill collapse onto itself making volume scaling possible. Inthe perturbative treatment, we explicitly find that the topo-logical entanglement entropy remains zero at low orders. Weexpect that the excess energy and the finite density of defectsshould prevent the formation of any topological order in a sud-den quench even if one breaks the integrability.The paper is organized as follows. In section II we discussquantum quenches for spin-polarized initial states in the x or z direction. We calculate the entanglement entropy analyti-cally and discuss the effects of disorder in the coupling con-stants. In section III, we study the role of the sharp corners inthe subsystem boundary in the generation of entanglement byexactly calculating the second order Renyi entropy for an ar-bitrary direction of the initial magnetization. We then discussthe effects of breaking the integrability, which we treat pertur-batively, in section IV. We conclude the paper in section V. II. INTEGRABLE QUENCH: INITIAL MAGNETIZATIONIN THE x OR z DIRECTION
Consider the toric code Hamiltonian in a field H = − X p λ p B p − X s µ s A s − h X i ~n · ~σ i . (1)The Hilbert space of the above Hamiltonian is given by spinsliving on the bonds of a square lattice on a torus and the sumsare over all plaquettes, stars and bonds respectively. The pla-quette and star operators are defined as follows B p = Y i ∈ p σ zi , A s = Y i ∈ s σ xi where σ s are the Pauli matrices. For the usual Kitaev modelthe coupling constants λ and µ are independent of the plaque-ttes and the stars but here we keep the p and s indeces to allowfor disorder in the couplings. Note that for the model in zeromagnetic field ( h = 0 ), the ground state is independent of theactual values of the coupling constants λ and µ as long as theyare all positive.One can study the dynamics of this model by changing thecoefficients λ , µ and h in time with some protocol. In thissection, we consider a quantum quench where the system isinitially in a fully magnetized state and then evolves with theintegrable toric code Hamiltonian, namely, we are initially inthe ground state for h = 0 and λ p = µ s = 0 for all starsand plaquettes. We then turn off the field h at t = 0 (whichdoes not change the quantum state) and turn on the toric codepart of the Hamiltonian as follows: λ = λ ( t ) and µ = µ ( t ) ( λ (0) = µ (0) = 0 ).Consider for now, the case where ~n = ˆ x ( ~n · ~σ = σ x ).The ground state of the usual model (without disorder in thecouplings) exhibits a topological phase transition: for h > h c the system is a paramagnet and for h < h c , it is topologi-cally ordered. This topological phase is robust against weakdisorder in the couplings. Let us call the initial state | Ω i . Inthe case where the initial magnetization is in the x direction | Ω x i = | ... i in the | σ x σ x ... i basis. The initial state is a ten-sor product of single spin states and thus has no entanglement.In the paramagnetic initial state we have h σ x i = 1 for allspin operators. Before addressing the question of the entan-glement entropy evolution, let us consider the simple problemof what happens to the σ x expectation value. As a functionof time, we have h σ x ( t ) i = h Ω x | U † σ x U | Ω x i where U is theevolution operator given by U = exp i X p Z t λ p ( t ′ ) dt ′ B p + i X s Z t µ s ( t ′ ) dt ′ A s ! . (2)Since all star and plaquettes commute with one another, wehave dropped the time ordering operator. If we consider a sud-den quench, we have R t λ p ( t ′ ) dt ′ = λ p t and R t µ s ( t ′ ) dt ′ = µ s t . Let us focus on sudden quenches for a more compactnotation. At the end of the calculation we can replace λt and µt with the corresponding integrals over time-dependent cou-plings for a more general protocol.The expectation value of σ xpp ′ , a spin between two nearest-neighbor plaquettes p and p ′ is then given by h σ xpp ′ ( t ) i = h Ω x | e − itλ p B p e − itλ p ′ B p ′ σ xpp ′ e itλ p B p e itλ p ′ B p ′ | Ω x i . Since B p =
1, we have e ± itλ p B p = cos( λ p t ) ± i sin( λ p t ) B p . We have B p,p ′ σ xpp ′ = − σ xpp ′ B p,p ′ . Therefore we get h σ xpp ′ ( t ) i = cos (2 λ p t ) cos (2 λ p ′ t ) . In the completely clean case (a theoretical abstraction), theaverage of h σ xpp ′ ( t ) i over time is equal to but any amount ofdisorder will eventually lead to a state which unlike the initialparamagnet has a vanishing time-averaged expectation valuefor σ x .We next consider the main problem of this section, namelythe evolution of the entanglement entropy. The density matrixof the system can be written as ρ ( t ) = U | Ω ih Ω | U † where U is the time evolution operator and | Ω i is a generic initial state.By inserting the resolution of identity twice, we can write, ρ ( t ) = X α,σ h α | U | Ω ih Ω | U † | σ i| α ih σ | . (3)We would like to calculate the entanglement entropy be-tween two subsystems A and B . Any star or plaquette opera-tor acts on four spins which may all belong to A , all belong to B or some belong to A and some to B . In this last case, thestar or plaquette operator must lie at the boundary of the twosubsystems. We can then write U = U A U B U ∂ (4)where U A and U B only act on spins in A and B respec-tively and U ∂ only depends on star and plaquette operatorsat the boundary of the two subsystems. The three opera-tors U A , U B and U ∂ commute with one another. Each ba-sis state can be written as the tensor product of states forspins in A and B , | α i = | α A β B i . We now integrate out thespins in B to find the reduced density matrix for subsystem A , ρ A ( t ) = P β B h β B | ρ ( t ) | β B i and obtain, ρ A ( t ) = X α A ,σ A ,β B h α A β B | U | Ω ih Ω | U † | σ A β B i| α A ih σ A | . (5)Notice that since the trace of an operator is independent of thebasis, we can replace | β B i by U B | β B i in the above expres-sion. This amounts to calculating the trace in a rotated basis.We then substitute Eq. (4) into Eq. (5) above and obtain ρ A ( t ) = X α A ,σ A ,β B h α A β B | U A U ∂ | Ω ih Ω | U † A U † ∂ | σ A β B i| α A ih σ A | . The part of the Hamiltonian that only acts on B spins has noeffect on the time evolution of the reduced density matrix forsubsystem A . In order to calculate the entanglement entropybetween the two subsystems we need to find tr( ρ nA ) . We givea replica index to α , β and σ and drop the subsystem index forbrevity. States | α i i are in subsystem A and | β i i in B . Using h σ i | α i +1 i = δ σ i ,α i +1 we obtain, tr( ρ nA ) = X { α i } , { β i } n Y i =1 h α i β i | U A U ∂ | Ω ih Ω | U † A U † ∂ | α i +1 β i i (6)with α n +1 = α . Once again because traces are basis-independent we can replace | α i i by U A | α i i in the above equa-tion and we obtain an expression for tr( ρ nA ) that only dependson U ∂ and not the whole U . Let us define the following com-muting operators U ± p = e ± i P p ∈ ∂p λ p B p t , U ± s = e ± i P s ∈ ∂s µ s A s t (7)where ∂ p ( ∂ s ) is the set of all boundary plaquettes (stars)which act on both A and B spins. We then obtain tr( ρ nA ) = X { α i } , { β i } n Y i =1 h α i β i | U + p U + s | Ω ih Ω | U − s U − p | α i +1 β i i (8)with α n +1 = α . We now focus on the simple case whereinitially all spins are aligned in the x direction. In this case, | Ω i = | Ω x i is an eigenstate of U s and we can therefore sub-stitute U + s | Ω x ih Ω x | U − s = | Ω x ih Ω x | in Eq. (8) above. Also | α i and | β i are chosen in the σ x basis so the action of B p ona basis state is to flip the four spins around the plaquette p .For concreteness, we assume that the partition B contains allspins which lie inside a simply connected M × N rectangularregion. We then have M + N ) − boundary plaquettes.Let us consider the action of U ± p = Q p ∈ ∂ p e λ p B p t on a state | αβ i . Since B P =
1, we have U ± p = Y p ∈ ∂ p (cos λ p t ± iB p sin λ p t ) . (9)Notice that h αβ | U + p | Ω x i 6 = 0 only if the state the state | αβ i is obtained by acting with a certain number of ∂ p plaquette flips on | Ω x i . Let us consider then Ising spin variables livingat the center of ∂ p plaquettes. A plaquette is flipped for θ = − and not flipped for θ = 1 . Notice that, on a torus, flippinga certain set of plaquettes { θ } and flipping all other plaquettesthat do not belong to { θ } have the same effect on the spins,i.e. they describe the same element in the group of plaquetteflips. However, since we have eliminated U A and U B fromthe formulation of the problem by a basis rotation, we canuniquely specify a state with nonvanishing h αβ | U + p | Ω x i bya set of flipped boundary plaquettes. Let the state | α β i belabeled by { θ p } for p ∈ ∂ p . If the partitions are not too small,it is not possible for a set of boundary plaquette flips acting on | Ω x i to create a state with the same B spins and different A spins or vice versa. Therefore | α β i = | α β i = |{ θ p }i andwe obtain using Eq. (8) tr( ρ nA ) = X { θ p } ,p ∈ ∂ p (cid:0) h{ θ p }| U + p | Ω x i h Ω x | U − p |{ θ p }i (cid:1) n . (10)Using Eq. (9) we find that h{ θ p }| U + p | Ω x i is a product of termsthat are equal to cos λ p t for each unflipped plaquette with θ p = 1 and i sin λ p t for each flipped plaquette with θ p = − .We can then write h{ θ p }| U + p | Ω x i = Y p ∈ ∂ p (cid:18) θ p λ p t + i − θ p λ p t (cid:19) . (11)Switching the order of the sum and product after pluggingEq. (11) into Eq. (10), we obtain tr( ρ nA ) = Y p ∈ ∂ p X θ p (cid:18) θ p λ p t + 1 − θ p λ p t (cid:19) n . We have used the fact that for θ = ± , ± θ = (cid:0) ± θ (cid:1) .Summing over θ p = ± for each boundary plaquette sim-ply gives tr( ρ nA ) = Y p ∈ ∂ p (cid:2) (cos λ p t ) n + (sin λ p t ) n (cid:3) . (12)By analytic continuation of n , we write the entanglement en-tropy as S = − tr ( ρ A log( ρ A )) = lim n → ∂∂n tr( ρ nA ) and we finally obtain, S = − X p ∈ ∂ p (cid:2) cos λ p t ln (cid:0) cos λ p t (cid:1) + sin λ p t ln (cid:0) sin λ p t (cid:1)(cid:3) . (13)For the clean Kitaev model, the entanglement entropy is anoscillatory function of time whose maximum is proportionalto the number of boundary plaquettes or the perimeter of thesubsystem. For this initial condition, the star term in the toriccode Hamiltonian has no effect on the evolution of the entan-glement entropy because the initial state is an eigenstate ofthat term. Similarly if initially all spins are in the z direction,the plaquette term will have no effect. The completely cleansystem will never reach a steady state.If there is some disorder in the coupling constants however,the entanglement entropy per unit boundary length will satu-rate to a steady state value in the thermodynamic limit. Sup-pose the number of boundary plaquettes P is very large. As t → ∞ , due to the presence of some disorder in couplings λ p ,the phases λ p t are equally likely to take any value between and π modulo π . Therefore, taking the thermodynamiclimit first we can write S ( { φ i } ) = − P Pi =1 ζ ( φ i ) with ζ ( φ ) ≡ cos φ ln (cid:0) cos φ (cid:1) + sin φ ln (cid:0) sin φ (cid:1) where φ i is a random variable between and π drawn froma uniform distribution. The average entanglement entropy isthen given by ¯ S = − P π Z π dφ ζ ( φ ) = P ln(4 /e ) (14)We have ( S − ¯ S ) ∝ P which implies that for any realizationof disorder | S ( { φ i } ) P − ¯ SP | = O ( 1 √ P ) indicating that in the thermodynamic limit, the entanglementper unit area length saturates to a steady-state value ln(4 /e ) .Notice that these results, namely the survival of area law atall times and the oscillatory behavior in the entanglement en-tropy in the absence of disorder are generic features of inte-grable quenches where the integrals of motion are operatorswith local support. The entanglement entropy per unit bound-ary length relaxes in the thermodynamic limit to a steady statevalue in the presence of disorder in the coupling constants.Because the system can reach a steady state, it is natural toask if the steady state of a subsystem has the properties of athermalized state described by a Gibbs or a generalized Gibbsensemble. For such thermalized states however, we expectthe entanglement entropy of a small subsystem to scale as itsvolume in the limit of long times. The survival of the area lawat all times which follows trivially from having local integralsof motion indicates the absence of thermalization even to ageneralized Gibbs ensemble.
III. SUBSYSTEM CORNERS AND THE DYNAMICALENTANGLEMENT ENTROPY
Let us now consider a more generic case where the initialdensity matrix is not diagonal in either σ x or σ z basis and bothterms in the Kitaev Hamiltonian affect the entanglement dy-namics. We find that even though in the thermodynamic limitthe area law still applies, the physics of the out of equilibriumentanglement is much richer. We find that here there are order-one corrections to the entanglement entropy from the concaveand convex corners of the subsystem.Subleading corrections to the entanglement entropy aris-ing form the sharp corners in the boundary were found be-fore for static entanglement entropy in the logarithmic correc-tions to the area law in a class of conformal quantum critical FIG. 1. Subsystem B consists of spins inside the painted region,while subsystem A includes spins outside and on the boundary ofthis region. The convex (concave) corners for subsystem B are rep-resented by squares (circles). points. Here we find that corners can also make sublead-ing contributions to the dynamical entanglement. We find thatwhen subsystem B consists of spins inside a region with C > convex and C < concave corners as seen in Fig. 1, the entan-glement entropy after a quantum quench behaves as S = f P + f > C > + f < C < (15)where P is the perimeter of the subsystem and f , f > and f < are functions of time that do not depend on the subsystemgeometry. For a subsystem with a mostly smooth boundaryand only a few sharp corners, the leading term is the f P butfor a subsystem with a corrugated boundary, the entanglemententropy can be different in the thermodynamic limit.Notice that the number of concave and convex corners arenot independent. Suppose subsystem B consists of spins in-side of an area with m connected regions, each with a cer-tain number of handles such that the total number handles is g . As an example in Fig. 1, we have C > = C < = 11 and m = g = 2 . We argue that C > = C < + 4 ( m − g ) . This is because for a rectangular simply connected region, wehave C > = 4 and C < = 0 and transforming the shape byadding or subtracting plaquettes on the square lattice does notchange C > − C < . Therefore each simply connected regionmakes a contribution of +4 to C > − C < and by a similarargument, each hole contributes − . A. Initial magnetization in the y direction To see how the contribution from sharp corners comes in,consider an initial state with all the spins in the positive y di-rection. Let us work in the the σ y basis with σ y |±i = ±|±i , σ x |±i = |∓i and σ z |±i = ∓ i |∓i . The action of a star orplaquette operator in this basis is then to flip the correspond-ing four σ y spins modulo a phase factor of ± from B p . Wewill see later that this sign has no effect on the entanglemententropy.We assume for simplicity that there is a minimum distanceof times the lattice spacing between all corners. Consider aconvex corner as in the right hand side of Fig. 2. The β spinslive in the dashed blue (color online) bonds and the α spins FIG. 2. A convex (concave) corner of the B subsystem is shown onthe right (left) hand side. The dashed blue bonds belong to subsystem B and the solid black bonds to subsystem A . A corner boundaryplaquette (star) is represented by a red square (cross). on the solid black bonds. Flipping the corner plaquette repre-sented by a red square has the same effect on β spins as flip-ping the two boundary stars represented by red crosses. Theseoperations of course have a different effect on α spins. Be-cause we assumed a minimum distance between the corners,these corner star and plaquette operators uniquely correspondto one corner. Similarly as in the left hand side of Fig. 2, flip-ping the corner star at a concave corner has the same effect onthe α spins as flipping the two corresponding plaquettes. Wesee from this discussion that the effects of stars and plaquetteflips mix at the corners, which leads to a different entangle-ment generation than along a smooth boundary.The details of the derivation will be given in appendix A.Let us write for now the result for the entanglement entropyfor a subsystem B consisting of the bonds inside an N × M rectangular region. A boundary star or plaquette that is notlocated at a corner is refereed to as a regular star or plaquette.As seen in Fig. 2, at each convex corner, one plaquette andtwo stars are special. For our N × M subsystem, we have fourconvex corners, M + N − regular boundary plaquettes and M + N − regular boundary stars and the entanglemententropy for a clean system after a sudden quench is given by S = 2( M + N − s p + 2( M + N − s s + 4 s c (16)where the contributions from regular plaquettes ( s p ), regularstars ( s s ) and convex corners ( s c ) are given by s p = (sin λt ln sin λt + cos λt ln cos λt ) s s = (sin µt ln sin µt + cos µt ln cos µt ) s c = [ r ln r + 2 r ln r + r ln r ] with r = sin λt sin µt + cos λt cos µtr = sin µt cos µtr = sin λt cos µt + cos λt sin µt. The above solution conforms to the generic form Eq. (15) for P = 2( M + N − , C > = 4 , C < = 0 , f = s p + s s and f > = s c − s s − s p .In appendix A, we give a derivation of the above result inthe more general setting of an arbitrary subsystem geometry.If instead of a sudden quench we turn on the coupling con-stants with some protocol λ ( t ) and µ ( t ) , we just need to re-place λt and µt by their integrals over time. ℓ ℓ I ℓ IIIII IV
FIG. 3. The partitions used for calculating the topological entangle-ment entropy.
Using the partitions of Ref. 21, shown in Fig. 3, the topo-logical entanglement entropy is given by S topo = S I − S II − S III + S IV . (17)From the form of Eq. (15), we can then see that although thereare order-one corrections to the area law, the topological en-tanglement entropy remains exactly zero at all times indepen-dently of how the toric code Hamiltonian is turned on.It seems plausible that for any initial state, as long as theintegrability is not broken (with a field h for example), evo-lution with the Kitaev Hamiltonian does not change the topo-logical entanglement entropy for any time dependence of thecoupling constants. This is because this integrable evolutiondoes not allow propagation of information. B. Arbitrary initial magnetization
The von Neumann entropy is a special case (limit of n → )of the order n Renyi entropy defined as S n ( ρ A ) = 11 − n tr( ρ nA ) . Hereafter we refer to the second order ( n = 2 ) Renyi entropyas the Renyi entropy. The Renyi entropy S ( ρ A ) = − tr( ρ A ) is simpler to calculate than the von Neumann entropy and isan equally valid measure of entanglement. It has been shownthat the topological entanglement entropy calculated by anyorder Renyi entropy contains the same information in the abe-lean topological phases such as the toric code. For an in-tegrable quench from a paramagnetic state with magnetiza-tion in an arbitrary direction, we exactly calculate the Renyientropy S ( ρ A ) . We consider the clean case with position-independent coupling constants. With disorder in the cou-plings, the entropy saturates to the time-average of the cleancase similarly to the simpler initial conditions considered be-fore.In the σ x basis, an arbitrary paramagnetic state can be writ-ten as | Ω i = (cid:0) c + | + i + c − |−i (cid:1) ⊗ N . The calculation is rather involved and is done by a transfermatrix method explained in appendix B. We find the samegeneral features for an arbitrary initial magnetization as forthe initial magnetization in the y direction. Let us summarizethese features below. FIG. 4. The dependence of the saturation value of the generatedRenyi entropy per unit boundary length on the two angles parame-terizing the magnetization direction of the initial state.
The Renyi entropy can be written as a sum of contributionsfrom the perimeter of the subsystem as well as the convex andconcave corners. These contributions are oscillatory functionsof time in the absence of disorder in the couplings which atlong times, saturate to the values given by the time-averageof the clean case if we have some disorder in the couplings.The saturation values are independent of the final couplingconstants and the protocol used for turning them on. They dohowever depend on the initial state.For a large subsystem with only a few sharp corners, theleading term for the Renyi entropy is the first term in Eq. (15),namely S ≈ f P . The saturation value of the Renyi entropyper unit boundary length is plotted in Fig. 4 as a function ofthe Θ / and Φ angles which determine the initial state through c + = cos Θ / e i Φ and c − = sin Θ / . Notice that this satura-tion value is a unique function of h B p i = (sin Θ cos Φ) and h A s i = (cos Θ) . IV. BREAKING THE INTEGRABILITY:ENTANGLEMENT ENTROPY IN THE INTERACTIONPICTURE
To study the effects of breaking the integrability, we fo-cus on one simple case. Namely, we start with a state thathas all spins aligned in the x direction. As opposed to theintegrable case where the magnetic field is turned off beforeturning on the integrable toric code Hamiltonian, the the toriccode Hamiltonian is turned on in the presence of a small mag-netic field h . We cannot perform exact analytical calculationsin this case. We can however discuss the effects of the mag-netic field by calculating a formal expansion for the (Renyi)entanglement and topological entropies in powers of h . Theexpansion can only be considered as a perturbative treatmentif the toric code coupling constants are quickly pushed above h and the time scales considered are not too large.Two main insights are derived from this analysis. First we find indications as to how breaking the integrability can even-tually lead to volume scaling in the entanglement entropy. Thephysical picture is that in the expansion in powers of h , thecontribution to the entanglement entropy comes from ring-shaped areas around the subsystem boundary whose thick-ness grows at higher and higher orders in h . Eventually, ata high enough order, these correlated regions meet and theterms in the expansion no longer respect the area law. Evenfor a small h , these high order terms become important atlong times making volume scaling of the entanglement en-tropy possible as expected of thermalized states. Secondly,we find that in this expansion in powers of h , the topologi-cal entanglement entropy remains zero at low orders until wereach an order proportional to the system size, suggesting thata small integrability-breaking term is not enough for generat-ing topological entanglement.A useful notion that we introduce in this section is the en-tanglement entropy in the interaction picture. We are used tothe Shcr¨odinger, Heisenberg and interaction pictures for cal-culating the expectation values of observables, which are ofcourse independent of the picture. If one defines the entangle-ment entropy as a property of the wave function alone, thensince the wave function does not change in the Heisenbergpicture, calculating the entanglement entropy in that pictureleads to the nonsensical result that time evolution does notchange the entanglement of a quantum state. An interestingalternative is to try and define the entanglement entropy en-tirely in terms of the correlation functions of appropriate op-erators. We do not pursue this path here. We argue insteadthat as explained below, the entanglement entropy calculatedonly from the wave function in the interaction picture can bea useful concept.Let us first consider a non-interacting Hamiltonian H suchas a paramagnet. Consider a partitioning of the system intosubsystems A and B . The evolution operator for H can bewritten as U ( t ) = U A ( t ) U B ( t ) where U A ( U B ) only acts on the degrees of freedom in A ( B ).Notice that acting with U ( t ) on any wave function does notchange its entanglement properties. Now consider a Hamilto-nian H = H + V with the corresponding evolution operator U ( t ) .The wave functions at time t in the Schr¨odinger and inter-action pictures are respectively given by | ψ S ( t ) i = U ( t ) | ψ (0) i , | ψ I ( t ) i = U † ( t ) U ( t ) | ψ (0) i . Now since acting with U or U † on any wave function does notchange its entanglement, the interaction picture wave functionhas the same entanglement entropy as the Schr¨odinger picture.Note however the calculation may be much simpler in the in-teraction picture.In the discussion above, we made use of the interaction pic-ture entanglement entropy because we had a Hamiltonian H whose evolution operator did not change the entanglement en-tropy. At the next level of complexity, where using the in-teraction picture wave function may be useful for calculatingthe entanglement entropy, the evolution operator of the unper-turbed Hamiltonian ( U ) may have another known property.For example it may only create entanglement entropy scalingwith the subsystem area or it could be incapable of changingthe topological entanglement. This is the case for the inte-grable toric code Hamiltonian. Now if we are interested infinding whether a perturbation can generate entanglement en-tropy with volume scaling or whether it can generate topolog-ical entanglement entropy, we can perform the calculation forthe wave function | ψ I ( t ) i in the interaction picture.Let us go back to our quench problem and write the wavefunction in the interaction picture with H ( t ) = − X p λ ( t ) B p − X s µ ( t ) A s and V = h X i σ xi . Starting with the initial state | i ≡ | Ω x i with all the spinsin the positive x direction, the wave-function at time t in theinteraction picture is given by | ψ I ( t ) i = e i R t dt ′ H ( t ′ ) | ψ S ( t ) i = ˆ U V ( t ) | i (18)where ˆ U V ( t ) is given by ˆ U V ( t ) = 1 − ih Z t dt ˆ V ( t )+ ( − ih ) Z t dt Z t dt ˆ V ( t ) ˆ V ( t ) + · · · for ˆ V ( t ) = h e − i R t dt ′ λ ( t ′ ) P p B p X i σ xi ! e i R t dt ′ λ ( t ′ ) P p B p . For the simplicity of notation, we focus on the sudden quenchand write R t dt ′ λ ( t ′ ) as λt , bearing in mind however that fora more general protocol, we can substitute the integral at theend of the calculation.We can explicitly calculate ˆ V ( t ) and the result is ˆ V ( t ) = cos 4 λt + 12 C + sin 4 λt P + cos 4 λt − D ≡ c ( t ) C + p ( t ) P + d ( t ) D (19)with the operators C , P and D defined as follows C = X i σ xi P = X p ( zz zy + zz yz + yz zz + zy zz ) D = X d zx zzzzz . (20) The sums are over all bonds, plaquettes and dominos andthe graphical notation represents a product of Pauli matrices.For example if the bonds in a plaquette are numbered clock-wise from to with the top bond being number we have zz zy = σ y σ z σ z σ z , zz yz = σ z σ y σ z σ z , · · · The closed form in Eq. (19) is simply derived from the re-peated substitution of the following commutation relations "X p B p , C = 2 iP, "X p B p , P = − i ( C + D ) "X p B p , C + D = 4 iP in the Baker-Hausdorff expansion of ˆ V ˆ V ( t ) = h X i σ xi + ( − iλt ) "X p B p , X i σ xi + ( − iλt ) "X p B p , "X p B p , X i σ xi + · · · ! and resumming the series.The Renyi entropy for the interaction picture wave function | ψ I ( t ) i is given by S ( t ) = − log h α β | ˆ U V | i h | ˆ U † V | α β i ×h α β | ˆ U V | i h | ˆ U † V | α β i (21)where summation over repeated indices is implied and as be-fore, the α and β spins belong to the A and B subsystemsrespectively. First let us consider the first order term in the tr( ρ A ) (the expression for the Renyi entropy without the log).We use the shorthand notation V V · · · V | {z } n times = Z t dt · · · Z t n − dt n ˆ V ( t ) · · · ˆ V ( t n ) , and similarly for the Hermitian conjugate. The first order termvanishes as seen below tr( ρ A ) = 1 − ih h | V | i + 2 ih h | V † | i + · · · = 1 + O ( h ) . The second order contribution to tr( ρ A ) is given by tr( ρ A ) (cid:12)(cid:12) o ( h ) = h (cid:2) |h β | V | i| + 2 |h α | V | i| − h | V | i − h | V † | i − h | V V | i − h | V † V † | i (cid:3) . Let us assume that the whole system is on a torus and thetotal number of bonds (spins) is N b . The number of plaquettes N p = N b and the number of dominos in the system is N d = N b . Let us introduce another shorthand notation for integralsinvolving the functions c , p and d , ( f , f , · · · , f n ) = Z t dt · · · Z t n − dt n f ( t ) · · · f n ( t n ) so we have for example ( c ) = R t dt c ( t ) and ( c, p ) = R t dt R t dt c ( t ) p ( t ) . We denote the number of plaque-ttes (dominos) with all bonds in subsystem A by n pA ( n dA )and similarly by n pB ( n dB ) for subsystem B . We then have |h β | V | i| = N b ( c ) + n pB ×
16 ( p ) + n dB ( d ) |h α | V | i| = N b ( c ) + n pA ×
16 ( p ) + n dA ( d ) h | V | i = N b ( c ) h | V V | i = N b ( c, c ) + N p ×
16 ( p, p ) + N d ( d, d ) . The factors of
16 = − i × i come from the fact that eachterm in P is a sum of four operators that flip the spins arounda plaquette. Noting that because c ( t ) c ( t ) is obviously sym-metric under the exchange t ↔ t , we have ( c, c ) = ( c ) / and similarly for ( p, p ) and ( d, d ) . We then obtain tr( ρ A ) (cid:12)(cid:12) = 1 − h ( N p − n pA − n pB )( p ) − h ( N d − n dA − n dB )( d ) , which using log(1 + x ) = x + O ( x ) gives S = h ( N p − n pA − n pB )( p ) (22) + h ( N d − n dA − n dB )( d ) + O ( h ) . This expression demonstrates the physical picture for buildingvolume entanglement which we alluded to before. We observefrom the above expression that at this order the contribution tothe entropy is from a ring-shaped area of thickness two latticespacings around the subsystem boundary. At higher orders,the thickness of this correlated region, which contributes tothe entanglement entropy, grows until the ring closes in onitself, making volume scaling possible.Notice that in contrast to exact area law which is possiblein a finite system, since the Renyi entropy is symmetric underexchanging the subsystems, volume scaling can only hold ex-actly for a finite subsystem in an infinite system and as a goodapproximation in the limit when a subsystem is much smallerthan the whole system.It is a simple exercise in geometry to count n p and n d ofthe two subsystems for partitions shown in Fig. 3 and verifydirectly that no topological entanglement is generated at order h . The details of this calculation will be given in appendix C.We will also show in that Appendix that at the next order witha non-vanishing contribution to the Renyi entropy ( O ( h ) ),the topological term will still vanish due to an exact cancella-tion.Similarly to the volume scaling, topological entanglementis non-perturbative and does not show up in an expansion inpowers of h in the thermodynamic limit. For a finite system,we expect non-vanishing contributions to the topological en-tropy to appear at orders that scale with the system size. Let uscomment here that on physical grounds, in the case of a sud-den quench, we do not expect any topological entanglementto form even at infinitely long times. A trend toward formingvolume entanglement observed at low orders suggests that inthe limit of large times, the generic thermalization behavior isexpected. Although we cannot give a mathematical proof, thefragility of topological order to thermal fluctuations and the extensive excess energy in the system suggests that thenon-equilibrium state, after a sudden quantum quench, willlikely bear no signature of the equilibrium topological orderof the phase.We emphasize that the expansion used here cannot betreated as a perturbative treatment in the adiabatic limit wherethe toric code coupling constants are turned on slowly becausein that case λ < h for short times. We know that in the caseof adiabatic evolution, we approximately remain in equilib-rium, i.e. in the ground state manifold, and topological orderemerges after a time of order system size. V. CONCLUSIONS
We studied the quench dynamics of the entanglement en-tropy in the Kitaev toric code model. In the case of the inte-grable quench we found the entropy exactly as a function oftime. The leading order contribution was found to scale asthe subsystem perimeter with an oscillatory behavior that getssuppressed in the presence of some disorder in the couplingconstants leading to the saturation of the entanglement en-tropy to disorder-independent values respecting the area law.We obtained results for the generated entanglement entropyas a function of an initial paramagnetic state with arbitrarymagnetization. We studied the effect of the sharp corners inthe subsystem boundary and found that in analogy with thestatic entanglement discussed in the literature, they make sub-leading contributions to the dynamically generated entangle-ment entropy as well. We showed that the topological entan-glement entropy remains zero in this quantum quench. Theresults of the integrable quench as well as the analysis of theintegrability-breaking perturbations strongly suggest that thenon-equilibrium state after a quench into the topological phasedoes not develop any of the topological structure present in theequilibrium phase.We used the concept of the entanglement entropy in the in-teraction picture which simplifies perturbative approaches tothe calculation of the entanglement entropy. This allows usto tease apart in this particular problem the effects of the timeevolution with an unperturbed Hamiltonian that cannot gen-erate either entanglement with volume scaling or topologicalentanglement and focus only the part of the evolution that canpossibly lead to volume scaling or topological entanglement.By studying the effects of an integrability-breaking perturba-tion, we found a mechanism where volume scaling of the gen-erated entanglement entropy can emerge at high orders in per-turbation theory as expected of thermalized states. The notionof thermalization after a quench naturally ties in with topolog-ical order through the fragility of topological order to thermalfluctuations.
ACKNOWLEDGMENTS
We are grateful to C. Castelnovo, A. Hamma, M. Haque,A. Polkovnikov and M. Rigol for helpful discussions. Thiswork was supported by in part by the DOE Grant DE-FG02-06ER46316.
Appendix A: Derivation of entropy for y direction initialmagnetization We explicitly calculate in this appendix the contributionsto the entanglement entropy from regular boundary stars andplaquettes as well as the convex and concave corners when theinitial magnetization is in the y direction.First let us find the number of regular boundary stars ( S r )and plaquettes ( P r ) in terms of the subsystem perimeter andthe number of convex and concave corners. We define theperimeter P as the total number of boundary plaquettes orboundary stars (they are equal). We then have P r = P − C > − C < , S r = P − C > − C < . Once again we use Eq. (8) which generally holds for anyinitial state | Ω i . The only intermediate states (in the σ y ba-sis) which contribute to the tr( ρ nA ) are the ones connected to | Ω y i by some star or plaquette flips. We label these states byvariable θ = ± and ψ = ± living on the boundary plaque-ttes and stars respectively. A negative θ ( ψ ) indicates that thecorresponding plaquettes (star) is flipped with respect to | Ω i y . | αβ i = |{ θ p r , ψ s r , θ p > , ψ jp > , ψ s < , θ js < }i (A1)where p r , s r , p > and s < label regular plaquettes, regularstars, convex corner plaquettes and concave corner stars re-spectively. For j = 1 , , { ψ jp > } ( { θ js < } ) are variables liv-ing on the two stars (plaquettes) corresponding to the convex(concave) corners. Notice unlike regular star variables the twostars at a convex corner are labeled by the corner plaquette andsimilarly the two plaquettes at a concave corner are labeled bya corner star since a convex (concave) corner uniquely corre-sponds to a corner plaquette (star).Now suppose two such states have the same β spins. Allregular and concave corner plaquette and star flips must bethe same for these two states. As seen in Fig. 2 however, at aconvex corner the θ or ψ variables for these two states couldbe different if both the convex corner plaquette and the corre-sponding two corner stars simultaneously have opposite flips.Then for two states | αβ i and | α ′ β i (which can be labeled asEq. ( A1) because they contribute tr( ρ nA ) ) all their θ and ψ variables must be the same except for the ones correspond-ing to convex corners, namely { θ p > , ψ jp > } and { θ p > , ψ jp > } ′ which follow the relation { θ p > , ψ jp > } ′ = { γ p > θ p > , γ p > ψ jp > } with auxiliary variable γ p > = ± . Similarly we can introducevariables ρ s < = ± for concave corner and for two states | αβ i and | αβ ′ i we must have { ψ s < , θ js < } ′ = { ρ s < ψ s < , ρ s < θ js < } . Introducing a replica index l on γ l and ρ l , we can write | α β i = |{ θ p r , ψ s r , θ p > , ψ jp > , ψ s o , θ js o }i| α β i = |{ θ p r , ψ s r , γ p > θ p > , γ p > ψ jp > , ψ s < , θ js < }i| α β i = |{ θ p r , ψ s r , γ p > θ p > , γ p > ψ jp > , ρ s < ψ s < , ρ s < θ js < }i| α β i = |{ θ p r , ψ s r , γ p > θ p > , γ p > ψ jp > , ρ s < ψ s < , ρ s < θ js < }i· · · Let us for the moment focus on the Renyi entropy − tr( ρ A ) for simplicity. The calculation for tr( ρ nA ) will be very similar.First we can see from the form of the the states above that thesign factors from B p acting on states in the σ y basis have noeffect. A B p plaquette flip can only give a different sign forthe matrix elements involving | α β i and | α β i at a corner.But then the plaquette flip at that corner would also give adifferent sign for the matrix elements involving | α β i and | α β i and the overall sign for the product of the four matrixelements appearing in tr( ρ A ) will be positive.The calculation is then similar to the simpler cases with theinitial magnetization in the x and z directions except for eachconvex or concave corner there is an additional auxiliary vari-able to sum over. Again, each matrix element can be writtenas a product similar to Eq. (11). By simply bringing togetherterms we can write tr( ρ A ) in the following form: tr( ρ A ) = X Y p r f p r Y s r f s r Y p > f p > Y s < f s < (A2)where the summation is over all θ , ψ , γ and ρ variables forboth replicas. We can take the sum Eq. A2 inside the productsby switching their order and perform the sum for individualterms. This gives the factorized form below tr( ρ A ) = Y p r F p r Y s r F s r Y p > F p > Y s < F s < . (A3)where we will compute the F s below. Assuming no disorderand taking the log of both sides, we immediately obtain ourmain result for the entropy S = ( P − C > − C < ) ln F p r + ( P − C > − C < ) ln F s r + C > ln F p > + C < ln F s < which conforms to the general expression Eq. (15).For regular stars and plaquttes, the terms appearing in theabove Eq. (A2) are identical to what we obtained when theinitial spins were in the x or z direction. For example f p r = (cid:18) θ p r λ p r t + 1 − θ p r λ p r t (cid:19) . Let us now define the following function for a more com-pact notation g ( ϑ, x ) ≡ ϑ x + 1 − ϑ x. We can then write f p r = [ g ( θ p r , λ p r t )] , f s r = [ g ( ψ s r , µ s r t )] . (A4)0At the corners, we have a more complicated structure andfor example, f p > , in addition to θ p > , also depends on γ p > and ψ jp > for j = 1 , . Notice that for the Renyi entropy α = α and therefore γ p > = 1 . We can explicitly write f p > = g ( θ p > , λ p > t ) g ( γ p > θ p > , λ p > t ) × Y j g ( ψ jp > , µ jp > t ) g ( γ p > ψ jp > , µ jp > t ) . It can be easily shown that g ( ϑ, x ) = | sin x cos x | e ϑ ln | cot x | . Using the above identity we can write f p > = ( 18 sin 2 λ p > t sin 2 µ p > t sin 2 µ p > t ) × e (1+ γ p> ) ( θ p> ln | cot λ p> t | + P j ψ jp> ln | cot µ jp> t | ) . A similar expression with θ ↔ ψ , λ ↔ µ and γ → ρ gives f s < . We can now perform the summations and we obtain F p r = X θ pr f p r = (cos λ p r t ) n + (sin λ p r t ) n (A5)with n = 2 for the Renyi entropy and similarly for F s r with λ ↔ µ . To calculate F p > = X θ p> ,ψ p> ,ψ p> ,γ p> f p > , we first perform the sum over θ s and ψ s and then over γ s andwe obtain after some algebra F p > = X k =1 (cid:0) ξ k sin λ p > t sin µ p > t sin µ p > t (A6) + ξ k cos λ p > t cos µ p > t cos µ p > t (cid:1) n with n = 2 for the Renyi entropy and ξ = 1 , ξ = tan λ p > t , ξ = tan µ p > t and ξ = tan µ p > t . We can similarly find F s < for concave corners by switching λ ↔ µ . The general-ization to n > is similar and leads to the appearance of adifferent n in Eq. (A6) and Eq. (A5). Appendix B: Transfer matrix method for an arbitrary initialmagnetization
In this appendix, we discuss the transfer matrix methodused in calculating the Renyi entropy for an arbitrary initialmagnetization in the positive ˆ n direction. The initial state canbe written as | Ω i = (cid:0) c + | + i + c − |−i (cid:1) ⊗ N where | + i and |−i are the eigenstates of σ x . Let us begin by writing an explicit equation for tr( ρ A ) . Us-ing n = 2 in Eq. (8) we can write tr( ρ A ) = X α,β h α β | U + p U + s | Ω ih Ω | U − s U − p | α β i (B1) h α β | U + p U + s | Ω ih Ω | U − s U − p | α β i . Consider the first of the four matrix elements appeasing in theequation above, namely h α β | U + p U + s | Ω i . We represent the | α β i spin configuration in a basis where boundary spins, i.e.spins belonging to a boundary star or plaquettes, are in the σ x basis and all other spins in the σ ~n basis. In a schematic no-tation explained below, we can then write the matrix elementas h α β | U + p U + s | Ω i = X { θ } h e iµt P α ′ βαα (B2) (cos λt ) p ↑ ( i sin λt ) p ↓ c + b ↑ c − b ↓ i where in cos λt and sin λt we can more generally make thesubstitution λt → R t λ ( t ′ ) dt ′ .The notation in the above equation needs some explana-tion. As before we have introduced Ising variables θ p on the P boundary plaquettes. For a given configuration of thesevariables { θ } , p ↓ and p ↑ are the number of flipped and un-flipped boundary plaquettes respectively and b ↓ and b ↑ are thenumber of up and down boundary spins in the state obtainedfrom | α β i after flipping the plaquettes represented by { θ } .We have introduced a new notation where α denotes only theboundary spins in subsystem A that belong to a boundary pla-quette and the other boundary spins in A, which only belong toa boundary star, are denoted by α ′ . Similarly boundary spinsin subsystem B that belong to a boundary star are labeled β and the rest β ′ as seen in Fig. 5. Assuming that the total num-ber of boundary spins is B , we can relate p ↑ , ↓ and b ↑ , ↓ to thespins variables above as follows p ↑ + p ↓ = P, p ↑ − p ↓ = X θ and b ↑ + b ↓ = Bb ↑ − b ↓ = X α ′ + X θα + X θβ ′ + X θβθ. In the sums above each α and β ′ spin is multiplied bythe Ising variable θ of the boundary plaquette it belongsto. Each β spin belongs to two boundary plaquettes and istherefore multiplied by the two corresponding θ variables.In the notation explained above, we now write the term h α β | U + p U + s | Ω i as follows h α β | U + p U + s | Ω i = ( c + c − ) B ( i sin λt cos λt ) P e iµt P α ′ βαα × A P α ′ X { θ } A ( P θα + P θβ ′ + P θβθ ) H P θ where H ≡ − i cot λt, A ≡ c + c − . (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1) (cid:0)(cid:0)(cid:1)(cid:1) (cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1) α β β ′ α ′ FIG. 5. The notation used in Eq. B2. The boundary spins in subsys-tem A are divided into α ′ spins that do not belong to any boundaryplaquette and α spins that do belong to a boundary plaquette. Sim-ilarly, the boundary spins in subsystem B are decided into β ′ spinsthat do not belong to any boundary star and β spins that do belong toa boundary star. We have represented the α ′ spins by a black circleon a dashed black bond, the α spins by a black circle on a solid blackbond, the β ′ spins by a blue square on a dashed blue bond and the β spins by a blue square on a solid blue bond. Like the one-dimensional classical Ising model in a mag-netic field, we can now integrate out the θ Ising variables usingthe transfer matrix method. The terms containing θ are pro-moted to × transfer matrices where each transfer matrix T × is between two adjacent boundary plaquettes and h α β | U + p U + s | Ω i = ( c + c − ) B ( i sin λt cos λt ) P tr (cid:16)Y T × (cid:17) . (B3)The other degrees of freedom can be systematically integratedout by repeatedly using the following trace identity tr( A A · · · )tr( B B · · · ) = tr [ ( A ⊗ B ) ( A ⊗ B ) · · · ] . (B4)Repeated use of the identity Eq. (B4) leads to enlarged trans-fer matrices as discussed below. Tracing out all degrees offreedom will lead to × transfer matrices the trace ofwhose product gives the Renyi entropy.Notice that at concave and convex corners, we get a differ-ent × transfer matrix. The × transfer matrices at thecorners will be given at the end of this appendix for complete-ness. We discuss below the procedure for obtaining the final × transfer matrices from the × ones along a smoothboundary. At the sharp corners, we can obtain the × transfer matrices from the corresponding convex and concavecorner × matrices essentially in the same way.Let us now explicitly write T × in Eq. (B3) along a smoothboundary. For a bond β i between two adjacent boundary pla-quettes as show in the top left corner of Fig. 6, we have T × = A α ′ i exp h iµtα ′ i β i α i α i +1 i × H A ( β ′ i + α i ) H − A − ( β ′ i + α i ) ! A β i A − β i A − β i A β i ! . (cid:0)(cid:0)(cid:1)(cid:1) (cid:0)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1) (cid:0)(cid:0)(cid:1)(cid:1) (cid:0)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1) (cid:0)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1) (cid:0)(cid:0)(cid:1)(cid:1) (cid:0)(cid:0)(cid:1)(cid:1) α ′ i β i α i a α i b β ′ i T L> : α ′ i α i +1 α i β i T : β i α i a β ′ i +1 β ′ i β ′ i β i +1 α i b T < : T R> : β i α i b α i a α ′ i FIG. 6. Top left: the transfer matrix T × along a smooth boundaryand the corresponding spin variables. Top right: the transfer matrix T L> × for entering a convex corner. Bottom left: the transfer matrix T R> × for exiting a convex corner. Bottom right: the transfer matrix T < × for a concave corner. The next step is to consider R , ≡ X β h α β | U + p U + s | Ω ih Ω | U − s U − p | α β i . We can now use the trace identity Eq. B4 to bring togetherterms with the same β and β ′ and trace out these degrees offreedom to get R , = ( | c + c − | ) B (sin λt cos λt ) P tr (cid:16)Y T × (cid:17) where the × transfer matrix above is given by T × = exp h i β µt (cid:16) α ′ i α i α i +11 − α ′ i α i α i +12 (cid:17)i t t t t with the following t i matrices t = (cid:18) H H − (cid:19) ⊗ H ∗ H ∗− ! ,t = A α i A − α i ! ⊗ A ∗ α i A ∗− α i ! ,t = X β ′ = ± (cid:18) A β ′ A − β ′ (cid:19) ⊗ A ∗ β ′ A ∗− β ′ ! ,t = X β = ± (cid:18) A β A − β A − β A β ′ (cid:19) ⊗ A ∗ β A ∗− β A ∗− β A ∗ β ! . Using the trace identity above to bring together the degreesof freedom via the tensor product of matrices is the main strat-egy for calculating the Renyi entropy. These calculations looktedious but are greatly simplified using Mathematica. Thenext step is to bring together the same degrees of freedomin R , R , using the trace identity which results in × matrices that only depend on α and α ′ degrees of freedom forboth replicas. Notice that so far we have integrated out all β and β ′ degrees of freedom (for both replicas).2Now since the α ′ degrees of freedom do not interact wecan do the sum directly. The α degrees of freedom have anIsing type interaction so the transfer matrix method can beapplied again promoting functions of α (which appear in × matrices) to × matrices (because there are replicas).This operation can be done very simply in Mathematica. Atthe end of the day the Renyi entropy will be given by the traceof the product of P (the number of boundary plaquettes) × transfer matrices.If the boundary of the subsystem is largely smooth, i.e. onlyhas a few sharp corners, the different transfer matrices at thecorners will give an order-one contribution to the Renyi en-tropy and the leading order area law term is S ( ρ A ) ≈ − P log Λ where Λ is the largest eigenvalue of the smooth boundary × transfer matrix constructed above.An exact calculation of the entanglement entropy is neededfor evaluating the topological entropy. Therefore correctionsfrom concave and convex corners must be taken into ac-count by using different transfer matrices at the corners in theboundary. There is a small complication due to the the factthat a convex (concave) corner plaquette has an additional α ( β ) spin. Let us review our method: we first consider two ma- trix elements each given by the trace of a product of × matrices, we then bring together the same degrees of freedomand explicitly sum over β s to get the × transfer matrices.Now having two β s to sum over does not change anything ex-cept we have to do a few more summations to get T × . Withour method of first integrating out the B spins and then the A spins, the additional α spin at a convex corner gives a jump inthe indices of the transfer matrices at the end of the day butthis can be easily accounted for by inserting a matrix C ij = 1 between the transfer matrices and calculating the trace of theirproduct as usual.Let us now write the × matrices for the corners. Firstconsider the transfer matrix T L> × for entering a convex corneras seen in the top right corner of Fig. 6. Each transfer matrixis written for a β bond and the α spins before ( α b ) and after α a that bond. As mentioned above the discontinuity in the α index ( α a for a transfer matrix is not equal to α b for the nexttransfer matrix) can be taken care of by inserting the matrix × ⊗ C × ( × because we have two replicas) at the endof the day between the × transfer matrices. For exitinga convex corner plaquette, with the notation of the bottom leftcorner of Fig. 6, we can write a × transfer matrix T R> × andsimilarly a transfer matrix T < × for a concave corner (bottomright corner of Fig. 6). In terms of the spin variables shownon the corresponding corner of Fig. 6, these matrices are givenby T L> × = A α ′ i exp h iµt α ′ i β i α i b α i a i × H A ( β ′ i + α i b ) H − A − ( β ′ i + α i b ) ! A β i A − β i A − β i A β i ! A α i a A − α i a ! T R> × = A α ′ i exp h iµt α ′ i β i α i b α i a i × H A α i b H − A − α i b ! A β i A − β i A − β i A β i ! T < × = exp (cid:2) iµt α i b α i a β i β i +1 (cid:3) × H A ( β i + β ′ i + α i b ) H − A − ( β i + β ′ i + α i b ) ! A β i +1 A − β i +1 A β i +1 A − β i +1 ! . Notice that for a concave corner, the transfer matrix is writ-ten between two neighboring boundary plaquettes which inthis case do not have a common edge because the plaquette inthe middle lies entirely in subsystem B and is therefore not aboundary plaquette. The procedure for going from the initial × transfer matrix to the final × is analogous to thesmooth boundary case discussed above and a detailed discus-sion would not be illuminating. We have explicitly verifiedusing the transfer matrix method described in this appendixthat the order-one correction to the area law is entirely fromthe convex and concave corners and the topological entangle-ment entropy remains identically zero. Appendix C: Renyi entropy at order h We stated in section IV that Eq. (22) gives zero topologicalentanglement at order h . Consider the partitions of Fig. 3 with the spins inside (but not on the boundary of) the blueregions belonging to subsystem B and all other spins to A .For partition I for example, the outer and inner boundaries ofthe subsystem B are squares of side ℓ and ℓ respectively.The whole system is on an ℓ × ℓ torus. In terms of these di-mensions, we can explicitly find the quantities n p A , n p B , n d A and n d B for the four different partitions and we observe thatthe value of each quantity for partitions II and III, which areobviously equal, is the average of those for partitions I andIV. This immediately implies that the topological entangle-ment entropy identically vanishes at second order in h . Forcompleteness we include the results of this simple geometry3 FIG. 7. The system for o ( h ) calculation. The dashed green bondsbelong to subsystem B and the solid black bonds to the subsystem A . exercise for partitions I and IV below n Ip A = ℓ − ℓ + ℓ n Ip B = ( ℓ − − ℓ − ℓ n Id A = 2( ℓ − ℓ − ℓ − ℓ + ℓ ) n Id B = 2 (cid:0) − ℓ + ℓ − ℓ − ℓ (cid:1) n IVp A = ℓ − ℓ ( ℓ − ℓ ) n IVp B = ( ℓ − ℓ − ℓ − n IVd A = 2 ℓ − ℓ + ℓ + ℓ (2 ℓ − n IVd B = 24 + 2 ℓ + 5 ℓ − ℓ (15 + 2 ℓ ) where we have assumed ℓ − ℓ > .Next, we discuss the fourth order term. The combinatoricsget more complicated at higher orders but we can calculatethe terms in the expansion by direct enumeration. This givessemi-analytical results with the time-dependence entering an-alytically through the multiple integrals depending on func-tions c , p and d and the geometry dependence computed fora system of given dimensions by computer enumeration. Wehave performed these calculation for several systems. Herewe present the details for an × system shown in Fig. 7.If tr( ρ A ) = 1 + f h + gh + · · · , we have S = − f h −
12 (2 g − f ) h + · · · Let us assume we have coefficients f i and g i for partition i = 1 · · · . We can then write the fourth order term in thetopological Renyi entropy as S topo (cid:12)(cid:12) o ( h ) = h (cid:2) g − g − g − (2 f − f − f ) (cid:3) . To find the fourth order term in S topo , in addition to the f coefficients that we found before, we need to calculatethe coefficients g , g and g for the different partitions.Looking at the structure of the Eq. (B1), we observe thatat fourth order we need to insert four V s into the four bins h· · · ih· · · ih· · · ih· · · i which give a total of terms of fivedifferent categories. There is one term of type (1 , , , (with one V in each h· · · i ) and
12 (2 , , , , , , , ,
12 (3 , , , and , , , terms. First consider the (1 , , , term. We define a matrix A such that A α,β ≡ h αβ | V | i and obtain g (1 , , , = tr h(cid:0) AA † (cid:1) i . Now since one V term acting on | i can only change the stateby flipping either one plaquette or one domino, the | αβ i stateswith a non-vanishing A α,β can be enumerated and the matrix A can be constructed. The elements of the matrix A containintegrals ( c ) , ( p ) and ( d ) depending on which term in V con-tributes to the matrix element. We then obtain for example thefollowing expression for partition 1 of the × system g , (1 , , , = 8 (cid:2) c )( c )( c )( c ) + 299( d )( d )( d )( d )+ 6336( d )( d )( p )( p ) + 34176( p )( p )( p )( p )+ 196608( c )( c )( d )( d ) + 2359296( c )( c )( p )( p ) (cid:3) . We can similarly calculate g (1 , , , for other partitions. Thenext two categories are (2 , , , and (2 , , , . The expres-sion for g (2 , , , and g (2 , , , are given below. g (2 , , , = 2 ( − h β | V V | ih | V † | αβ ih α | V | i− h α | V V | ih | V † | αβ ih β | V | i + h αβ | V V | ih | V † | β ih | V † | α i + c . c . (cid:1) and g (2 , , , = (cid:0) h β | V V | ih | V † V † | β i + h α | V V | ih | V † V † | α i + h | V V | ih | V V | i + c . c . ) . Computing these expressions in terms of the functions c , d and f with direct enumeration yields the following expressionfor the topological Renyi entropy. S topo (cid:12)(cid:12) o ( h ) = − (cid:8) d )( p ) [( p, d ) + ( d, p )]+( d ) − d, d ) −
16 [( p, d ) + ( d, p )] − p ) + 256( p ) ( p, p ) − p, p ) (cid:9) . Notice that the last category g (4 , , , is independent ofthe partitions and does not contribute to the topological en-tanglement entropy at fourth order in h . We can show thatthe expression in Eq. (C1) above identically vanishes using ( p, p ) = ( p ) / , ( d, d ) = ( d ) / and the following integralidentity ( p, d ) + ( d, p ) = ( d )( p ) . The fourth order calculation also confirms the notion thata growing ring-like (with the expansion order in h ) regionaround the system boundary contributes to the entanglemententropy. M. Greiner et al ., Nature (London) , 51 (2002); A. K. Tuch-man et al ., Phys. Rev. A , 051601 (2006); T. Kinoshita, T. Wenger, and D. S. Weiss, Nature (London) , 900 (2006); S. Hofferberth et al ., ibid . , 324 (2007). L. E. Sadler, J. M. Higbie, S. R. Leslie, M. Vengalattore and D.M. Stamper-Kurn, Nature , 312 (2006). I. Bloch, J. Dalibard and W. Zwerger, Rev. Mod. Phys. , 885(2008). M. Rigol, V. Dunjko and M. Olshanii, Nature , 854 (2008). N. Linden, S. Popescu, A. J. Short and A. Winter, Phys. Rev. E , 061103 (2009). K. Sengupta, S. Powell and S. Sachdev, Phys. Rev. A , 053616(2004). W. H. Zurek, U. Dorner, and P. Zoller, Phys. Rev. Lett. , 105701(2005). A. Polkovnikov, Phys. Rev. B , 161201(R) (2005). R. W. Cherng and L. S. Levitov, Phys. Rev. A , 043614 (2006). M. A. Cazalilla, Phys. Rev. Lett. , 156403 (2006). P. Calabrese and J. Cardy, Phys. Rev. Lett. , 136801 (2006). V. Gritsev, E. Demler, M. Lukin, A. Polkovnikov, Phys. Rev. Lett. , 200404 (2007). A. Polkovnikov and V. Gritsev, Nat. Phys. , 477 (2008); R.Barankov and A. Polkovnikov, Phys. Rev. Lett. , 076801(2008); C. De Grandi, R. A. Barankov, and A. Polkovnikov, Phys.Rev. Lett. , 230402 (2008). K. Sengupta, D. Sen and S. Mondal, Phys. Rev. Lett. , 077204(2008); D. Sen, K. Sengupta and S. Mondal, Phys. Rev. Lett. ,016806 (2008). S. Sotiriadis, P. Calabrese and J. Cardy, EPL (2009) 20002. P. Calabrese and J. Cardy, J. Stat. Mech. 2007, P06008. M. Fagotti and P. Calabrese, Phys. Rev. A , 010306(R) (2008). A. Y. Kitaev, Ann. Phys. (N.Y.) , 2 (2003). D. I. Tsomokos, A. Hamma, W. Zhang, S. Haas and R. Fazio,Phys. Rev. A , 060302 (2009). A. Kitaev and J. Preskill, Phys. Rev. Lett. , 110404 (2006). M. Levin and X.-G. Wen, Phys. Rev. Lett. , 110405 (2006). S. Trebst, P. Werner, M, Troyer, K Shtengel and C. Nayak, Phys.Rev Lett. , 070602 (2007). J. Vidal, S. Dusuel and K. Schmidt, Phys. Rev. B , 033109(2009); J. Vidal, R. Thomale, K. Schmidt and S. Dusuel, Phys.Rev. B , 081104 (2009). I. S. Tupitsyn, A. Kitaev, N. V. Prokofev, and P. C. E. Stamp, 82,085114 (2010). M. Rigol, V. Dunjko, V. Yurovsky and M. Olshanii, Phys. Rev.Lett. J. L. Cardy and I. Peschel, Nucl. Phys.
B300 , 377 (1988). E. Fradkin and J. E. Moore, Phys. Rev. Lett. , 050404 (2006). S. Papanikolaou, E. Luijten, and E. Fradkin, Phys. Rev. B ,134514 (2007). S.T. Flammia, A. Hamma, T.L. Hughes and X.-G. Wen, Phys. Rev.Lett. , 261601 (2009). C. Castelnovo and C. Chamon, Phys. Rev. B , 174416 (2007);C. Castelnovo and C. Chamon, Phys. Rev. B , 184442 (2007). Z. Nussinov and G. Ortiz, PNAS , 16944 (2009); Z. Nussinovand G. Ortiz, Ann. Phys. , 977 (2009). A. Hamma and D. A. Lidar, Phys. Rev. Lett.100