Exact results on the quench dynamics of the entanglement entropy in the toric code

# 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
July 26, 2019
###### Abstract

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 initial state with the toric code Hamiltonian. We find that the area law survives at all times. Adding disorder to the toric code couplings makes the entanglement entropy per unit boundary length saturate to disorder-independent values at long times and in the thermodynamic limit. There are order-one corrections to the area law from the corners in the subsystem boundary but the topological entropy remains zero at all times. We argue that breaking the integrability with a small magnetic field could change the area law to a volume scaling as expected of thermalized states but is not sufficient for forming topological entanglement due to the presence of an excess energy and a finite density of defects.

## I introduction

Cold atomic systems have provided an ideal playground for the experimental studies of the unitary evolution of thermally isolated quantum systems. Greiner02 (); Sadler06 (); Bloch08 () Apart from the experimental interest, studies of quantum evolution have shed light on some fundamental questions such as thermalization. Rigol08 (); Linden09 ()

A widely-studied scenario is the unitary evolution of a quantum state when the Hamiltonian is suddenly changed – a “quantum quench”, Sengupta04 (); Zurek05 (); Polkovnikov05 (); Cherng06 (); Cazalilla06 (); Calabrese06 () which has been the subject of many recent studies (see Refs. [Gristev07, ; Polkovnikov08, ; Sengupta08, ; Sotiriadis09, ] for a few examples). A challenging question, which has been addressed theoretically in one dimensional systems, is the evolution of the entanglement entropy following the quench. Calabrese07 (); Fagotti08 () More generally, such studies are concerned with the fate of an out-of-equilibrium quantum state with a Hamiltonian that changes with some time-dependent protocol.

The case where the coupling constants are changed across a quantum phase transition is of special interest. In recent experiments with spinor condensates, for example, the formation of ferromagnetic domains was observed in a quench to a Hamiltonian with a symmetry-breaking ferromagnetic ground state. Sadler06 () Even though the interactions after the quench favor a certain type of order, it is not at all obvious that the order should emerge out of equilibrium. Since symmetry-breaking order has a local order parameter, the non-equilibrium state can support the formation of ordered domains. But what about quenches across a topological phase transition? Would any of the topological characteristics of the equilibrium phase show up in the non-equilibrium wave-function?

Another interesting aspect of the non-equilibrium dynamics is the issue of thermalization and whether a large system can act as an effective heat bath for its subsystems. If thermalization occurs, due to the presence of some excess energy, the non-equilibrium state at long times should be in some sense close to a thermal state at a finite temperature, which implies volume scaling for the entanglement entropy.

The qualitative discussions above motivate a detailed analysis of a quantum quench across a topological phase transition into the topological phase. We study the specific case of the Kitaev toric code model Kitaev03 () with spin-polarized initial states. A related quantum quench problem with topologically ordered initial states was recently studied numerically in Ref. [Tsomokos09, ]. We focus on the time evolution of the entanglement entropy because it is intimately related to topological order through the topological entanglement entropy. Kitaev06 (); Levin06 () Entanglement entropy is also related to thermalization since, as mentioned above, volume scaling of the entropy as opposed to the area law is generically expected of thermal states.

We compute exactly the bipartite entanglement entropy as a function of time. This is an example where exact calculation can be done in 2+1D (thus far, analytical results for the time-evolution of entanglement had focused on 1+1D systems). We show that the entanglement entropy respects the area law at all times and if one adds disorder to the coupling constants, the recurrences or revivals are suppressed at long times and the entanglement entropy per unit “area” (boundary length) equals

 S/P=ln(4/e)

in the thermodynamic limit when the initial polarization is in the or direction.

We find the dependence of this saturation value on the initial state for an arbitrary direction of the initial magnetization. A new ingredient that comes in when considering these initial states is order-one corrections to the entanglement entropy from the convex and concave corners in the subsystem boundary. We find that the general form of the generated entanglement entropy is given by

 S(t)=f(t)P+f>(t)C>+f<(t)C<

where is the perimeter of the subsystem and () the number of convex (concave) corners in its boundary. The functions of time , and are independent of the subsystem geometry. From the form of the entanglement entropy above, we find that the topological entropy vanishes at all times.

We study the effects of breaking the integrability in a perturbative approach. We identify the mechanism by which an integrability-breaking magnetic field can change the area law to the volume law as expected of thermalized states. The physical picture is that the area law survives at low orders in an expansion in the integrability-breaking field. Going to higher orders however increases the thickness of a ring-like region around the subsystem boundary that contributes to the entanglement entropy. Eventually at high enough orders, the ring will collapse onto itself making volume scaling possible. In the perturbative treatment, we explicitly find that the topological entanglement entropy remains zero at low orders. We expect that the excess energy and the finite density of defects should prevent the formation of any topological order in a sudden quench even if one breaks the integrability.

The paper is organized as follows. In section II we discuss quantum quenches for spin-polarized initial states in the or direction. We calculate the entanglement entropy analytically and discuss the effects of disorder in the coupling constants. In section III, we study the role of the sharp corners in the subsystem boundary in the generation of entanglement by exactly calculating the second order Renyi entropy for an arbitrary direction of the initial magnetization. We then discuss the effects of breaking the integrability, which we treat perturbatively, in section IV. We conclude the paper in section V.

## Ii Integrable quench: initial magnetization in the x or z direction

Consider the toric code Hamiltonian in a field Trebst07 (); Vidal09 (); Tupitsyn10 ()

 H=−∑pλpBp−∑sμsAs−h∑i→n⋅→σi. (1)

The Hilbert space of the above Hamiltonian is given by spins living on the bonds of a square lattice on a torus and the sums are over all plaquettes, stars and bonds respectively. The plaquette and star operators are defined as follows

 Bp=∏i∈pσzi,As=∏i∈sσxi

where s are the Pauli matrices. For the usual Kitaev model the coupling constants and are independent of the plaquettes and the stars but here we keep the and indeces to allow for disorder in the couplings. Note that for the model in zero magnetic field (), the ground state is independent of the actual values of the coupling constants and as long as they are all positive.

One can study the dynamics of this model by changing the coefficients , and in time with some protocol. In this section, we consider a quantum quench where the system is initially in a fully magnetized state and then evolves with the integrable toric code Hamiltonian, namely, we are initially in the ground state for and for all stars and plaquettes. We then turn off the field at (which does not change the quantum state) and turn on the toric code part of the Hamiltonian as follows: and ().

Consider for now, the case where (). The ground state of the usual model (without disorder in the couplings) exhibits a topological phase transition: for the system is a paramagnet and for , it is topologically ordered. This topological phase is robust against weak disorder in the couplings. Let us call the initial state . In the case where the initial magnetization is in the direction in the basis. The initial state is a tensor product of single spin states and thus has no entanglement.

In the paramagnetic initial state we have for all spin operators. Before addressing the question of the entanglement entropy evolution, let us consider the simple problem of what happens to the expectation value. As a function of time, we have where is the evolution operator given by

 U=exp(i∑p∫t0λp(t′)dt′Bp+i∑s∫t0μs(t′)dt′As). (2)

Since all star and plaquettes commute with one another, we have dropped the time ordering operator. If we consider a sudden quench, we have and . Let us focus on sudden quenches for a more compact notation. At the end of the calculation we can replace and with the corresponding integrals over time-dependent couplings for a more general protocol.

The expectation value of , a spin between two nearest-neighbor plaquettes and is then given by

 ⟨σxpp′(t)⟩=⟨Ωx|e−itλpBpe−itλp′Bp′σxpp′eitλpBpeitλp′Bp′|Ωx⟩.

Since , we have

 e±itλpBp=cos(λpt)±isin(λpt)Bp.

We have . Therefore we get

 ⟨σxpp′(t)⟩=cos(2λpt)cos(2λp′t).

In the completely clean case (a theoretical abstraction), the average of over time is equal to but any amount of disorder will eventually lead to a state which unlike the initial paramagnet has a vanishing time-averaged expectation value for .

We next consider the main problem of this section, namely the evolution of the entanglement entropy. The density matrix of the system can be written as where is the time evolution operator and is a generic initial state. By inserting the resolution of identity twice, we can write,

 ρ(t)=∑α,σ⟨α|U|Ω⟩⟨Ω|U†|σ⟩|α⟩⟨σ|. (3)

We would like to calculate the entanglement entropy between two subsystems and . Any star or plaquette operator acts on four spins which may all belong to , all belong to or some belong to and some to . In this last case, the star or plaquette operator must lie at the boundary of the two subsystems. We can then write

 U=UAUBU∂ (4)

where and only act on spins in and respectively and only depends on star and plaquette operators at the boundary of the two subsystems. The three operators , and commute with one another. Each basis state can be written as the tensor product of states for spins in and , . We now integrate out the spins in to find the reduced density matrix for subsystem , and obtain,

 ρA(t)=∑αA,σA,βB⟨αAβB|U|Ω⟩⟨Ω|U†|σAβB⟩|αA⟩⟨σA|. (5)

Notice that since the trace of an operator is independent of the basis, we can replace by in the above expression. This amounts to calculating the trace in a rotated basis. We then substitute Eq. (4) into Eq. (5) above and obtain

 ρA(t)=∑αA,σA,βB⟨αAβB|UAU∂|Ω⟩⟨Ω|U†AU†∂|σAβB⟩|αA⟩⟨σA|.

The part of the Hamiltonian that only acts on spins has no effect on the time evolution of the reduced density matrix for subsystem . In order to calculate the entanglement entropy between the two subsystems we need to find . We give a replica index to , and and drop the subsystem index for brevity. States are in subsystem and in . Using we obtain,

 tr(ρnA)=∑{αi},{βi}n∏i=1⟨αiβi|UAU∂|Ω⟩⟨Ω|U†AU†∂|αi+1βi⟩ (6)

with . Once again because traces are basis-independent we can replace by in the above equation and we obtain an expression for that only depends on and not the whole . Let us define the following commuting operators

 U±p=e±i∑p∈∂pλpBpt,U±s=e±i∑s∈∂sμsAst (7)

where () is the set of all boundary plaquettes (stars) which act on both and spins. We then obtain

 tr(ρnA)=∑{αi},{βi}n∏i=1⟨αiβi|U+pU+s|Ω⟩⟨Ω|U−sU−p|αi+1βi⟩ (8)

with . We now focus on the simple case where initially all spins are aligned in the direction. In this case, is an eigenstate of and we can therefore substitute in Eq. (8) above. Also and are chosen in the basis so the action of on a basis state is to flip the four spins around the plaquette . For concreteness, we assume that the partition contains all spins which lie inside a simply connected rectangular region. We then have boundary plaquettes. Let us consider the action of on a state . Since , we have

 U±p=∏p∈∂p(cosλpt±iBpsinλpt). (9)

Notice that only if the state the state is obtained by acting with a certain number of plaquette flips on . Let us consider then Ising spin variables living at the center of plaquettes. A plaquette is flipped for and not flipped for . Notice that, on a torus, flipping a certain set of plaquettes and flipping all other plaquettes that do not belong to have the same effect on the spins, i.e. they describe the same element in the group of plaquette flips. However, since we have eliminated and from the formulation of the problem by a basis rotation, we can uniquely specify a state with nonvanishing by a set of flipped boundary plaquettes. Let the state be labeled by for . If the partitions are not too small, it is not possible for a set of boundary plaquette flips acting on to create a state with the same spins and different spins or vice versa. Therefore and we obtain using Eq. (8)

 tr(ρnA)=∑{θp},p∈∂p(⟨{θp}|U+p|Ωx⟩⟨Ωx|U−p|{θp}⟩)n. (10)

Using Eq. (9) we find that is a product of terms that are equal to for each unflipped plaquette with and for each flipped plaquette with . We can then write

 ⟨{θp}|U+p|Ωx⟩=∏p∈∂p(1+θp2cosλpt+i1−θp2sinλpt). (11)

Switching the order of the sum and product after plugging Eq. (11) into Eq. (10), we obtain

 tr(ρnA)=∏p∈∂p∑θp(1+θp2cos2λpt+1−θp2sin2λpt)n.

We have used the fact that for , .

Summing over for each boundary plaquette simply gives

 tr(ρnA)=∏p∈∂p[(cosλpt)2n+(sinλpt)2n]. (12)

By analytic continuation of , we write the entanglement entropy as

 S=−tr(ρAlog(ρA))=limn→1∂∂ntr(ρnA)

and we finally obtain,

 S=−∑p∈∂p[cos2λptln(cos2λpt)+sin2λptln(sin2λpt)]. (13)

For the clean Kitaev model, the entanglement entropy is an oscillatory function of time whose maximum is proportional to the number of boundary plaquettes or the perimeter of the subsystem. For this initial condition, the star term in the toric code Hamiltonian has no effect on the evolution of the entanglement entropy because the initial state is an eigenstate of that term. Similarly if initially all spins are in the direction, the plaquette term will have no effect. The completely clean system will never reach a steady state.

If there is some disorder in the coupling constants however, the entanglement entropy per unit boundary length will saturate to a steady state value in the thermodynamic limit. Suppose the number of boundary plaquettes is very large. As , due to the presence of some disorder in couplings , the phases are equally likely to take any value between and modulo . Therefore, taking the thermodynamic limit first we can write with

where is a random variable between and drawn from a uniform distribution. The average entanglement entropy is then given by

 ¯S=−P2π∫2π0dϕζ(ϕ)=Pln(4/e) (14)

We have which implies that for any realization of disorder

 |S({ϕi})P−¯SP|=O(1√P)

indicating that in the thermodynamic limit, the entanglement per unit area length saturates to a steady-state value . Notice that these results, namely the survival of area law at all times and the oscillatory behavior in the entanglement entropy in the absence of disorder are generic features of integrable quenches where the integrals of motion are operators with local support. The entanglement entropy per unit boundary length relaxes in the thermodynamic limit to a steady state value in the presence of disorder in the coupling constants.

Because the system can reach a steady state, it is natural to ask if the steady state of a subsystem has the properties of a thermalized state described by a Gibbs or a generalized Gibbs ensemble. Rigol07 () For such thermalized states however, we expect the entanglement entropy of a small subsystem to scale as its volume in the limit of long times. The survival of the area law at all times which follows trivially from having local integrals of motion indicates the absence of thermalization even to a generalized Gibbs ensemble.

## Iii Subsystem corners and the dynamical entanglement entropy

Let us now consider a more generic case where the initial density matrix is not diagonal in either or basis and both terms in the Kitaev Hamiltonian affect the entanglement dynamics. We find that even though in the thermodynamic limit the area law still applies, the physics of the out of equilibrium entanglement is much richer. We find that here there are order-one corrections to the entanglement entropy from the concave and convex corners of the subsystem.

Subleading corrections to the entanglement entropy arising form the sharp corners in the boundary were found before for static entanglement entropy in the logarithmic corrections to the area law in a class of conformal quantum critical points. Cardy88 (); Fradkin06 (); Papanikolaou07 () Here we find that corners can also make subleading contributions to the dynamical entanglement. We find that when subsystem consists of spins inside a region with convex and concave corners as seen in Fig. 1, the entanglement entropy after a quantum quench behaves as

 S=fP+f>C>+f

where is the perimeter of the subsystem and , and are functions of time that do not depend on the subsystem geometry.

For a subsystem with a mostly smooth boundary and only a few sharp corners, the leading term is the but for a subsystem with a corrugated boundary, the entanglement entropy can be different in the thermodynamic limit.

Notice that the number of concave and convex corners are not independent. Suppose subsystem consists of spins inside of an area with connected regions, each with a certain number of handles such that the total number handles is . As an example in Fig. 1, we have and . We argue that

 C>=C<+4(m−g).

This is because for a rectangular simply connected region, we have and and transforming the shape by adding or subtracting plaquettes on the square lattice does not change . Therefore each simply connected region makes a contribution of to and by a similar argument, each hole contributes .

### iii.1 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 direction. Let us work in the the basis with , and . The action of a star or plaquette operator in this basis is then to flip the corresponding four spins modulo a phase factor of from . We will see later that this sign has no effect on the entanglement entropy.

We assume for simplicity that there is a minimum distance of times the lattice spacing between all corners. Consider a convex corner as in the right hand side of Fig. 2. The spins live in the dashed blue (color online) bonds and the spins on the solid black bonds. Flipping the corner plaquette represented by a red square has the same effect on spins as flipping the two boundary stars represented by red crosses. These operations of course have a different effect on spins. Because we assumed a minimum distance between the corners, these corner star and plaquette operators uniquely correspond to one corner. Similarly as in the left hand side of Fig. 2, flipping the corner star at a concave corner has the same effect on the spins as flipping the two corresponding plaquettes. We see from this discussion that the effects of stars and plaquette flips mix at the corners, which leads to a different entanglement 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 entropy for a subsystem consisting of the bonds inside an rectangular region. A boundary star or plaquette that is not located at a corner is refereed to as a regular star or plaquette. As seen in Fig. 2, at each convex corner, one plaquette and two stars are special. For our subsystem, we have four convex corners, regular boundary plaquettes and regular boundary stars and the entanglement entropy for a clean system after a sudden quench is given by

 S=2(M+N−4)sp+2(M+N−6)ss+4sc (16)

where the contributions from regular plaquettes (), regular stars () and convex corners () are given by

 sp = (sin2λtlnsin2λt+cos2λtlncos2λt) ss = (sin2μtlnsin2μt+cos2μtlncos2μt) sc = [r1lnr1+2r2lnr2+r3lnr3]

with

 r1 = sin2λtsin4μt+cos2λtcos4μt r2 = sin2μtcos2μt r3 = sin2λtcos4μt+cos2λtsin4μt.

The above solution conforms to the generic form Eq. (15) for , , , and .

In appendix A, we give a derivation of the above result in the more general setting of an arbitrary subsystem geometry. If instead of a sudden quench we turn on the coupling constants with some protocol and , we just need to replace and by their integrals over time.

Using the partitions of Ref. Levin06, , shown in Fig. 3, the topological entanglement entropy is given by

 Stopo=SI−SII−SIII+SIV. (17)

From the form of Eq. (15), we can then see that although there are order-one corrections to the area law, the topological entanglement entropy remains exactly zero at all times independently of how the toric code Hamiltonian is turned on.

It seems plausible that for any initial state, as long as the integrability is not broken (with a field for example), evolution with the Kitaev Hamiltonian does not change the topological entanglement entropy for any time dependence of the coupling constants. This is because this integrable evolution does not allow propagation of information.

### iii.2 Arbitrary initial magnetization

The von Neumann entropy is a special case (limit of ) of the order Renyi entropy defined as

 Sn(ρA)=11−ntr(ρnA).

Hereafter we refer to the second order () Renyi entropy as the Renyi entropy. The Renyi  entropy is simpler to calculate than the von Neumann entropy and is an equally valid measure of entanglement. It has been shown that the topological entanglement entropy calculated by any order Renyi entropy contains the same information in the abelean topological phases Flammia09 () such as the toric code. For an integrable quench from a paramagnetic state with magnetization in an arbitrary direction, we exactly calculate the Renyi entropy . We consider the clean case with position-independent coupling constants. With disorder in the couplings, the entropy saturates to the time-average of the clean case similarly to the simpler initial conditions considered before.

In the basis, an arbitrary paramagnetic state can be written as

 |Ω⟩=(c+|+⟩+c−|−⟩)⊗N.

The calculation is rather involved and is done by a transfer matrix method explained in appendix B. We find the same general features for an arbitrary initial magnetization as for the initial magnetization in the direction. Let us summarize these features below.

The Renyi entropy can be written as a sum of contributions from the perimeter of the subsystem as well as the convex and concave corners. These contributions are oscillatory functions of time in the absence of disorder in the couplings which at long times, saturate to the values given by the time-average of the clean case if we have some disorder in the couplings. The saturation values are independent of the final coupling constants and the protocol used for turning them on. They do however depend on the initial state.

For a large subsystem with only a few sharp corners, the leading term for the Renyi entropy is the first term in Eq. (15), namely . The saturation value of the Renyi entropy per unit boundary length is plotted in Fig. 4 as a function of the and angles which determine the initial state through and . Notice that this saturation value is a unique function of and .

## Iv Breaking the integrability: entanglement entropy in the interaction picture

To study the effects of breaking the integrability, we focus on one simple case. Namely, we start with a state that has all spins aligned in the direction. As opposed to the integrable case where the magnetic field is turned off before turning on the integrable toric code Hamiltonian, the the toric code Hamiltonian is turned on in the presence of a small magnetic field . We cannot perform exact analytical calculations in this case. We can however discuss the effects of the magnetic field by calculating a formal expansion for the (Renyi) entanglement and topological entropies in powers of . The expansion can only be considered as a perturbative treatment if the toric code coupling constants are quickly pushed above 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 eventually lead to volume scaling in the entanglement entropy. The physical picture is that in the expansion in powers of , the contribution to the entanglement entropy comes from ring-shaped areas around the subsystem boundary whose thickness grows at higher and higher orders in . Eventually, at a high enough order, these correlated regions meet and the terms in the expansion no longer respect the area law. Even for a small , these high order terms become important at long times making volume scaling of the entanglement entropy possible as expected of thermalized states. Secondly, we find that in this expansion in powers of , the topological entanglement entropy remains zero at low orders until we reach an order proportional to the system size, suggesting that a small integrability-breaking term is not enough for generating topological entanglement.

A useful notion that we introduce in this section is the entanglement entropy in the interaction picture. We are used to the Shcrödinger, Heisenberg and interaction pictures for calculating the expectation values of observables, which are of course independent of the picture. If one defines the entanglement entropy as a property of the wave function alone, then since the wave function does not change in the Heisenberg picture, calculating the entanglement entropy in that picture leads to the nonsensical result that time evolution does not change the entanglement of a quantum state. An interesting alternative is to try and define the entanglement entropy entirely in terms of the correlation functions of appropriate operators. We do not pursue this path here. We argue instead that as explained below, the entanglement entropy calculated only from the wave function in the interaction picture can be a useful concept.

Let us first consider a non-interacting Hamiltonian such as a paramagnet. Consider a partitioning of the system into subsystems and . The evolution operator for can be written as

 U0(t)=U0A(t)U0B(t)

where () only acts on the degrees of freedom in (). Notice that acting with on any wave function does not change its entanglement properties. Now consider a Hamiltonian

 H=H0+V

with the corresponding evolution operator .

The wave functions at time in the Schrödinger and interaction pictures are respectively given by

 |ψS(t)⟩=U(t)|ψ(0)⟩,|ψI(t)⟩=U†0(t)U(t)|ψ(0)⟩.

Now since acting with or on any wave function does not change its entanglement, the interaction picture wave function has the same entanglement entropy as the Schrödinger picture. Note however the calculation may be much simpler in the interaction picture.

In the discussion above, we made use of the interaction picture entanglement entropy because we had a Hamiltonian whose evolution operator did not change the entanglement entropy. At the next level of complexity, where using the interaction picture wave function may be useful for calculating the entanglement entropy, the evolution operator of the unperturbed Hamiltonian () may have another known property. For example it may only create entanglement entropy scaling with the subsystem area or it could be incapable of changing the topological entanglement. This is the case for the integrable toric code Hamiltonian. Now if we are interested in finding whether a perturbation can generate entanglement entropy with volume scaling or whether it can generate topological entanglement entropy, we can perform the calculation for the wave function in the interaction picture.

Let us go back to our quench problem and write the wave function in the interaction picture with

 H0(t)=−∑pλ(t)Bp−∑sμ(t)As

and

 V=h∑iσxi.

Starting with the initial state with all the spins in the positive direction, the wave-function at time in the interaction picture is given by

 |ψI(t)⟩=ei∫t0dt′H0(t′)|ψS(t)⟩=^UV(t)|0⟩ (18)

where is given by

 ^UV(t) = 1−ih∫t0dt1^V(t1) + (−ih)2∫t0dt1∫t10dt2^V(t1)^V(t2)+⋯

for

 ^V(t)=he−i∫t0dt′λ(t′)∑pBp(∑iσxi)ei∫t0dt′λ(t′)∑pBp.

For the simplicity of notation, we focus on the sudden quench and write as , bearing in mind however that for a more general protocol, we can substitute the integral at the end of the calculation.

We can explicitly calculate and the result is

 ^V(t) = cos4λt+12C+sin4λt2P+cos4λt−12D (19) ≡ c(t)C+p(t)P+d(t)D

with the operators , and defined as follows

 C = ∑iσxi P = D = (20)

The sums are over all bonds, plaquettes and dominos and the graphical notation represents a product of Pauli matrices. For example if the bonds in a plaquette are numbered clockwise from to with the top bond being number we have

The closed form in Eq. (19) is simply derived from the repeated substitution of the following commutation relations

 [∑pBp,C]=2iP,[∑pBp,P]=−4i(C+D) [∑pBp,C+D]=4iP

in the Baker-Hausdorff expansion of

 ^V(t) = h(∑iσxi+(−iλt)[∑pBp,∑iσxi] +(−iλt)22![∑pBp,[∑pBp,∑iσxi]]+⋯)

and resumming the series.

The Renyi entropy for the interaction picture wave function is given by

 S2(t) = −log⟨α1β1|^UV|0⟩⟨0|^U†V|α2β1⟩× (21) ⟨α2β2|^UV|0⟩⟨0|^U†V|α1β2⟩

where summation over repeated indices is implied and as before, the and spins belong to the and subsystems respectively. First let us consider the first order term in the (the expression for the Renyi entropy without the log). We use the shorthand notation

 VV⋯Vn times=∫t0dt1⋯∫tn−10dtn^V(t1)⋯^V(tn),

and similarly for the Hermitian conjugate. The first order term vanishes as seen below

 tr(ρ2A)=1−2ih⟨0|V|0⟩+2ih⟨0|V†|0⟩+⋯=1+O(h2).

The second order contribution to is given by

 tr(ρ2A)∣∣o(h2) = h2[2|⟨0β|V|0⟩|2+2|⟨α0|V|0⟩|2−⟨0|V|0⟩2 − ⟨0|V†|0⟩2−2⟨0|VV|0⟩−2⟨0|V†V†|0⟩].

Let us assume that the whole system is on a torus and the total number of bonds (spins) is . The number of plaquettes and the number of dominos in the system is . Let us introduce another shorthand notation for integrals involving the functions , and ,

 (f1,f2,⋯,fn)=∫t0dt1⋯∫tn−10dtnf1(t1)⋯fn(tn)

so we have for example and . We denote the number of plaquettes (dominos) with all bonds in subsystem by () and similarly by () for subsystem . We then have

 |⟨0β|V|0⟩|2=Nb2(c)2+npB×16(p)2+ndB(d)2 |⟨α0|V|0⟩|2=Nb2(c)2+npA×16(p)2+ndA(d)2 ⟨0|V|0⟩=Nb(c) ⟨0|VV|0⟩=Nb2(c,c)+Np×16(p,p)+Nd(d,d).

The factors of come from the fact that each term in is a sum of four operators that flip the spins around a plaquette. Noting that because is obviously symmetric under the exchange , we have and similarly for and . We then obtain

 tr(ρ2A)∣∣ = 1−h2(Np−npA−npB)(p)2 − h2(Nd−ndA−ndB)(d)2,

which using gives

 S2 = h2(Np−npA−npB)(p)2 + h2(Nd−ndA−ndB)(d)2+O(h4).

This expression demonstrates the physical picture for building volume entanglement which we alluded to before. We observe from the above expression that at this order the contribution to the entropy is from a ring-shaped area of thickness two lattice spacings around the subsystem boundary. At higher orders, the thickness of this correlated region, which contributes to the entanglement entropy, grows until the ring closes in on itself, making volume scaling possible.

Notice that in contrast to exact area law which is possible in a finite system, since the Renyi entropy is symmetric under exchanging the subsystems, volume scaling can only hold exactly for a finite subsystem in an infinite system and as a good approximation in the limit when a subsystem is much smaller than the whole system.

It is a simple exercise in geometry to count and of the two subsystems for partitions shown in Fig. 3 and verify directly that no topological entanglement is generated at order . The details of this calculation will be given in appendix C. We will also show in that Appendix that at the next order with a non-vanishing contribution to the Renyi entropy (), the topological term will still vanish due to an exact cancellation.

Similarly to the volume scaling, topological entanglement is non-perturbative and does not show up in an expansion in powers of in the thermodynamic limit. For a finite system, we expect non-vanishing contributions to the topological entropy to appear at orders that scale with the system size. Let us comment here that on physical grounds, in the case of a sudden quench, we do not expect any topological entanglement to form even at infinitely long times. A trend toward forming volume entanglement observed at low orders suggests that in the limit of large times, the generic thermalization behavior is expected. Although we cannot give a mathematical proof, the fragility of topological order to thermal fluctuations Castelnovo07 (); Nussinov09 () and the extensive excess energy in the system suggests that the non-equilibrium state, after a sudden quantum quench, will likely bear no signature of the equilibrium topological order of the phase.

We emphasize that the expansion used here cannot be treated as a perturbative treatment in the adiabatic limit where the toric code coupling constants are turned on slowly because in that case for short times. We know that in the case of adiabatic evolution, we approximately remain in equilibrium, i.e. in the ground state manifold, and topological order emerges after a time of order system size. Hamma08 ()

## V conclusions

We studied the quench dynamics of the entanglement entropy in the Kitaev toric code model. In the case of the integrable quench we found the entropy exactly as a function of time. The leading order contribution was found to scale as the subsystem perimeter with an oscillatory behavior that gets suppressed in the presence of some disorder in the coupling constants leading to the saturation of the entanglement entropy to disorder-independent values respecting the area law.

We obtained results for the generated entanglement entropy as a function of an initial paramagnetic state with arbitrary magnetization. We studied the effect of the sharp corners in the subsystem boundary and found that in analogy with the static entanglement discussed in the literature, they make subleading contributions to the dynamically generated entanglement entropy as well. We showed that the topological entanglement entropy remains zero in this quantum quench. The results of the integrable quench as well as the analysis of the integrability-breaking perturbations strongly suggest that the non-equilibrium state after a quench into the topological phase does not develop any of the topological structure present in the equilibrium phase.

We used the concept of the entanglement entropy in the interaction picture which simplifies perturbative approaches to the calculation of the entanglement entropy. This allows us to tease apart in this particular problem the effects of the time evolution with an unperturbed Hamiltonian that cannot generate either entanglement with volume scaling or topological entanglement and focus only the part of the evolution that can possibly lead to volume scaling or topological entanglement. By studying the effects of an integrability-breaking perturbation, we found a mechanism where volume scaling of the generated entanglement entropy can emerge at high orders in perturbation theory as expected of thermalized states. The notion of thermalization after a quench naturally ties in with topological order through the fragility of topological order to thermal fluctuations.

###### Acknowledgements.
We are grateful to C. Castelnovo, A. Hamma, M. Haque, A. Polkovnikov and M. Rigol for helpful discussions. This work was supported by in part by the DOE Grant DE-FG02-06ER46316.

## Appendix A Derivation of entropy for y direction initial magnetization

We explicitly calculate in this appendix the contributions to the entanglement entropy from regular boundary stars and plaquettes as well as the convex and concave corners when the initial magnetization is in the direction.

First let us find the number of regular boundary stars () and plaquettes () in terms of the subsystem perimeter and the number of convex and concave corners. We define the perimeter as the total number of boundary plaquettes or boundary stars (they are equal). We then have

 Pr=P−C>−2C<,Sr=P−2C>−C<.

Once again we use Eq. (8) which generally holds for any initial state . The only intermediate states (in the basis) which contribute to the are the ones connected to by some star or plaquette flips. We label these states by variable and living on the boundary plaquettes and stars respectively. A negative () indicates that the corresponding plaquettes (star) is flipped with respect to .

 |αβ⟩=|{θpr,ψsr,θp>,ψjp>,ψs<,θjs<}⟩ (23)

where , , and label regular plaquettes, regular stars, convex corner plaquettes and concave corner stars respectively. For , () are variables living on the two stars (plaquettes) corresponding to the convex (concave) corners. Notice unlike regular star variables the two stars at a convex corner are labeled by the corner plaquette and similarly the two plaquettes at a concave corner are labeled by a corner star since a convex (concave) corner uniquely corresponds to a corner plaquette (star).

Now suppose two such states have the same spins. All regular and concave corner plaquette and star flips must be the same for these two states. As seen in Fig. 2 however, at a convex corner the or variables for these two states could be different if both the convex corner plaquette and the corresponding two corner stars simultaneously have opposite flips. Then for two states and (which can be labeled as Eq. ( 23) because they contribute ) all their and variables must be the same except for the ones corresponding to convex corners, namely and which follow the relation

 {θp>,ψjp>}′={γp>θp>,γp>ψjp>}

with auxiliary variable . Similarly we can introduce variables for concave corner and for two states and we must have

 {ψs<,θjs<}′={ρs<ψs<,ρs<θjs<}.

Introducing a replica index on and , we can write

 |α1β1⟩ = |{θpr,ψsr,θp>,ψjp>,ψso,θjso}⟩ |α2β1⟩ = |{θpr,ψsr,γ1p>θp>,γ1p>ψjp>,ψs<,θjs<}⟩ |α2β2⟩ = |{θpr,ψsr,γ1p>θp>,γ1p>ψjp>,ρ1s<ψs<,ρ1s<θjs<}⟩ |α3β2⟩ = |{θpr,ψsr,γ2p>θp>,γ2p>ψjp>,ρ1s<ψs<,ρ1s<θjs<}⟩ ⋯

Let us for the moment focus on the Renyi entropy for simplicity. The calculation for will be very similar. First we can see from the form of the the states above that the sign factors from acting on states in the basis have no effect. A plaquette flip can only give a different sign for the matrix elements involving and at a corner. But then the plaquette flip at that corner would also give a different sign for the matrix elements involving and and the overall sign for the product of the four matrix elements appearing in will be positive.

The calculation is then similar to the simpler cases with the initial magnetization in the and directions except for each convex or concave corner there is an additional auxiliary variable to sum over. Again, each matrix element can be written as a product similar to Eq. (11). By simply bringing together terms we can write in the following form:

 tr(ρ2A)=∑∏prfpr∏srfsr∏p>fp>∏s

where the summation is over all , , and variables for both replicas. We can take the sum Eq. 24 inside the products by switching their order and perform the sum for individual terms. This gives the factorized form below

 tr(ρ2A)=∏prFpr∏srFsr∏p>Fp>∏s

where we will compute the s below. Assuming no disorder and taking the log of both sides, we immediately obtain our main result for the entropy

 S = (P−C>−2C<)lnFpr + (P−2C>−C<)lnFsr+C>lnFp>+C

which conforms to the general expression Eq. (15).

For r