SciPost Submission Page

Quenching the Kitaev honeycomb model

by Louk Rademaker

Submission summary

As Contributors: Louk Rademaker
Arxiv Link:
Date submitted: 2019-07-08
Submitted by: Rademaker, Louk
Submitted to: SciPost Physics
Domain(s): Theoretical
Subject area: Condensed Matter Physics - Theory


I studied the non-equilibrium response of an initial N\'{e}el state under time evolution with the Kitaev honeycomb model. With isotropic interactions ($J_x = J_y = J_z$) the system quickly loses its antiferromagnetic order and crosses over into a steady state valence bond solid, which can be inferred from the long-range dimer correlations. There is no signature of a dynamical phase transition. Upon including anisotropy ($J_x = J_y \neq J_z$), an exponentially long prethermal regime appears with magnetization oscillations. The prethermal regime can be interpreted using the toric code.

Current status:
Editor-in-charge assigned

Author comments upon resubmission

\section{Referee 1}

I will answer point by point the question from the referee. However, since some points have clear overlap (for example (1), (A) and 1.), I’ve grouped them together to make the discussion as streamlined as possible.

\subsection{Questions: (1) / (A) / 1.}

It is easy to see that $\langle \omega_p \omega_{p’} \rangle = 0$ for $p,p’$ two different plaquettes, because the plaquette operator $\omega_p$ flips the spins on the 1, 2, 4 and 5 sites of the plaquette $p$. To get a non-zero expectation value you need to ‘repair’ such a flip and that cannot be done with a finite number of plaquette operators.
Now $\langle \omega_p \omega_{p’} \rangle = 0$ implies that the ‘equally weighted superposition of a vortex-free and a vortex-lattice state’ is not the correct way to describe the Neel state.
There is one minor issue, that is actually addressed at the bottom of Page 3: for a finite size system with periodic boundary conditions not all flux configurations contribute. On page 3 I emphasize that the Neel state is an eigenstate of $\sigma^z$ operators. Each such $\sigma^z$ operator changes the gauge/matter configuration, and by systematically applying $\sigma^z$ operators one can probe all the relevant gauge/matter configurations that compose the Neel state. Notice that this leads, in the end, to a superposition of $N_c = 2^{3L_x L_y - L_y}$ different matter-gauge field configurations, and thus to almost all flux configurations. I have extended the discussion on the Neel flux configurations in the final manuscript.

\subsection{Questions: (2) / (3) / (B) / 2. / 3. }

There are two main questions posed here.

1. The first is regarding the meaning of the phrase ‘prethermalization’. My wording here is completely based on References 15-20, and in particular Ref. 18: "A Rigorous Theory of Many-Body Prethermalization" by Abanin et al. The key ‘strangeness’ of prethermalization is that timescales scale exponentially in the model anisotropy.

Naively, the speed at which some excitation relaxes is some power of the relevant Hamiltonian parameter. The simplest example is just a single-parameter Hamiltonian like the Heisenberg model: double $J$, and the speed of relaxation is halved. Another simple example is the dynamics close to a quantum phase transition, which scales with some (anomalous) power of the separation from the critical point. To have any kind of observable or measurable quantity change exponentially when I change a Hamiltonian parameter linearly is quite unique.

So when does this happen? Well, if you are close to an integrable system, but not just to any integrable system. This is mentioned explicitly on Page 5: the integrable Hamiltonian $H_0$ itself should be a projector model, that is the sum of locally commuting terms. Deviations from such form, so $H = H_0 + \lambda H’$ lead to an exponential time-scale $\mathcal{O}(\lambda)$ in which the integrable part $H_0$ is approximately conserved. I think it’s pretty clear that the Kitaev model in the anisotropic limit satisfies this set-up. It’s a bit strange that the referee asks for “some references” even though I explicitly mention on page 5: "We thus find the emergence, for the anisotropic model, of a distinct prethermalized regime. This can be understood using the framework of Ref. [18].”
(Below, at point C, I will also discuss the commonplace confusion about ‘integrability’ and naive preclusion of thermalization itself.)

2. The second question is whether these results are not very trivial since we project onto, quote, “a nearly flat energy band”. The Kitaev model is exactly solvable, yes, but it is not a model of free fermions and to speak of “flat energy bands” is not correct. I guess the referee means that the distribution of energy eigenvalues of the initial Neel state is quite narrow, and indeed, this width scales as $Var(E) \sim J_{xy}^2$. However, that still begs the question why the expected slower dynamics are not scaling as some power of $J_{xy}$? To reiterate, such an exponential dependence on the model parameters is almost a textbook definition of prethermalization.

Small comment: I agree that choosing an initial ferromagnetic state will lead to the same behavior. I can add this as a comment to the paper.

\subsection{Question: (C) }

On the one hand, the Kitaev model is not a free-fermion model, and it is the interference between different gauge sectors that makes up the time dependence of physical observables. This is made explicit several times throughout the paper, for example in Eqn. (4) where the magnetization (aka a physical observable) is the result of looking at the overlap of time-evolved states in different gauge sectors.

On the other hand, it is a common but very wrong assumption that free systems do not thermalize. Every reasonable translationally invariant system eventually thermalizes to its appropriate generalized Gibbs ensemble or to the relevant diagonal ensemble. Now in the case of the Kitaev model, this means that the final state can be written as a GGE with additional generalized inverse temperatures for the localized plaquette operators - since they are the relevant integrals of motion. The thermalization towards the Diagonal ensemble, see appendix B, is shown in the paper.

\subsection{Question (D) }
They don’t ‘disappear’, they are just not included. In Fig. 1 we focus on the long-time behavior, so how the system goes initially from $m=1$ to $m=0.6$ or so is not really the point of the paper. Notice that the data for $J_{xy} = 0.2 J_z$, that is for a higher anisotropy, is included with a finite size scaling in Fig. 5.

On a related note: I noticed an error in the paper, I forgot to include the fact that I used the algorithms of M. Wimmer (\verb+ to compute the pfaffians in Mathematica. I've included this information in the appendix.

\subsection{Question (E) }
I agree, basically the whole paper is a gem. I’m happy to rephrase “Kitaev’s genius was his realization” into “Kitaev's key insight was” but in general I’d like to avoid dry writing.

\subsection{Question (F) }
I know, I mentioned the brick-wall representation explicitly (Ref. 14). However, for the current purposes I found the method of Kitaev’s majoranas more amenable.

\subsection{Question (G)}
I respectfully disagree with the Referee here. Yao and Qi showed that if a state can be written as “gauge-fermion product states” $|u \rangle \otimes | \phi(u) \rangle$, where $|u \rangle$ is the gauge configuration and $|\phi(u) \rangle$ a matter fermion configuration, the resulting entanglement entropy can be split between the entropy generated by the gauge configurations and entropy generated from the fermions. Obviously, this is then true for all eigenstates. However, all bets are off once you include superpositions of such ‘gauge-fermion product states’: on the one extreme, the state suggested by the referee (equally weighted superposition of a vortex-free and a vortex-lattice state) clearly still has topological entanglement entropy, whereas the Neel state obviously has no entanglement whatsoever.

Also, the suggested paper from Bray-Ali does not discuss fermion-gauge models so I don’t understand its relevance for the gauge degrees of freedom.

In conclusion, I addressed the two major critiques:
\item As for prethermalization, I follow the definition of Ref 18. It is clear that the anisotropic Kitaev model satisfies the conditions layed out in that paper. This can explain the otherwise unexpected exponential dependence of the timescales on the anistropy parameter.
\item As for the initial Neel state, I showed that also expectation values of products of plaquette operators are zero in the Neel state, thus effectively providing almost all the possible flux configurations that one can have. However, as can be seen on page 3, for systems with periodic boundary conditions not all gauge configurations contribute and the new version will reflect this insight.

\section{Referee 2}

I will address point by point the issues raised in the section ‘Requested changes’:

\subsection{Question 1.}

The Referee is absolutely right in that in the original manuscript there was no clear discussion of the possible topological aspects of the quench I studied. I thought about this for a long time, and I have decided that it is better to present the work without much focus on the topological aspects. After all, there is no clear measure of topology that I have studied, neither in terms of anyons nor in terms of entanglement. I feel that such a study, while interesting, is more something for a follow-up work.

Therefore, I have rewritten large parts of the text (including the introduction and abstract) to focus on the dynamical transition from a magnetic state to a non-magnetic state. The discussion of a possible dynamical phase transition or crossover, the prethermal regime and the final steady state are better captured in terms of generic quench dynamics. A discussion of the topological aspects is reserved for the final Discussion section.

\subsection{Question 2.}
I removed the last sentence and the comment about the method from the abstract, and put more emphasis on the explicit results.

\subsection{Question 3 / Question 9}
Following the overall change of emphasis as explained in the answer to Q1, I have added a longer discussion on topological aspects to the Discussion section. I hope I thereby also answered Question 9 satisfactory.

\subsection{Question 4.}
The Referee asks here about a comparison between the entanglement structure expected in 'simple' quenches such as the XY model and the quench in the Kitaev model. As answered by Question 1, I have decided to drop the emphasis on entanglement as I could not say anything decisive about it. I feel that therefore there is no also need to discuss entanglement in XY models. I did, however, add some notes on timescales in the transverse field Ising model in response to Question 6.

\subsection{Question 5.}
To clarify the notion of a dynamical phase transition, I have added a separate subsection 3.1 focusing completely on the question of whether there is a phase transition when quenching the Kitaev model.

\subsection{Question 6.}
I have added a discussion, in section 3.3, showing that for quenches in the transverse field Ising model the typical timescale diverges as a power-law, as is expected in most systems. This emphasises the special nature of the exponential long prethermal regime.

\subsection{Question 7.}
It has been shown in Ref.~[26] that a valence-bond solid (VBS) has long-range dimer-dimer correlations. Therefore I have computed this quantity to test whether VBS order existed in the steady state. I have changed the text in the new section 3.3 to reflect the special role of the dimer correlations. Note that measuring a correlation function as a test of long-range order is not restricted to zero temperature.

\subsection{Question 8.}
The dynamic two-time spin correlation function as a function of frequency is similar to AC conductivity, which for charged systems can be expressed in terms of the dynamic density-density correlation function. I understand, however, how references to a "Drude peak" can be confusing, so I have removed such mentions and made it clearer what the $\omega = 0$ peak of the correlation function means.

List of changes

This version is a major revision, see the answers to questions from Referees in the "Author Comments" section.

Reports on this Submission

Anonymous Report 2 on 2019-8-8 Invited Report


1- The problem of quenching and time evolving under a Hamiltonian known to host topological and complex spin phases and obtaining concrete results on such non-equilibrium dynamics is an excellent one to look at.
2- Employing the Kitaev honeycomb model and its exactly solvable nature makes for a good approach. The methods used are very nice; expressing the product state in terms of variables used to diagonalize the Kitaev model to solve the quench problem is an innovative addition to the existing literature.

3- The care taken to address issues, such as doing a more careful study and explaining then notion of the dynamic phase transition, is apparent in the resubmission. A clearer message has been given on quench results between the Neel and valence bond solid phases.


1. The main draws of focus on topological and entanglement properties in the quenching procedure have now been removed. But given that they were nebulously stated anyway and now the author does a much more dedicated study of aspects that can be rigorously derived, this weakness is compensated for.


It is a pity that in this resubmission the aspects that I was most excited by – effects of quench on topological and entanglement properties – have now been removed. However, the author has put significant efforts into improving the aspects where he had concrete statements, namely, the behavior of dynamic quantities such as non-equilibrium free energy density and magnetization. The results are non-trivial with regards to non-equilibrium features, though sans dynamic phase transition, and pre-thermalization properties.

There are still some outstanding issues that are to be addressed:

1. The other referee has brought up excellent points on the statements made on the toric code and pre-thermalization. This would require an involved discussion.

2. There has been considerable discussion in the reports and response with regards to free fermions, etc. I believe that if the initial state chosen were the zero-flux sector, one could get away with an effectively free fermion treatment (since the Kitaev Hamiltonian conserves flux)? But now, the novelty here is the initial Neel state, which has superpositions of different flux sectors, making the gauge coupling important.
The paper would benefit greatly from having a discussion on this point in connection to quenches, given the subtle discussions. Specifically, would there be a scenario, given the right choice of initial state, which would be associated with free fermion evolution? If so, what would the expectation be? How exactly does this deviate from the free fermion problem (there is some nice discussion on what the gauge field is doing but it is buried in the text)? And what striking difference does this give rise to in the quench behavior?

3. The other referee brought up an interesting comment on Kitaev’s treatment and relating fermions and spin. I too am a fan of Kitaev’s paper as well as the Majorana basis. I am not sure that the referee comment was addressed though, namely, expressing things in the spin basis is a non-trivial task. I think the specific choice of operators, i.e. S^z correlators, and perhaps the magnetization too, overcomes this problem in getting rid of the string in the Majorana basis? Could the author comment on this in the manuscript?

Requested changes

See Report – Please address all three points

  • validity: high
  • significance: high
  • originality: high
  • clarity: good
  • formatting: good
  • grammar: excellent

Anonymous Report 1 on 2019-7-27 Invited Report


Technically difficult numerical simulations of the Kitaev honeycomb model.

Numerical results that are consistent with what would expect and so I believe the methodolgy is sound and may inspire other studies.


The author is making a claim of exponentially long time scales which are not backed up by numerical evidence or by what is predicted by the underlying effective models.


The author has made a real effort to expand and explain some of the points I was unsure about. In particular I take the point regarding the growth of entanglement and can see that I was in error here. I also think the numerical method the author is using is giving sensible results and so commend the paper on this technical aspect.

However the paper still suffers from problems. The most serious being the claim that there is some exponentially long time scale with respect to thermalisation.

In response to the reply:

“The second question is whether these results are not very trivial since we project onto, quote, “a nearly flat energy band”. The Kitaev model is exactly solvable, yes, but it is not a model of free fermions and to speak of “flat energy bands” is not correct. I guess the referee means that the distribution of energy eigenvalues of the initial Neel state is quite narrow, and indeed, this width scales as Var(E)∼J2xy. However, that still begs the question why the expected slower dynamics are not scaling as some power of Jxy? To reiterate, such an exponential dependence on the model parameters is almost a textbook definition of prethermalization.”

I point out that in the paper the author actully says the observed dynamics can be explained by the Toric code, …

“This prethermal regime is governed by an effective high-temperature toric code”

…so I’m not suite sure what to make of some of the sentences above.

I actually agree with the statement that the Toric Code (TC) Hamiltonian governs the so-called "prethermal" regime, which was what I was attempting to pointing out my previous comments.

Lets examine what should happen. In the TC limit (A phases) Kitaev showed that this band of states it is split on the fourth order into $ -C \sum Q_p$ where $C=J_x^2 J_y^2/16 |J_z|^3$ and the Q’s are plaquettes formed out of effective spins on the z-bonds. With respect to the mass gap which is ~ $2J_z$, I would say this set of states is flat, certainly for the system sizes used here, and especially when we are working in the limit of small $J_x$ and $J_y$.

Irrespective of our differing notions of what flat actually means, if the Neel states projected entirely onto the TC band (which with the sign convention used here is actually the lowest energy band in the model) then the dynamics is governed by the relative splittings $4C$ and we would be guaranteed exact recurrences after the relative dynamical phases of all states completed an integer number of cycles e.g. at $T^*=2 \pi /4 C$. The next non-trivial order is the 6th at which case the relative positions of the gauge configurations play some small role (see Phys. Rev. B 78, 245121 (2008)). These and higher order terms probably prevents exact recurrences (within the TC manifold) from occurring.

The T* recurrences are clearly visible in magnetisation plot of Figure 4, with first occurring at pi/2 x 10^4. This makes it clear also that the so-called exponentially long plateaus are simply a consequence of plotting a sinusoidal (e.g cos (t/F)) (for some large F) on a log scale.

Away from the anisotropic limit the dynamics of the matter sector would play a more prominent role, and I would argue that the sudden initial changes in overlap amplitude are due to the much more rapid phase mixing that occurs to the part of the state that is not projected onto the anti-ferromagnetic band.

With respect to how the Neel state is supported by different gauge sectors the argument in the comments is convincing. Indeed it looks like exactly half of the gauge sectors support the Neel state, where for any of the allowed vortex configurations (those consistent with anti-ferromagnetic configurations) only two of four possible topological sectors are supported.

One important point that occurs here is that (on a torus) there are only 2^(N -3) independent Z_2 gauge configurations in a given topological sector (choice of boundary condition) with the -3 coming from the two topological degrees of freedom and the last because \prod w_p=I. ) While I understand that potentially one can use all u_ij to encode the vortex sector, many of these encodings lead to the same physical scenario. It seems to me that the N_c \approx 2^3N number here seems to be a fairly significant over-counting.

Is it a problem that you haven’t fixed to one gauge convention?

There is a typo in the definition of the BdG Hamiltonian in the appendix. The way it is written is not Hermitian.

Requested changes

The paper still requires a significant rewrite. In particular the paper cannot suggest that the the prethermal regime is governed by the Toric code when in fact the Toric code predicts something else entirely.

  • validity: ok
  • significance: ok
  • originality: high
  • clarity: ok
  • formatting: acceptable
  • grammar: acceptable

Login to report or comment