# Dynamical Quantum Phase Transitions in the Kitaev Honeycomb Model

###### Abstract

The notion of a dynamical quantum phase transition (DQPT) was recently introduced in [Heyl et al., Phys. Rev. Lett. 110, 135704 (2013)] as the non-analytic behavior of the Loschmidt echo at critical times in the thermodynamic limit. In this work the quench dynamics in the ground state sector of the two-dimensional Kitaev honeycomb model are studied regarding the occurrence of DQPTs. For general two-dimensional systems of BCS-type it is demonstrated how the zeros of the Loschmidt echo coalesce to areas in the thermodynamic limit, implying that DQPTs occur as discontinuities in the second derivative. In the Kitaev honeycomb model DQPTs appear after quenches across a phase boundary or within the massless phase. In the 1d limit of the Kitaev honeycomb model it becomes clear that the discontinuity in the higher derivative is intimately related to the higher dimensionality of the non-degenerate model. Moreover, there is a strong connection between the stationary value of the rate function of the Loschmidt echo after long times and the occurrence of DQPTs in this model.

###### pacs:

64.70.Tg, 05.70.Ln, 05.30.Rt## I Introduction

Recent advances in experimental techniques allow to realize closed quantum systems with cold atomic gases in optical traps. Greiner et al. (2002); Kinoshita et al. (2006) These setups are precisely controllable and the unitary time evolution of the systems can be resolved such that the dynamics are experimentally accessible under well-known conditions. Motivated by the new experimental possibilities a lot of theoretical research on the non-equilibrium dynamics of quantum systems has been conducted in the past years. In these theoretical investigations a common protocol for driving a system out of equilibrium is called quantum quench. Considering a parametrised Hamiltonian , where the parameter typically corresponds to some external field strength in the experimental setup, the system is initially assumed to be in equilibrium with regard to some value of the parameter. Then, the parameter is suddenly quenched to a different final value driving the system out of equilibrium and inducing a non-trivial time evolution.

Studying the quench dynamics of a quantum many-body system, Heyl et al. Heyl et al. (2013) pointed out the close formal similarity of the canonical partition function of an equilibrium system, , and the return amplitude

(1) |

of a time-evolved state, suggesting the possibility of critical behavior in the time evolution in analogy to equilibrium phase transitions. It is known that in the thermodynamic limit the zeros of a partition function coalesce to lines in the complex temperature plane and the equilibrium phase transition is marked by the intersection of the zero-line with the real temperature axis.Fisher (1965) Heyl et al. found that in the case of the transverse field Ising model the boundary partition function

(2) |

has zeros in the complex time plane, which accordingly coalesce to lines in the thermodynamic limit. These lines cross the real time axis after quenching the external field across the quantum critical point inducing nonanalyticities in the rate function of the Loschmidt echo

(3) |

at equidistant critical times . Heyl et al. denote this non-analytic behavior at critical times in the thermodynamic limit as a dynamical quantum phase transition (DQPT). They showed that in experiment the DQPT would be observable by measuring the work distribution function of a double quench; in particular, the Loschmidt echo is the probability of performing no work.

These findings triggered further work aiming at a better understanding of the phenomenon. By considering an additional integrability-breaking interaction in the transverse field Ising chain it was demonstrated that DQPTs are not a peculiarity specific to integrable models, but are stable against some non-integrable perturbations.Karrasch and Schuricht (2013); Kriel et al. (2014) Moreover, the signature of DQPTs was found in higher-dimensional systems, namely, in two-dimensional topological insulatorsVajna and Dóra (2015) and effectively infinite dimensions using DMFT.Canovi et al. (2014) It was observed in various cases that DQPTs are not necessarily connected to quenching across a quantum critical point;Fagotti (2013); Canovi et al. (2014); Vajna and Dóra (2014); Andraschko and Sirker (2014) however, there seems to be a strong connection to topological phase transitions.Vajna and Dóra (2015); Budich and Heyl (2015) Canovi et al.Canovi et al. (2014) detected coexisting solutions for so called generalized expectation values in post quench dynamics and, therefore, they introduced the notion of a first order dynamical phase transition. This could be a way to classify dynamical phase transitions. Furthermore, a close connection between the analytic behavior of in the complex plane and its long time limit is conjectured.Heyl and Vojta (2013)

In this work we study quench dynamics in the Kitaev honeycomb model Kitaev (2006) regarding dynamical quantum phase transitions. The model features a rich phase diagram comprising an extended gapless phase, anyonic excitations, and topological order. Moreover, it is a rare example of a Jordan-Wigner-solvable model in two dimensions. Chen and Hu (2007); Chen and Nussinov (2008) As such it has been studied extensively under various aspects. In this paper we restrict the discussion to the dynamics in the ground state sector.

The dynamics of gapped two-dimensional two-band systems were already studied by Vajna and Dóra with focus on a connection between DQPTs and topological phases.Vajna and Dóra (2015) In the presence of a magnetic field the Kitaev model acquires topological order and becomes a system of the same family, albeit being a spin model. In that case we find the behavior in accordance with their results, namely, DQPTs occur after quenches across the boundary between phases with different Chern number. However, the focus of this work lies on quenching between the topologically trivial phases in the absence of a magnetic field. Similar to their results we find DQPTs as discontinuities in the second derivative, which is inherent to the higher dimensionality of the system, and we elaborate on the relevance of the complex zeros of the dynamical partition function in this context. Moreover, we discuss a remarkable observation regarding the long time behavior of the Loschmidt echo. If no DQPTs occur in the post-quench dynamics, then, although the approached stationary state is always an excited state, the long time limit of the Loschmidt echo is given by the fidelity, i.e., the overlap of the initial state with the ground state of the quenched Hamiltonian. This, however, does not hold if the dynamics exhibit DQPTs.

The rest of this paper is organized as follows: in section II the way of solving the model using Jordan-Wigner transformation is sketched and the phase diagram is introduced. Furthermore, the expressions for the dynamical free energy is derived. In section III the zeros of the dynamical partition function in the complex time plane are treated assuming a general BCS-type Hamiltonian, yielding the criteria for the occurrence of DQPTs and the order of the corresponding nonanalyticity. Finally, the zeros of the partition function and the real time evolution are studied explicitly for the Kitaev model in section IV, and two interesting limits are taken into account as well as ramping as an alternative protocol and the quenching with an additional magnetic field.

## Ii The Kitaev honeycomb model

### ii.1 The model

The Kitaev honeycomb model is defined by the Hamiltonian

(4) |

which describes a spin-1/2 system with the spins located on the vertices (labeled by ) of a honeycomb lattice as depicted in Fig. 1.Kitaev (2006) In this paper we assume the lattice spacing to equal unity.

It has been
shownChen and Nussinov (2008) that for the above Hamiltonian one can find a
Jordan-Wigner contour, which after identifying a conserved
operator^{1}^{1}1Here, we choose (cf. Ref.
Chen and Nussinov, 2008), such that the result for the spectrum agrees with
the result in Ref. Kitaev, 2006 and switching to momentum space yields a
BCS-type Hamiltonian

(5) |

with

(6) |

This Hamiltonian can be diagonalised by a Bogoliubov transformation

(7) |

where

(8) |

(see appendix A for details). Plugging the transformation into eq. (5) yields the diagonal Hamiltonian

(9) |

with spectrum

(10) |

The Hamiltonian splits into a sum over -sectors, i.e., a sum over
, where is a subset of the Brillouin-zone such that
, . A possible choice is one half of the
Brillouin-zone, e.g., all with .^{2}^{2}2To be precise one
would have to decide how to deal with the -axis; however, this will
not play any role in the later calculations.

The spectrum has roots at

(11) |

if , where is any permutation of . We will in the following only consider non-negative on the plane. In this section the above result on the gappedness of the spectrum corresponds to a phase diagram as depicted in Fig. 2. The gapless phase at the center of the diagram is surrounded by three distinct gapped phasesKitaev (2006) , and , where , or , respectively.

In the presence of a magnetic field the spin Hamiltonian (4) becomes

(12) |

and the additional term opens a gap also in the -phase. Moreover, the -phase becomes topologically non-trivial with Chern-number , whereas the phases remain trivial with .Kitaev (2006) At there exists a diagonal form of the Hamiltonian even with non-zero magnetic fieldKitaev (2006) and the spectrum reads

(13) |

with , and

(14) |

where . Through a Bogoliubov transformation (see appendix A) this maps to

(15) |

This case will be studied in section IV.6. Before, we will stick to the case without magnetic field, i.e., real valued and .

### ii.2 Post-quench dynamics

In order to study the dynamics in the Kitaev honeycomb model we will consider situations where the system is initially, at , prepared in the ground state of , i.e., . In terms of the free fermion degrees of freedom the initial state is the vacuum: . At the parameter is quenched to its final value , such that for the time-evolved state is . Making use of the Bogoliubov transformation the initial state can be expressed in terms of the final free fermions, which diagonalise :

(16) |

Here,

(17) |

and the normalization constant

(18) |

were introduced. A more detailed derivation is given in appendix A.

For the sake of brevity and lucidity we will in the following refrain from dragging along the dependencies on and explicitly, i.e., identify and .

Using (16) to compute the dynamical partition function we obtain

(19) | ||||

(20) |

The dynamical partition function has large deviation form , where is the system size. Thus, in the thermodynamic limit only the rate function, or dynamical free energy density,

(21) |

is well defined.

## Iii Zeros of the partition function and critical times

### iii.1 General aspects

From the study of equilibrium phase transitions it is known that a very insightful approach is to consider the zeros of the partition function in the complex temperature or complex magnetization plane, respectively.Yang and Lee (1952); Fisher (1965); Bena et al. (2005) In the thermodynamic limit the zeros of the partition function coalesce to lines or areas in the complex plane, which mark the critical points when approaching the real temperature (magnetization) axis. Analogous reasoning has proven useful in the study of dynamical quantum phase transitions.Heyl et al. (2013); Vajna and Dóra (2015)

In particular, an interesting analogy to electrodynamics allows to characterize the dynamical phase transition through the density of zeros of the dynamical partition function in the complex time plane. The starting point is the observation that the dynamical partition function (19) is an entire function of and can as such, according to the Weierstrass factorization theorem, be written as

(22) |

where is some discrete index set, are the zeros, and is an entire function Heyl et al. (2013). With this, the dynamical free energy density reads

(23) |

From this expression it becomes clear, that any non-analytic behavior of the dynamical free energy can only occur at or in the vicinity of the zeros of the dynamical partition function . Since we are interested in nonanalyticities, we will in the following ignore the contribution of and only consider the singular part

(24) |

In the thermodynamic limit the sum becomes an integral over some continuous variable , where is a region corresponding to the previously used index set , and the zeros become a function of this variable , such that

(25) |

A transformation of the integration variable yields

(26) |

where the Jacobian determinant can be interpreted as the density of zeros in the complex plane.van Saarloos and Kurtze (1984) Moreover, setting for allows to extend the integration domain to the full complex plane. We will now discuss the real part

(27) |

We will later see that the Loschmidt echo on the real time axis is directly given by . For with is the Green’s function of the Laplacian , i.e.,

(28) |

In other words, the real part of the dynamical free energy density can be interpreted as the electrostatic potential produced by a charge density in two dimensions and the question of the behavior of the free energy at critical points becomes the question of the behavior of the electrostatic potential at surfaces. If the zeros form lines in the complex plane, this allows to deduce the order of the phase transition directly from the density of zeros at or in the vicinity of the physically relevant . Bena et al. (2005)

Although the zeros can form areas in the complex plane these areas do not cover the physical axis in the case of equilibrium phase transitions. However, as we will see in the following section, this is possible in the case of dynamical phase transitions. Thus, consider the situation as depicted in Fig. 3. The density of zeros is given by and in area and , respectively, and at the boundary there is a discontinuous change in the density of zeros. Assume the electric potentials , , solve the Laplace equation (28) with the corresponding density . With this a global solution is

(29) |

and we demand continuity of at the boundary. Let us now focus on the behavior of at the intersection of the boundary with the real time axis. If we choose the --coordinate frame as indicated in Fig. 3, the continuity of implies that on the boundary

(30) |

Since we are interested in the behavior of in the real time axis, we transform to --coordinates,

(31) |

With the Laplace equation (28) this yields

(32) |

and, consequently,

(33) |

This means, if an area of zeros of the partition function overlaps the real time axis, the second derivative of the free energy is discontinuous.

### iii.2 2d BCS-type models

In the Kitaev honeycomb model (and general BCS-type models) the partition function is given by (20), i.e., the zeros in the complex plane are given by

(34) |

At this point it becomes obvious, that in the thermodynamic limit the double product over and in (20) leads to dense areas of zeros in the complex plane, since, generally, . These areas of zeros cover parts of the real time axis () if , i.e., if

(35) |

Dubbing the isoline , there are intervals

(36) |

on the real time axis, which are covered by areas of zeros. If the spectrum is gapped, the beginnings and end points of two consecutive intervals and are equidistant with and , where are momenta minimizing/maximizing on the domain given by . The length of the single intervals increases linearly with . However, if the spectrum is gapless, all those intervals extend to infinity.

The condition (35) allows for a physical interpretation; namely, the occurrence of DQPTs is through a continuity argument related to non-thermal mode occupation.Heyl et al. (2013) In BCS-type models the mode occupation is given by

(37) |

This
means for all modes , where the condition (35) is satisfied,
the mode occupation is . Let us assume that for
any two points in there exists a path connecting both points, along which
is continuous,
^{3}^{3}3Assuming conventional continuity is too strong in this case, since in
the Kitaev model the mode occupation number is not necessarily
continuous after quenching (see Fig. 8).
and the existence of modes with . Both assumptions
should be true for physically relevant models and were found to hold in all cases
considered in the Kitaev model. Then, we can set up the following chain of
consequences: through the continuity condition, the existence of non-thermally
occupied modes with implies the existence of modes with
. This in turn is equivalent to the fulfilling of the
condition (35), which implies the occurrence of DQPTs in the time
evolution. The mode occupation
is non-thermal in the sense that it cannot be realized by Fermi-Dirac
statistics with positive temperature.

An equivalent formulation of condition (35) is

(38) |

where is the quench parameter of the BCS-type Hamiltonian, (see appendix B). From this it becomes clear, that after quenching to a gapless phase there are zeros of the dynamical partition function on the real time axis, since . In the mode occupation picture this can be interpreted as follows: when quenching to a gapless phase, excitations cost no energy; thus, any quench produces inverted mode occupation.

As discussed in the previous section, areas of zeros covering the real time axis result in jumps in the second time derivative of the dynamical free energy density if there is a jump in the density of zeros. Eq. (34) gives a “layer” of zeros for every . Therefore, our total density of zeros is a sum of the densities of the individual “layers”, . The single layer densities are given the Jacobi determinant of the change of variables ,van Saarloos and Kurtze (1984)

(39) |

At this point a more technical view of the zeros of the partition function is useful: the zeros correspond to intersections of the isolines

(40) |

in the momentum plane. This means that the density of zeros diverges at the boundary, since there . Thus, when approaching the boundary of an interval from the inside of the interval, the second time derivative of will diverge.

## Iv Dynamical phase transitions in the Kitaev honeycomb model

### iv.1 Zeros of the dynamical partition function in the Kitaev model

In the Kitaev model not only quenches to the massless phase create inverted mode occupation. Also quenches across phase boundaries with final parameter in a massive phase induce critical points in the real time evolution. It is physically reasonable to assume, that the mode occupation number is sufficiently well behaved, namely, that for any two , there exists a path with and such that , , is continuous. We found this to be true for all considered cases. Under this prerequisite, the existence of a fully occupied mode , , implies that somewhere, because . corresponds to and this happens when

(41) |

One possibility to fulfill this is

(42) |

Now, consider a quench ending in the -phase () and . Then and . We find that at this point both quenches, starting from another massive phase,

(43) |

and from the massless phase,

(44) |

lead to non-analytic behavior because (42) is fulfilled in both cases. The same can be shown for quenches ending in the other massive phases, only needs to be chosen appropriately. This shows that in the Kitaev model occupation inversion is produced by quenches within the massless phase or quenches crossing phase boundaries.

Figure 4 displays locations of the zeros of the Loschmidt echo in the complex plane given by eq. (34) for two quenches, one within the phase and one from the phase to the massless phase. Both panels include a phase diagram with an arrow indicating the quench parameters . The numerical values for the parameters for this figure and all following figures are listed in Tab. 1 in the appendix. The zeros do indeed form areas, which are restricted to the left half-plane for the quench within the massive phase but cover parts of the real time (imaginary ) axis when and lie in different phases.

### iv.2 Real time evolution

On the real time axis the rate function of the Loschmidt echo reads

(45) |

and the time derivative is

(46) |

Figure 5 shows the time evolution of the rate function and its time derivative for various quenches obtained by numerical evaluation of the corresponding integrals. The gray-shaded areas in the plots indicate the intervals (cf. eq. (36)) of vanishing partition function. The results exhibit the properties expected from the previous considerations. The rate function is smooth for quenches within the gapped phases; however, nonanalyticities occur when phase boundaries are crossed in a quench or after a quench within the gapless phase. The beginning and end points of the critical intervals are equidistant, respectively, and for quenches ending in the massless phase the intervals extend to infinity. As anticipated, nonanalyticities only show up at the boundaries of the critical intervals. Moreover, the nonanalyticities emerge as discontinuities of , i.e., kinks in .

Note that the two plots in panels a) and c) show the time evolution of the rate function after the two quenches for which Fig. 4 shows the location of the zeros of the partition function.

### iv.3 The 1d limit

In the limit for any the 2d Kitaev model (4) degenerates and becomes a set of separate 1d spin chains. Let us consider the case . The vanishing of one of the other two parameters will give the same result due to the threefold symmetry. For the condition for nonanalyticities (35) is fulfilled at with

(47) |

Along this line also the spectrum is constant,

(48) |

Thus, the critical intervals defined in eq. (36) become critical points

(49) |

on the real time axis. Figure 6 shows the rate functions for two different quenches with . The quench in Fig. 6a does not cross a phase boundary and therefore the rate function is analytic. However, in Fig. 6b and lie in different phases, and according to eq. (49) there are critical times at which the rate function becomes singular. In particular, these singularities are discontinuities in the first time derivative of the rate function, not in the second as in the genuinely two-dimensional cases. This observation underlines the fact that the continuity of the first derivative is inherent to the higher dimensionality of the non-degenerate Kitaev model.

### iv.4 The long time limit

In a recent workHeyl and Vojta (2013) it was stated that if the Loschmidt echo supports analytic continuation, then the long real time limit and the long imaginary time limit coincide. In the context of dynamical phase transitions this gives rise to the conjecture, that the occurrence of DQPTs is closely related to the long time behaviour of . Dynamical quantum phase transitions occur if the zeros of the dynamical partition function cross the real time axis. This means, , the rate function of the Loschmidt echo, is non-analytic in the half plane; therefore, the long imaginary-time limit and the long real time limit do not necessarily have to coincide.

On the imaginary time axis , with the eigenbasis of the quenched Hamiltonian , corresponding energies , and ,

(50) |

Now, shifting the energy such that and assuming the ground state to be non-degenerate,

(51) |

where is the fidelity. Thereby, the Loschmidt echo is connected to the fidelity in the large imaginary time limit.

As the Loschmidt echo is not well defined in the thermodynamic limit, one should rather formulate eq. (51) in terms of the rate function:

(52) |

According to the previous considerations this yields