Defect production and quench dynamics in the three-dimensional Kitaev model

TL;DR

This paper investigates the quench dynamics of the 3D Kitaev model, revealing how defect density and correlations depend on quench rate and identifying key factors influencing non-adiabatic defect production.

Contribution

It provides an analytical and numerical study of defect production in the 3D Kitaev model, highlighting the role of level crossings over the energy spectrum's gapless contour.

Findings

Defect density scales linearly with quench rate for slow quenches.

Defect correlation decays as τ^{-1} e^{-A/τ}, linking it to entropy.

Zeros of level coupling determine non-adiabatic defect production.

Abstract

We study the quench dynamics of the three-dimensional Kitaev (spin) model under a linear drive using both exact numerical calculations and analytical "independent crossing approximation". Unlike the two-dimensional Kitaev model, the three-dimensional Kitaev model reduces to a multilevel Landau-Zener problem for each momentum. We show that for the slow quench, the defect density is proportional to the quench rate $1/\tau$. We find that the zeros of the relevant coupling between the levels determine the non-adiabatic condition for the production of defects. The contour on which the energy spectrum becomes gapless does not play an active role. The asymptotic behavior of the defect density crucially depends on the way the system reaches the non-adiabatic regime during the quenching process. We analytically show that defect correlation varies as $\tau^{-1} e^{-A/\tau}$, where $A$ is a constant independent of $\tau$. For the slow quench, the qualitative dependence of the entropy (produced during the quenching process) on the quench time is the same as that of the defect correlation, indicating a close connection between the defect correlation and the entropy content of the final state. Possible experimental realization of such quench dynamics is also described briefly.

Figures

Equations

$$H=-J_x\sum _{x-link}{\sigma }_i^x{\sigma }_j^x-J_y\sum _{y-link}{\sigma }_i^y{\sigma }_j^y-J_z\sum _{z-link}{\sigma }_i^z{\sigma }_j^z.$$

$$\displaystyle H=i\sum_{\mathbf{r}}\Big{[}\sum^{1,2}_{j}J_{z}c_{ja}(\mathbf{r})c_{jb}(\mathbf{r})+J_{x}c_{1a}(\mathbf{r})c_{2b}(\mathbf{r_{3}})+$$

$$\displaystyle J_{y}c_{1a}(\mathbf{r})c_{2b}(\mathbf{r_{13}})+J_{y}c_{1b}(\mathbf{r})c_{2a}(\mathbf{r_{2}})+J_{x}c_{1b}(\mathbf{r})c_{2a}(\mathbf{r})\Big{]},~{}~{}$$

$$h(k)$$

$$\tilde{h}(k)$$

$$\displaystyle\mathcal{U}_{k}=\left(\begin{array}[]{cccc}\epsilon_{2k}&0&\gamma_{k}&-g_{k}\\ 0&\epsilon_{1k}&g_{k}&-\gamma_{k}\\ \gamma_{k}&-g_{k}&-\epsilon_{1k}&0\\ g_{k}&-\gamma_{k}&0&-\epsilon_{2k}\end{array}\right),$$

$$\displaystyle\tilde{\mathcal{U}}_{k}=\left(\begin{array}[]{cccccc}4J_{3}&\gamma_{k}&g_{k}&-g_{k}&-\gamma_{k}&0\\ \gamma_{k}&0&0&0&0&-\gamma_{k}\\ +g_{k}&0&-2|\Delta_{2k}|&0&0&-g_{k}\\ -g_{k}&0&0&2|\Delta_{2k}|&0&+g_{k}\\ -\gamma_{k}&0&0&0&0&\gamma_{k}\\ 0&-\gamma_{k}&-g_{k}&g_{k}&\gamma_{k}&-4J_{3}\end{array}\right);~{}~{}$$

$$\hat{U}(k,t)=\left(\hat{I}+i\hat{h}(k,t)\frac{dt}{2}\right){\left(\hat{I}-i\hat{h}(k,t)\frac{dt}{2}\right)}^{-1},$$

$$p_k^{(6)}$$

$$n_d^{(i)}=\frac{1}{{\Omega }_{HBZ}}{\int }_{HBZ}{d}^{3}k{p}_k^{(i)},i=4,6,$$

$$|{\Delta }_{1\mathbf{k}}{|}^2$$

$$\tan (k_x)\tan (k_z)$$

$$|{\Delta }_{1\mathbf{k}}{|}^2$$

$$n_d=\frac{\pi }{{\Omega }_{HBZ}}\left(\frac{1}{(J\tau )\sqrt{||f_{ij}(\alpha )||}}\right),$$

$${\mathcal{E}}_k^{1-4}=\pm \sqrt{8J_z^2+2(|{\Delta }_{1k}{|}^2+|{\Delta }_{2k}{|}^2)\pm {\Delta }_k},$$

$${\left(4J_z^2+|{\Delta }_{1k}{|}^2-|{\Delta }_{2k}{|}^2\right)}^2$$

$$+4|{\Delta }_{1k}{|}^2|{\Delta }_{2k}{|}^2{\sin }^2{\theta }_{2k}=0.$$

$${\mathcal{O}}^{3D}(\mathbf{r})={\sigma }_{\mathbf{R}}^z{\sigma }_{\mathbf{R}+\mathbf{r}}^z,$$

$${\mathcal{C}}^{3D}(\mathbf{r})=i{c}_{1a}(\mathbf{R})c_{1b}(\mathbf{R}+\mathbf{r}).$$

$$⟨{\mathcal{C}}^{3D}(\mathbf{r})⟩=\pm \frac{1}{2N}\sum _{k\in HBZ}\cos (\mathbf{k}\cdot \mathbf{r}),$$

$$⟨{\mathcal{C}}^{3D}(\mathbf{r})⟩=-{\delta }_{\mathbf{r},0}+\frac{1}{{\Omega }_{HBZ}}\int d^{3}k{p}_k\cos (\mathbf{k}\cdot \mathbf{r}),$$

$$⟨{\mathcal{C}}^{3D}(\mathbf{r})⟩=\frac{\pi }{{\Omega }_{HBZ}}\left(\frac{\cos \left[k_x^0(N_x+N_y)+k_z^0{N}_z\right]}{(J\tau )\sqrt{||f_{ij}(\alpha )||}}\right)\times$$

$$\exp \left[-\frac{(N_x+N_y{)}^2{f}_{zz}+4N_z^2{f}_{xx}-4(N_x+N_y)N_z{f}_{xz}}{4\pi J\tau ||f_{ij}(\alpha )||}\right].$$

$$⟨{\mathcal{C}}^{3D}(\mathbf{r})⟩\sim {\tau }^{-1}\exp (-A{\tau }^{-1}),$$

$$S=-{\int }_{HBZ}\frac{d^{3}k}{{\Omega }_{HBZ}}\sum _i{A}_i(k,\infty )\ln [A_i(k,\infty )],$$

$$|{\Psi }_4(k,t){⟩}_{ICA}$$

$$S^{(ICA)}=-\frac{1}{{\Omega }_{HBZ}}{\int }_{HBZ}{d}^{3}k[\sum _i{P}_{0i}\ln P_{0i}],i=0,2,3.$$

$$\left(\begin{array}{c}{\psi }_k\\ {\varphi }_k\end{array}\right)$$

$$u_k=\left(\begin{array}{cc}-e^{-i{\theta }_{2k}}& e^{-i{\theta }_{2k}}\\ 1& 1\end{array}\right),e^{i{\theta }_{2k}}=\frac{{\Delta }_{2k}}{|{\Delta }_{2k}|}.$$

$$|\tilde{\varnothing }⟩$$

$$|\mathcal{I}⟩$$

$$|\tilde{0}⟩={\eta }_{0k}^{†}{\eta }_{1k}^{†}|\tilde{\varnothing }⟩,{\widetilde{ϵ}}_0=4J_3,$$

Full text

I Introduction

Properties of physical systems near a quantum phase transition are determined by the change of the symmetry of the ground state. In particular, at low temperatures, a quantum phase transition is observed when the inherent quantum fluctuations win over the thermal fluctuations and determine the ground state properties Sachdev (2008, 2011); Carr (2011). In this respect, it is pertinent to ask what happens at absolute zero temperature, when the parameters of a given Hamiltonian (or physical system) are driven or varied in time and the system is taken across the different phases. Recently such driven quantum systems have attracted a lot of interest Polkovnikov et al. (2011); Dziarmaga (2010). On a more fundamental note, it is appropriate to ask if the quantum system always remains in its instantaneous ground state or fail to do so due to the diverging length and time scales, when taken across the critical region very slowly (by varying the relevant system parameters linearly with time) Sachdev (2008, 2011); Carr (2011); Polkovnikov et al. (2011); Dziarmaga (2010); Heyl (2018); Sondhi et al. (1997). Consequently, as the quantum critical point (QCP) is approached the system gradually enters into a non-adiabatic regime by failing to remain in its instantaneous ground state. As a result, excitations are created which are also called the defects Polkovnikov et al. (2011); Dziarmaga (2010).

In a generic second-order phase transition, the defect density, produced during quench across the QCP, scales with the quench time as $\tau^{-\frac{d\nu}{(z\nu+1)}}$, $d$ being the dimension of the system, $\nu$ and $z$ are the correlation length exponent and the dynamical critical exponent respectively Kibble (1976, 1980); Zurek (1985, 1993, 1996); Dziarmaga (2010); Polkovnikov et al. (2011); Damski (2005). It is found explicitly that in low dimensions the response of the system to the slow changes of the Hamiltonian parameters could be non-analytic and non-adiabatic Polkovnikov (2005). It has been shown that for a sufficiently slow quench at a rate $1/\tau$ $(<<1)$, the density of the above mentioned defect states scales with the quench time as, $n_{d}\sim\tau^{-(\frac{m\nu}{z\nu+1})}$, when the process of quenching takes the system across a $d-m$ dimensional gapless (critical) hyper-surface instead of a QCP Polkovnikov and Gritsev (2008); Dziarmaga (2010); Sen et al. (2008); Mondal et al. (2008). Such a scaling behavior has been studied in the two-dimensional (2D) Kitaev model as well as in several other exactly solvable models Dziarmaga (2010); Sen et al. (2008); Mondal et al. (2008).

The Kitaev model is a rare example of a fully anisotropic spin model which is exactly solvable Kitaev (2006). The presence of Kitaev-like interactions has been found in a class of materials, viz., ‘iridate’ systems Jackeli and Khaliullin (2009); Singh and Gegenwart (2010); Singh et al. (2012); Plumb et al. (2014). Several aspects of the Kitaev model have been explored theoretically such as, topological order, fractionalization of spins, quantum spin liquid, entanglement generation, as well as in the context of quantum computation Nayak et al. (2008); Balents (2010); Hermanns et al. (2018); Anderson (1973); Moessner and Sondhi (2001); Wen (2004); Lhuillier and Misguich (2011); Baskaran et al. (2007); Mandal et al. (2012, 2016). A three-dimensional (3D) version of the Kitaev model has been proposed and studied theoretically, and has been realized experimentally in layered ‘iridate’ compounds Mandal and Surendran (2009); Okamoto et al. (2007); Modic et al. (2014); Takayama et al. (2015). The 3D Kitaev model contains four sub-lattices in an unit cell and its energy spectrum vanishes on a contour. Furthermore, it supports fermionic and bosonic non-local excitations as well Mandal and Surendran (2014).

Motivated by the above studies, we investigate the quench dynamics in the spin-1/2 3D Kitaev model on a hyper-honeycomb lattice. Due to the presence of a four-sublattice structure, the model, in its ground state flux configuration, exhibits a pair of negative energy and positive energy states for each momentum. This enables us to explore a multi-level Landau-Zener (L-Z) problem in the quench dynamics of the 3D Kitaev model. By filling up both the negative energy states, we end up in a six-level L-Z problem. On the other hand, if we start with the most negative energy state being filled, the model reduces to a four-level L-Z problem. Although such a multi-level L-Z problem (both the four-level Sinitsyn (2015a, 2013, b) and beyond Usuki (1997); Demkov and Ostrovsky (2000)) has been studied using model Hamiltonians, however, most of them are either of mathematical interest or lack a direct connection to the real materials. The question of the exact solvability of multi-level L-Z problems has also attracted a lot of interest Chernyak et al. (2020, 2018). It is therefore intriguing to study the quench dynamics of the 3D Kitaev model (in the multi-level set-up) which is both exactly solvable and experimentally realizable in a class of materials mentioned earlier.

In this paper, we investigate the defect density, defect correlation, and entropy in the final state of the quench dynamics in the 3D Kitaev model. We define the defect correlation as the two-point fermionic correlation in the defect state (time-evolved initial state) Polkovnikov (2005); Sen et al. (2008); Mondal et al. (2008). Furthermore, quenching pushes the system to make transitions to other energy levels via non-adiabatic points. Moreover, we calculate the von-Neumann entropy of the final state Nielsen and Chuang (2010); von Neumann and BEYER (2018); Cherng and Levitov (2006). Interestingly enough, we find that the 3D Kitaev model offers a solvable physical system where the semi-classical “Independent Crossing Approximation(ICA)” can be applied to obtain an analytical form for the defect density, thereby enlarging our scope of understanding of the multi-level L-Z dynamics. Remarkably, we observe that in the multi-level L-Z dynamics of 3D Kitaev model the non-adiabatic condition between the relevant levels is not determined entirely by the spectrum of the model, rather by the coupling between the adiabatic levels. This might lead to a scaling of the defect density beyond the conventional two-level L-Z dynamics. We further investigate the defect correlation and obtain an analytical expression for it which agrees with the corresponding numerical results. The entropy is found to be generated due to the non-local correlation developed during the quench.

The paper is organized as follows. In section II, we introduce the 3D Kitaev model as conceived in Ref Mandal and Surendran (2009). We explain in detail the methods to transform the Hamiltonian in momentum space into the one that is suitable for the quench dynamics study. Here we elaborate our scheme to reduce the 3D Kitaev model to a multi-level L-Z problem (both four and six levels) and describe the numerical scheme to solve the time dependent Schrodinger equations. In Sec. III, we present the analytical formula for the defect density within the semi-classical “Independent Crossing Approximation (ICA)”. Sections IV, V, and VI describe the results corresponding to the defect density, defect correlation and the entropy of the final state respectively, obtained numerically using the method outlined in II.3, and analytically as well using ICA (whenever applicable). We finally conclude and discuss our results in section VII. In the Appendices B and C we briefly discuss the issue of dynamical phases in the transition probabilities associated with semi-classical trajectories and justify the ICA. The other Appendices describe various relevant details related to our calculations.

II The Model and Formalism

The three-dimensional (3D) Kitaev model is defined on a hyper-honeycomb lattice with lattice coordination number ‘3’ (see Fig. 1 (a)). Each spin interacts with its three nearest neighbors through three different links, viz., $x$-link, $y$-link, and $z$-link which only exhibit interaction between the corresponding spin components Mandal and Surendran (2009). The spin Hamiltonian on the hyper-honeycomb lattice has the following form,

[TABLE]

In the above equation, ‘$i,j$’ represents a pair of nearest neighbor spins and the sum is over each type of bonds. The unit cell contains four sites with the basis vectors given by, $\mathbf{a}_{1}=2\hat{\mathbf{x}}$, $\mathbf{a}_{2}=2\hat{\mathbf{y}}$, and $\mathbf{a}_{3}=(\hat{\mathbf{x}}+\hat{\mathbf{y}}+2\hat{\mathbf{z}})$ Mandal and Surendran (2009).

Following Kitaev’s original prescription for solving the 2D Kitaev model, the above Hamiltonian can be fermionized to reduce to a non-interacting Majorana Fermion hopping problem in the presence of conserved $\mathbb{Z}_{2}$ gauge fields defined on every links/bonds. These $\mathbb{Z}_{2}$ gauge fields can take values $\pm 1$. For each configuration of these $\mathbb{Z}_{2}$ gauge fields, one obtains a fermionic spectrum. It is found that Mandal and Surendran (2009), the ground state sector corresponds to the configuration where the product of $\mathbb{Z}_{2}$ gauge fields over an elementary plaquette is one. We chose the $\mathbb{Z}_{2}$ gauge field to be equal to ‘$+1$’ for each bond, and this enables us to perform a Fourier transform and work in the momentum space. In the ground state sector (where all $\mathbb{Z}_{2}$ gauge fields are uniformly ‘$+1$’) the effective Majorana fermion Hamiltonian is given by,

[TABLE]

where $\mathbf{r}$ is a position vector of a lattice point Mandal and Surendran (2009) and $\mathbf{r}_{3}=\mathbf{r}-\mathbf{a}_{3},\mathbf{r}_{2}=\mathbf{r}-\mathbf{a}_{2}$ and $\mathbf{r}_{13}=\mathbf{r}+\mathbf{a}_{1}-\mathbf{a}_{3}$. Owing to the bipartite nature of the lattice, in the above Hamiltonian, we conveniently introduce two indices to label the sites within a unit cell; $\mu=1,2$, denote the dimer to which a site belongs and $\alpha=a,b$, denote the sub-lattice indices Mandal and Surendran (2009). The spectrum of the Hamiltonian in (II) is found to be gapless in the middle white region as shown in Fig. 1 (b).

In our quench study, we take $J_{z}$ to be increasing linearly with time $t$, i.e., $J_{z}(t)=Jt/\tau$ at a rate $1/\tau$. To implement the quench study in the L-Z set-up, we rewrite the Hamiltonian in (II) by introducing the complex fermions $\psi$ and $\phi$ by regrouping the two Majorana fermions at two z-bonds of the unit cell as $\psi_{i}=c_{i,1a}+ic_{i,1b},~{}~{}~{}\phi_{i}=c_{i,2a}+ic_{i,2b}$. A subsequent Fourier transform of (II) leads us to, $H=\sum_{k}^{{}^{\prime}}\Psi_{k}^{\dagger}h(k)\Psi_{k}$, where ‘$k$’ belongs to half Brillouin zone (HBZ, denoted by the primed summation). Here,

[TABLE]

and $\Psi_{k}=\begin{pmatrix}\psi_{k}^{\dagger}&\phi_{k}^{\dagger}&\psi_{-k}&\phi_{-k}\end{pmatrix}^{\dagger}$, where $\Delta_{1k}=e^{i\mathbf{k\cdot a_{3}}}\delta_{1k}+\delta_{2k}$, and $\Delta_{2k}=e^{-i\mathbf{k\cdot a_{3}}}\delta_{1k}^{*}-\delta_{2k}^{*}$, and $\delta_{j,k}=J_{x}+J_{y}e^{i(-1)^{j}k_{j}}$ with $j=1,2$. Diagonalizing the upper-left and lower-right $2\times 2$ blocks we get a couple of eigenvalues, $\epsilon_{1k}=2J_{z}+|\Delta_{2k}|$ and $\epsilon_{2k}=2J_{z}-|\Delta_{2k}|$ corresponding to the upper-left block, and another couple $-\epsilon_{1k}$ and $-\epsilon_{2k}$ corresponding to the lower-right block respectively. Using the unitary transformation corresponding to the above mentioned diagonalization (see Appendix A for explicit form), we re-write (3) as,

[TABLE]

where $g_{k}=-i|\Delta_{1k}|\sin\theta_{12k}$ and $\gamma_{k}=-|\Delta_{1k}|\cos\theta_{12k}$ with $e^{i\theta_{12k}}=\frac{\Delta_{1k}\Delta_{2k}^{*}}{|\Delta_{1k}||\Delta_{2k}|}$, and $\epsilon_{1k}=2J_{z}+|\Delta_{2k}|$ and $\epsilon_{2k}=2J_{z}-|\Delta_{2k}|$. We denote $\alpha_{i,k}$($i=0,1,2,3$) as the annihilation operators corresponding to the four states in the rotated basis. The relations between the $\alpha_{ik}$ to the $\psi_{k}$ and $\phi_{k}$ are given in Appendix A. We apply the quench protocol by keeping $J$ fixed for a set of values of $\alpha=\tan^{-1}\frac{J_{y}}{J_{x}}$ quantifying the ratio of $J_{x}$ and $J_{y}$ .

The Hamiltonian corresponding to (4) offers two ways of studying the quench dynamics depending on if we start with the initial state as the lowest negative energy state being occupied or start by filling up both the negative energy states. In the former case, we end up with a four-level L-Z problem, and in the latter, we find a six-level L-Z problem. It is worthwhile to emphasize that the ground state of the six-level problem belongs the true ground state of the spin-model at $t=-\infty$. Therefore, we denote the lowest energy state at $t=-\infty$ corresponding to the six-level problem as the ground state and the same corresponding to the four-level problem as the initial state. However, we shall henceforth use them synonymously and expect their meaning to be clear in the context. Below we discuss both of these two possibilities one by one. Furthermore, in connection to the original spin model corresponding to (1), this physically means that for the four-level L-Z problem, the initial state is the ferromagnetic state of the lower $z$-dimer of the unit cell. On the other hand, for the six-level problem, the initial state is a ferromagnetic state of both the $z$-dimers of the unit cell.

II.1 As a four-level problem

In the limit of $t\rightarrow\pm\infty$, the dominating terms in (4) are the diagonal ones where we retain a very small non-zero value of $|\Delta_{2k}|$, given the fact that $J_{x(y)}$ have infinitesimal but non-zero values. This helps us to avoid the obvious degeneracy (which otherwise originates from the fact that two of the diagonal elements of (4) become equal) in the limit $t\rightarrow\pm\infty$, in the same spirit of removing the ground state degeneracy of a spin-$S$ Ferromagnetic Heisenberg model with an infinitesimally small symmetry breaking field. In the limit of $t\rightarrow-\infty$ the initial state is gapped, and the corresponding eigenvalue is given by, $\epsilon_{2k}=-2J_{z}-|\Delta_{2k}|$ with $|\Delta_{2k}|$ very small (but non-zero) so that the diabatic levels are well separated at $t\pm\infty$, as plotted in Fig.2(a). According to our convention, at $t\rightarrow-\infty$, the eigenstates corresponding to $\epsilon_{2k},\epsilon_{1k},-\epsilon_{1k}$ and $-\epsilon_{2k}$ are denoted by $|0\rangle_{k},|1\rangle_{k},|2\rangle_{k}$, and $|3\rangle_{k}$ respectively. The diabatic ground state/initial state in this time limit is given by, $|G_{-\infty}\rangle=\begin{pmatrix}1&0&0&0\end{pmatrix}^{\dagger}$ whereas, in the limit $t\rightarrow+\infty$ the diabatic initial state is $|G^{{}^{\prime}}_{+\infty}\rangle=\begin{pmatrix}0&0&1&0\end{pmatrix}^{\dagger}$. The corresponding adiabatic levels are plotted in Fig.2(b).

The quench dynamics of various model systems have been studied in the past, and these, in essence, reduce to a two-level L-Z problem Polkovnikov et al. (2011); Polkovnikov (2005); Sen et al. (2008); Mondal et al. (2008) which has been analytically examined. However, the system we are interested in turns into a multi-level L-Z model for each momentum $\mathbf{k}$, which limits the possibility of a complete analytical solution. In these circumstances, the numerical analysis of the evolution of the system under quench turns out to be the most convenient to determine the long time behavior of the system. We need to solve a set of coupled first-order differential equations corresponding to the time-evolution, which is given by, $i\frac{d}{dt}\Phi_{4,k}=\mathcal{U}_{k}\Phi_{4,k}$, where $\Phi^{\dagger}_{4,k}=(|0,t\rangle,|1,t\rangle,|2,t\rangle,|3,t\rangle)$, and the $4\times 4$ matrix $\mathcal{U}_{k}$ is given below,

[TABLE]

is the time evolution operator which is obtained from (4) by an unitary transformation, U_{k}=\left(\begin{array}[]{cc}\sigma_{x}&0\\ 0&I_{2\times 2}\end{array}\right).

II.2 As a six-level problem

We now describe the quench dynamics when the initial state is constructed by filling up both the negative energy states corresponding to (4) at $t\rightarrow-\infty$. The possibility of transition from the ground state to any state having two particles results into a six-level L-Z problem. To this end one needs to define the two particle sector of the Hilbert space and it is given in Appendix A. Time evolution is governed by the usual Schrodinger equation $i\frac{d}{dt}\Phi_{6,k}(t)=\tilde{\mathcal{U}}_{k}\Phi_{6,k}(t)$, where $\tilde{\mathcal{U}}_{k}$ is given by,

[TABLE]

see Appendix A for further details.

II.3 Numerical scheme

Numerical evolution of the time-dependent quantum states needs a careful algorithm to avoid the large error accumulation, which might result in producing a non-physical outcome. We use a numerical algorithm that captures the time evolution starting from an arbitrarily large time in the past to an arbitrarily large time in the future and proceeds in discrete but small steps $dt$. The unitary operator describing the time evolution is given by,

[TABLE]

where $\hat{h}(k,t)$ is the Hamiltonian matrix governing the time evolution which is (9) in the four-level setup and is (16) in the six-level one, with the quench protocol mentioned earlier. The identity matrix $\hat{I}$ is a $4\times 4$ or a $6\times 6$ one depending on the set-up. Such a time evolution operator has been used to study a four-level L-Z model earlier in the context of the model Hamiltonian Sinitsyn (2015a, 2013). The evolution operator (17) is equivalent to a true evolution operator in the limit of $dt\rightarrow 0$, and reduces the error to $\frac{(dt)^{2}}{4}$ (reduction by a factor of 4) at $\mathcal{O}(dt^{2})$. In this paper, we calculate the following quantities, $\langle 0|\hat{U}(k,t)|0\rangle$ and $\langle f|\hat{U}(k,t)|0\rangle$, where $|f\rangle$ are the final states except the post-quench ground state (i.e., the ground state corresponding to $t=\infty$). The quantity $P_{00}=\langle 0|\hat{U}(k,t)|0\rangle$ describes the probability of the system to remain in the initial state. Using (9) (for the four-level problem), (16) (for the six-level problem), and (17) we determine the defect density $n_{d}=\frac{1}{\Omega_{HBZ}}\int_{HBZ}d^{3}k~{}P_{00}$, with $\Omega_{HBZ}=\pi^{3}/2$ being the volume of the HBZ.

Furthermore, we aim to determine the dependence of the defect density on the quench rate $1/\tau$ (or equivalently $1/J\tau$). Physically, the defect density refers to the number density of the created excitations Damski (2005); Polkovnikov (2005); Polkovnikov and Gritsev (2008). These excitations are generated once the system is taken across a QCP where the system fails to remain in the adiabatic regime. As explained above, depending on the choice of the initial state the 3D-Kitaev model reduces to a four-level or a six-level L-Z problem for each momentum $\mathbf{k}$.

III Analytical scheme- the independent crossing approximation

We use the $independent~{}crossing~{}approximation$ ($\rm{ICA}$) Sinitsyn (2015a) to analytically solve the multi-level L-Z problem. This treats the multi-level L-Z problem as a collection of independent two-level problems, each two-level problem being solved without considering the influence of the rest of the levels. We use the $\rm{ICA}$ to calculate the scattering matrix $S$ whose elements’ squared $P_{ij}=|S_{ij}|^{2}$ represent the transition probabilities from $i$’th state at $-\infty$ to $j$’th state at $+\infty$.

Within the $\rm{ICA}$, one follows the diabatic levels, as illustrated in Fig.2(a) and (c), and apply the two-level L-Z formula at each crossing point along the path Sinitsyn (2015a, 2013, b); Landau and Lifshitz (1965); Zener and Fowler (1932). To apply the $\rm{ICA}$, the avoided level crossing points must be far from each other, and this is indeed satisfied in our case as shown in Fig. 2(b) and (d). During the quenching process when the system encounters the avoided level crossings it exhibits maximum probability of making transitions to the higher excited states. The necessary criteria for the applicability of the $\rm{ICA}$ is the cancellation of the dynamic phase (gained along the path of evolution) in the expression of the transition probabilities Sinitsyn (2015a, 2013, b). It indeed turns out that for the four-level problem the dynamic phases gained along all the possible trajectories corresponding to the transitions from the initial state $|0\rangle$ get canceled in the expression of the probability density, as shown in Appendix B. However, in the case of trajectories forward in time, all the transitions from states $|1\rangle$ and $|2\rangle$ exhibit the effects of the dynamical phases. On the other hand, for the six-level problem the dynamic phase difference remains finite for a couple of paths corresponding to the transition from the state $|\tilde{0}\rangle$ (the ground state at $-\infty$), as shown in the Appendix C.

It is worthwhile to point out that the above approximation is semi-classical because this does not allow an evolution backward in time. Therefore, transitions from $|0\rangle$ to state $|1\rangle$, and from $|3\rangle$ to state $|2\rangle$ are not possible in the case of four-level problem. However, for the six-level problem, transitions from state $|\tilde{0}\rangle$ to all other states are possible without requiring any evolution backward in time. In this sense, the six-level problem does not distinguish between the semi-classical and the quantum nature of the system. Further details of the calculations of the transition amplitudes, evaluated within $\rm{ICA}$, corresponding to both the four-level and six-level problems are described in Appendices B and C, respectively.

To our surprise, we observe that the analytical expressions corresponding to the probability $P^{(i)}_{00}$ (with $i=4,6$ corresponding to the four and six-level problems obtained in Appendices B and C, respectively) for the system to remain in the ground/initial state (of $t=-\infty$), within the semi-classical approximation, is very similar for both the four-level and the six-level problems. If we denote $p^{(4)}_{k}=P^{(4)}_{00}$ and $p^{(6)}_{k}=P^{(6)}_{00}$ as the transition probabilities within the ICA, we find that,

[TABLE]

where $p^{(4)}_{k}$ and $p^{(6)}_{k}$ differ by a factor of $1/2$ in the exponential. The defect density can be obtained straightforwardly by using the following formula,

[TABLE]

where the integration is performed over the HBZ as indicated above. It is worthwhile to point out that the probabilities $P_{0j}^{(i)}$ for the system to reach a state $|j\rangle$ at $t=\infty$ (except the corresponding ground state) also represent the defect densities. However, we have checked numerically for both the four-level and the six-level problems that the defect densities corresponding to $P_{0j}^{(i)}$ exhibit the same scaling with the quench rate $1/\tau$ as that of (19), thereby providing no new physics. Therefore, in what follows we evaluate (19) both analytically and numerically. From (18) and (19), we observe that the determining factor for the overall behavior of defect density $n^{(i)}_{d}$ is $|\Delta_{1\mathbf{k}}|^{2}$, which has the following expression,

[TABLE]

where $k_{\pm}=k_{x}+k_{y}\pm 2k_{z}$. For sufficiently slow quench $J\tau>>1$, the quantity $p_{k}$ is exponentially small for all values of $\mathbf{k}$ except on a contour dictated by $|\Delta_{1\mathbf{k}}|^{2}=0$ (see Fig. 3).

We henceforth refer to $|\Delta_{1\mathbf{k}}|^{2}=0$, for our case, as the non-adiabatic condition Sen et al. (2008); Mondal et al. (2008); Polkovnikov (2005); Polkovnikov and Gritsev (2008), and the corresponding points (satisfying the above condition) in the HBZ as the non-adiabatic points. Physically, it means that for all the non-adiabatic points and the points near these, the system exhibits the maximum non-adiabatic nature, thereby maximizing the probability of formation of the defects.

Next we calculate the defect density using (18) and (19) within the stationary phase approximation. We note that if a generic point on the above mentioned contour is given by $(k_{x}^{0},k_{y}^{0},k_{z}^{0})$, the most dominant contributions to the momentum integrals come from the values of $\mathbf{k}$ close to these points. Also, without loss of generality if we substitute $k_{x}=k_{y}\pm 2n\pi$ (corresponding to $\cos(k_{x}-k_{y})=1$ in (20) with $n$ being integers), we obtain from $|\Delta_{1\mathbf{k}}|^{2}=0$,

[TABLE]

where $\alpha=J_{x}/J_{y}$. The above equation indeed represents a set of critical lines (or, $d-m=1$ dimensional critical lines) for each values of $\alpha$ indicating $m=2$ in our 3D ($d=3$) model. However, the above conditions are valid only for all the values of $\alpha$ except $\alpha=0,\;\text{and}\,\pi/2$. For these two values of $\alpha$ the non-adiabatic condition is needed to be obtained directly from (20). It is therefore possible to Taylor-expand $|\Delta_{1\mathbf{k}}|^{2}$ around $(k_{x}^{0},k_{y}^{0},k_{z}^{0})$ on the above mentioned set of lines.

By making the above mentioned expansion up to leading order (i.e., the stationary phase approximation), we obtain

[TABLE]

where we have rewritten $\frac{\partial^{2}(|\Delta_{1\mathbf{k}}|^{2})}{\partial k_{i}\partial k_{j}}=4J^{2}f_{ij}(\alpha)$. The quantity $\left[\frac{\partial^{2}(|\Delta_{1\mathbf{k}}|^{2})}{\partial k_{i}\partial k_{j}}\right]$ is the $ij$’th element of the precision matrix $(2\times 2)$ of a multivariate Gaussian. For $J\tau>>1$, the values of $p_{k}$ corresponding to (18) is turn out to be very small for all values of $\mathbf{k}$ except the points near $(k_{x}^{0},k_{y}^{0},k_{z}^{0})$. We can then extend the upper limit of $k$-integration to $\infty$, and apply Gaussian integral in (19). In deriving the (22), we have used the condition $k_{x}=k_{y}\pm 2n\pi$, i.e., $k_{x}$ and $k_{y}$ are related to each other by a constant shift. This translates into $k_{x}-k_{x}^{0}=k_{y}-k_{y}^{0}$ which reduces the precision matrix to a $2\times 2$ from a $3\times 3$ one. Evaluating the Gaussian integral we find the defect density as,

[TABLE]

where $||f_{ij}(\alpha)||=4[f_{xx}f_{zz}-(f_{xz})^{2}]$ is the determinant of the $2\times 2$ precision matrix, which must be positive definite. It turns out that for $\alpha=0\;\text{and}\;\pi/2$ the above mentioned determinant vanishes and we can’t evaluate the Gaussian integral. Therefore, we have evaluated the corresponding integral in (19) numerically within the HBZ. Our analytical calculation using $\rm{ICA}$ thus establishes $n_{d}\propto(J\tau)^{-1}$, in the limit of the very slow quenching.

IV Defect density

In Fig. 4(a) and 4(b), we have plotted the defect density with respect to $J\tau$ and $\alpha$ for the four-level and the six-level L-Z problems respectively. The dots in Fig. 4 are the numerically (following the scheme outlined in II.3) obtained defect density. The 3D plot is obtained by using (19). We find that both the plots agree with each other with high degree of accuracy implying the validity of $\rm{ICA}$ in this model for both the four-level and the six-level problem. For slow quenching in both cases we obtain $n_{d}\sim\tau^{-1}$. The Appendix D features the ‘log-log’ plot (see Fig. 15 (a)) for $n_{d}$ vs. $J\tau$, which shows that for $\alpha=\pi/8,~{}\pi/4$ and $3\pi/8$ the scaling is indeed $\tau^{-1}$. It is worthwhile to point out that for these values of $\alpha$’s the system exhibits a 3D structure. For $\alpha=0,~{}\text{and}~{}\pi/2$ the system becomes a set of independent Kitaev chains, and Fig. 15 (b) shows that $n_{d}$ scales as $1/\sqrt{\tau}$ as expected.

Next we compare the above results with the prediction of Ref. Sen et al. (2008) that defect density crucially depends on the ‘dimensionality’ of the critical hyper-surface through which the system passes during the quenching. The important quantities are $z$, which tells us how the energy dispersion of a given system depends on momentum at low energy, the dynamical exponent $\nu$, which determines the asymptotic dependence of two point correlation function at large distance, and $m$, which tells us how the system becomes gapless in momentum space. In the case of our 3D Kitaev model, it can be checked easily that $\nu=z=1$. The dispersion of the 3D Kitaev model is such that it vanishes on a contour which constitutes a 2D critical hyper-surface under the action of quench yielding $m=2$. Keeping in mind that our system is a 3D one with $d=3$, and substituting the above values of $(z,\nu,d,m)$ in the expression of defect density $n_{d}\propto\tau^{-\frac{m\nu}{z\nu+1}}$, we obtain $n_{d}\propto\tau^{-1}$. It is worthwhile to point out that the 2D Kitaev model exhibits a scaling $n_{d}\propto 1/\sqrt{\tau}$ Sen et al. (2008).

We now discuss the gapless condition of the spectrum and its relation to the scaling of defect density. Diagonalization of (4) yields the following eigenvalues,

[TABLE]

where $\Delta_{k}=4|\Delta_{2k}|\sqrt{4J_{z}^{2}+|\Delta_{1k}|^{2}\cos^{2}\theta_{2k}}$. The gapless condition of the spectrum is found to be,

[TABLE]

We note from the above equations that the condition $|\Delta_{1k}|^{2}=0$, i.e., the non-adiabatic condition is a subset of the gapless condition. Thus our analysis suggests that for a general multilevel problem with partial or complete filling of the negative energy states the scaling of the defect density is decided by the coupling of the relevant levels in a given system.

However, the reason that the scaling of the defect density follows the law proposed in Ref. Polkovnikov (2005) is that for the case of the 3D Kitaev model the conditions for the zero of the spectrum, and the condition $|\Delta_{1k}|^{2}=0$ both vanish on a contour in an identical manner, (i.e., having $m=2$). We note that this is purely coincidental for the case of 3D Kitaev model.

Remarkably, the asymptotic behavior of defect density (i.e., the behavior of $n_{d}$ in the limit $J\tau>>1$) corresponding to the 2D Kitaev model differs from that of the 3D Kitaev model with respect to it’s dependence on $\alpha$. In the case of the 2D Kitaev model, when $\alpha$ is varied from 0 to $\pi/2$, the defect density increases monotonically up to $\alpha=\pi/4$ and then decreases monotonically such that it is symmetric with respect to $\alpha=\pi/4$ Sen et al. (2008). For the 3D Kitaev model although the asymptotic behavior of the defect density for large $J\tau$ is symmetric with respect to $\alpha=\pi/4$ however, in contrast to the 2D Kitaev model, it exhibits the maxima at $\alpha=0$ and $\pi/2$, and decreases monotonically and subsequently increases to reach local maximum at $\alpha=\pi/4$ as shown in Fig. 5 (a).

To explain this unusual feature we look back to the expression of $p_{k}$ appearing in (18) and (19).

Therefore, it is reasonable to define an effective density of states (effective-DOS) corresponding to the variable $\Delta_{1k}$, which essentially quantifies the number of $\mathbf{k}$-points in HBZ satisfying the non-adiabatic condition $\Delta_{1k}=0$.

From the effective-DOS plot in Fig. 5(b) we see that the effective DOS corresponding to $\Delta_{1k}=0$ is maximum at $\alpha=0~{}\text{and}~{}\pi/2$, and decreases initially away from these two values of the $\alpha$. However, the effective DOS exhibits a local maximum at $\alpha=\pi/4$. The infinite effective-DOS at $\alpha=0$ and $\pi/2$ refer to the contribution from each of the decoupled chains which are infinite in number. This explains the variation of asymptotic values of $n_{d}$ with respect to $\alpha$ in the sense that the more the number of $\mathbf{k}$-points in HBZ satisfy the non-adiabatic condition, the more the defect is generated.

V Correlation functions

In this section, we discuss the effect of the quenching in the correlation function and its connection to defect production. We note that the true quantum ground state of the 3D Kitaev model is supposed to be a quantum spin liquid where the spin-spin correlation (two spin correlation) function is short-ranged and anisotropic Baskaran et al. (2007). The spin-spin correlation is non-zero only for the nearest neighbor spins for all values of the parameters of the system. It remains an intriguing question how the correlation function develops non-locality due to the quenching. For simplicity we consider the following two-spin correlation function operator,

[TABLE]

which is nothing but a product of the spin operators corresponding to a site $\mathbf{R}$ at $a$ sub-lattice and a site $\mathbf{R+r}$ at $b$ sub-lattice. The quantity $\mathcal{O}^{3D}(\mathbf{r})$ is non-zero only for $\mathbf{r}=0$, and this fact remains the same for all points in the parameter space. In the Majorana fermion representation, the above correlation function takes the form $\mathcal{O}^{3D}(\mathbf{r})=S_{\mathbf{R},\mathbf{R+r}}ic_{1a}(\mathbf{R})c_{1b}(\mathbf{R+r})$ where $S_{\mathbf{R},\mathbf{R+r}}$ is the product of the $\mathbb{Z}_{2}$ gauge fields defined on the bonds along the path. In the following, we study the evolution of this two-fermion correlation function operator $ic_{1a}(\mathbf{R})c_{1b}(\mathbf{R+r})$ rather than the spin-spin correlation. Thus the object of our interest is,

[TABLE]

Doing a Fourier transform, the ground state expectation value of the two fermion correlation takes the following form,

[TABLE]

where $+$ and $-$ signs refer to the ground state (or initial state in the case of the four-level problem) corresponding to $J_{z}=-\infty$ and $\infty$ respectively and $\langle\mathcal{C}^{3D}(\mathbf{r})\rangle$ is indeed $\pm\delta_{\mathbf{r},0}$. The derivation of Eq. 28 is similar to previous quench study in 1D and 2D Kitaev models Mondal et al. (2008), and has been briefly described in Appendix E.

The final state after quench is superposition of both the initial ground state and excited states such that $p_{k}$ represents the probability of finding the system in the ground state of $t=-\infty$ (which is an excited state corresponding to $t=\infty$). In the limit $J_{z}\rightarrow-\infty$, the system is in the ground state where nearest-neighbor spins are ferromagnetically aligned and, as the quenching takes place, the state of the system deviates from the ground state configuration (corresponding to $t=-\infty$) and the defects are generated. Physically, $p_{k}$ quantifies how much ferromagnetic component is present in the final state given the fact that the initial state at $t=-\infty$ was a ferromagnetic state. In this case, the correlation function takes the form,

[TABLE]

where the second term in the above equation is a measure of defect correlation (i.e., the two-spin correlation in the presence of defects) and henceforth will be analyzed. We note that the final state characterizes the excitations or defects generated in the system due to the quench and, therefore, we refer the correlation function in the quenched state as the defect correlation following the earlier usages Mondal et al. (2008).

We evaluate the second term of the above equation by expanding the defect probability $p_{k}$, obtained within the $\rm{ICA}$, in the limit of very slow quench corresponding to $J\tau>>1$, for both four and six-level problems. As usual. the dominant contribution comes from the contour determined by $|\Delta_{1\mathbf{k}}|^{2}=0$. Evaluating the relevant Gaussian integrals, as has been done for the defect density in the previous section, we find,

[TABLE]

The above equation shows that the asymptotic behavior of the defect correlation as a function of $J\tau$ is given by,

[TABLE]

where $A=\frac{(N_{x}+N_{y})^{2}f_{zz}+4N_{z}^{2}f_{xx}-4(N_{x}+N_{y})N_{z}f_{xz}}{4\pi J||f_{ij}(\alpha)||}$, and $N_{x}=2n_{1}+n_{3}$, $N_{y}=2n_{2}+n_{3}$, and $N_{z}=2n_{3}$. The above asymptotic behavior, therefore, remains valid for both the four and the six-level problem.

V.1 Comparison of Defect correlation in various planes

After obtaining the asymptotic analytic expression of defect correlation function in (V), we analyze the latter in different directions and planes in the real-space. We evaluate (29) numerically in various planes and specific directions by noting that any arbitrary point on the lattice is given by $\mathbf{r}=n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2}+n_{3}\mathbf{a}_{3}$. First we begin by discussing how the spatial correlation builds in certain planes.

In Figs. 6, and 7 we plot the spatial variation of the defect correlation as a function of $n_{1}$, $n_{2}$, and $n_{3}$ for $J_{x}=J_{y}=J=1$ (i.e., $\alpha=\pi/4$), and also $J\tau=1$ for both the four-level and the six-level problems respectively. For the four-level problem $\langle\mathcal{C}^{3D}(\mathbf{r})\rangle$ has been evaluated by exact numerical calculation and for the six-level problem the same has been evaluated by both exact numerical as well using ICA for comparison. In Fig. 7, in the left panel we present the results obtained by exact numerical evaluation and in the right panel we present the same obtained within ICA. For a detailed description we refer to the caption of Fig. 7. (We have not used ICA for the four-level problem to avoid redundancy.)

Physically, the difference between the four-level and six-level problems concerning the correlation function study is that in the former the ground state of the system is such that only the $\psi$ fermion defined on the lower z-bond of the unit cell is occupied; see Fig. 8 for schematic illustration. On the other hand, for the six-level problem the ground state corresponds to both the $\psi$ and $\phi$ fermions defined on the both the z-bonds of the unit cell being occupied. We find that the correlation functions obtained in these two cases do not differ qualitatively because though in the four-level problem initially the $\phi$ fermion is absent, it affects the evolution of the $\psi$ fermion through the non-zero coupling $\Delta_{2k}$.

V.1.1 Correlation function in the $\mathbf{a}_{1}-\mathbf{a}_{2}$ plane

In Figs. 6(a) and 6(b), we plot defect correlation in the $\mathbf{a}_{1}-\mathbf{a}_{2}$ plane for $n_{3}=0$ and $n_{3}=1$ respectively for the four-level problem. Corresponding plot for the six-level problem are presented in Figs. 7(a) and 7(c) respectively. These figures indicate that the nearest neighbor correlations are the dominant ones and the next nearest and other long range correlations with far away sites quickly die out. However, occasionally the correlation with the far away sites build up due to the lattice geometry and the connectivity. It is worth noticing that the spatial variation of the defect correlation is symmetric in the $\mathbf{a}_{1}-\mathbf{a}_{2}$ plane with respect to the line $n_{1}=n_{2}$. Figures 7(b) and 7(d) represent the same correlation functions (plotted in Figs. 7(a) and 7(c)) obtained using ICA. We notice that ICA is able to show more oscillations than the exact numerical evaluation as the latter method has some practical limitations.

V.1.2 Correlation function in the $\mathbf{a}_{1}-\mathbf{a}_{3}$ plane

In Figs. 6(c) and 6(d), we plot the defect correlation in the $\mathbf{a}_{1}-\mathbf{a}_{3}$ plane for $n_{2}=0$ and $n_{2}=1$, respectively, for the four-level problem. The corresponding plots for the six-level problem are presented in Figs. 7(e) and 7(f), respectively. When comparing with the defect correlation in the $\mathbf{a}_{1}-\mathbf{a}_{2}$ plane, we find that in the $\mathbf{a}_{1}-\mathbf{a}_{3}$ plane the defect correlation is mostly zero and also not symmetric with respect to the line $n_{1}=n_{3}$ line.

V.1.3 Correlation function in the $\mathbf{a}_{2}-\mathbf{a}_{3}$ plane

Furthermore, owing to the symmetric nature of it in the $\mathbf{a}_{1}-\mathbf{a}_{2}$ plane the spatial variation of the defect correlation in the $\mathbf{a}_{1}-\mathbf{a}_{3}$ plane is the same as that of the $\mathbf{a}_{2}-\mathbf{a}_{3}$ plane. It is worthwhile to point out that the symmetric nature of the correlation function arises due to the condition $\alpha=1$ (or $J_{x}=J_{y}$). Therefore, we expect the defect correlation to develop spatial anisotropy for arbitrary values of the coupling constants corresponding to $J_{x}\neq J_{y}$.

V.2 Defect correlation along various directions

To investigate the above mentioned anisotropic nature of the correlation function in real space when the ratio $J_{y}/J_{x}$ is varied from [math] to $\infty$, while maintaining $J_{x}^{2}+J_{y}^{2}=1$, we plot in Fig. 9(a) and Fig. 9(b), the correlation $\langle\mathcal{C}^{3D}(\mathbf{r})\rangle$ as a function of $\alpha=\tan^{-1}\frac{J_{y}}{J_{x}}$ for the four-level and six-level L-Z quenches respectively. Here the quantity $\langle\mathcal{C}^{3D}(\mathbf{r})\rangle$ has been obtained by numerically evaluating the second term of (29) where $p_{k}$ has been obtained using both the $\rm{ICA}$ and the exact numerical calculations explained in II.3. All the correlation functions are the nearest-neighbor ones and obtained by substituting ${\bf{r}}=n_{1}\mathbf{a_{1}}+n_{2}\mathbf{a_{2}}+n_{3}\mathbf{a_{3}}$. Below we represent $\mathbf{r}$ by simply $(n_{1},n_{2},n_{3})$ and discuss defect correlation for some representative combinations of $(n_{1},n_{2},n_{3})$.

V.2.1 Correlation function for $(1,0,0)$.

First we discuss the behavior for ${\bf{r}}=(1,0,0)$ which describe how the correlation function behaves along $x$-direction. It is plotted in Fig. 9(a) and (b) in brown for the four-level and six-level problems, respectively. This represents a nearest neighbor correlation function joined by the $z-y$ bonds along the $x$-direction. For $\alpha=0$ (i.e., $J_{y}=0$) these two bonds are disconnected and we obtain a zero correlation. In this limit the lattice become a set of disconnected $x-z$ zig-zag chains as shown in Fig. 10(b). However, when $J_{y}$ becomes non-zero, finite correlation develops and reaches the maximum at $\alpha=\pi/4$. For $\alpha$ above $\pi/4$, it starts decreasing and asymptotically reaches to zero at $\alpha=\pi/2$ because this limit corresponds to $J_{x}=0$ for which the two bonds in the discussion again get disconnected and the lattice become a set of disconnected $y-z$ zig-zag chains as shown in Fig. 10(a). Also we observe that in the present context there is no qualitative difference between the four-level and six-level L-Z quench dynamics as expected.

V.2.2 Correlation function in $(0,0,1)$ direction

Next we calculate the correlation function for ${\bf r}=(0,0,1)$. The two $z$-bonds are next-neighbor along a $z-x$ zig-zag chain extending in vertical direction. For $\alpha=0$, it is expected to have maximum correlation and should decrease monotonically as we increase $J_{y}$. The red graph in Figs. 9 (a) and (b) show exactly the expected behavior. The corresponding lattice structure is shown in Fig. 10 (b).

V.2.3 Correlation function in $(1,0,1)$ and $(1,1,1)$ direction

For completeness, we present the correlation functions for ${\bf r}=(1,0,1)$ and ${\bf r}=(1,1,1)$. They are represented by blue and black graphs, respectively in Fig. 9. We notice that unlike ${\bf r}=(1,0,0)$, the maximum (or minimum) value does not appear for $\alpha=\pi/4$, though they have zero value for $\alpha=0$ and $\pi/2$. This is attributed to the fact that to connect two $z$-bonds of type-1 in $(1,0,1)$ and $(1,1,1)$, one needs to traverse an unequal number of $x$ and $y$ bonds which results in shifting their respective maxima away from $\alpha=\pi/4$

We re-emphasize that the zero value of correlation function in all of the above-mentioned cases can be easily understood by noticing that for $\alpha=\pi/2$, the lattice reduces to disconnected $z-y$ chains as drawn in the right panel of Fig. 10 (a). Similarly, at $\alpha=0$ the lattice becomes disconnected $z-x$ chains as depicted in Fig. 10 (b).

VI Entropy

In this section, we calculate the von-Neumann entropy of the final quenched state. We note that the initial state is a product of the pure-states and remains so due to the unitary evolution. However, for a given momentum ‘$k$’, the quenched state at any instant is a linear superposition of the instantaneous ground state and the other excited states whose probabilities exhibit explicit time dependence. It causes the off-diagonal matrix elements of the density matrix to be time-dependent (and oscillatory too) and leads to intrinsic decoherence such that for long quench times, the final quenched state is equivalent to a decohered state Cherng and Levitov (2006). It is, therefore, pertinent to calculate the von-Neumann entropy of this quenched states at $t\rightarrow\infty$. The physical reason behind this entropy generation is the underlying non-equilibrium processes in the quench dynamics and subsequent production of defects.

As explained above, for the four-level L-Z problem, the instantaneous state of the system $|\Psi_{4}(k,t)\rangle$ can be expressed as $|\Psi_{4}(k,t)\rangle=\sum_{i}\sqrt{A_{i}(k,t)}|i(k,t)\rangle,~{}~{}i=0,1,2,3$, where $|i(k,t)\rangle$ represents the instantaneous diabatic states(see Fig. 2) and $A_{i}(k,t)$ are the respective probabilities.

Keeping in mind that the off-diagonal element of the density matrix constructed from $|\Psi_{4}(k,t)\rangle$ does not contribute to the expectation value of any thermodynamic quantity Cherng and Levitov (2006), the von-Neumann entropy of the quenched state at $t=+\infty$ is given by,

[TABLE]

where the above equation contains only the time independent diagonal elements of the density matrix. The above equation is evaluated numerically to obtain the von-Neumann entropy of the final decohered state. Moreover, one can write an equivalent expression of $|\Psi_{4}(k,t)\rangle$ within ICA as,

[TABLE]

where the eigenvalues $E_{n,k}$ are essentially corresponding to $t=\pm\infty$ and $P_{ij}$’s are the matrix elements of the transition matrix (B). However, the final density matrix of the system still remains diagonal as the rapidly oscillating off-diagonal terms do not affect any physical quantity  Cherng and Levitov (2006); Mukherjee et al. (2007). Therefore, the von-Neumann entropy can be evaluated in the quenched state by the following formula,

[TABLE]

In the limit of very small $\tau$, i.e., $\tau\rightarrow 0$, $P_{00}\rightarrow 1$ and $P_{02}=P_{03}=0$, the von-Neumann entropy is zero. A procedure similar to (32) can be implemented for the six-level L-Z problem, and the von-Neumann entropy can be computed numerically. However, contrary to the four-level problem, an expression of the quenched state equivalent to (33) can not be obtained for the six-level L-Z problem due to the appearance of the dynamic phases in some of the probability amplitudes (see Appendix C).

The von-Neumann entropy of the initial state is zero because in this case $P_{00}=1$ (see (34)). The non-equilibrium process caused by quenching is expected to generate finite entropy by inducing transitions to the other states such that $P_{00}<1$. The von-Neumann entropy is plotted as a function of $J\tau$ and $\alpha$ in Figs. 11(a) and 11(b) for the four-level and six-level problems, respectively. We notice from Fig. 11 that for a given $\alpha$, the entropy reaches the maximum at a certain value of $J\tau$, and then decreases monotonically and asymptotically saturates to a finite value. This finite value appearing at $t\rightarrow+\infty$ is due to the fact that the final state is a linear superposition of the ground state of $t=+\infty$, and the defect states.

To obtain a more physical insight to this final decohered state originating from the quenching, we plot the correlation function as a function of $J\tau$ for some representative values of $\alpha$ in Fig. 12. The correlation function exhibits the form, $\frac{1}{\tau}e^{-A\tau}$ derived in (31), as a function of quench time $J\tau$. Comparing Figs. 11 and 12, we see that for a given $\alpha$, the entropy behaves the same as the correlation function as a function of $J\tau$. It indicates that the quenching generates, through the defect production, the correlation between the Majorana fermions (and therefore, spins) beyond the nearest neighbors, which manifests itself in the increased entropy of the final state. Also, the similar dependence of the entropy (of the final state) and the correlation function as a function of $J\tau$ conform to the fact that the quenched state is an inherently decohered state whose density matrix can be represented by the time independent diagonal elements. These diagonal elements, in effect, determine the expectation values of any thermodynamic quantities (correlation function in this case)Cherng and Levitov (2006).

VII Conclusions and discussions

To summarize, we have studied the quench dynamics and defect production in the spin- 1/2 3D Kitaev model within the set-up of a multi-level L-Z problem. By mapping the spin Hamiltonian to a quadratic Hamiltonian of Majorana fermions hopping in the background of a uniform static $\mathbb{Z}_{2}$ gauge field corresponding to each bond, we have obtain a Hamiltonian with $4\times 4$ block diagonal form in the momentum space. We have constructed a four-level and a six-level L-Z problem by filling up the lowest negative energy state and both the negative energy states respectively. Our analytical calculations using the ICA agree with the corresponding exact numerical calculations. Our major findings are the following:

In the limit of very slow quench, viz., $\tau>>1$ the defect density $n_{d}$ scales as $\tau^{-1}$ with the quench time $\tau$ and satisfies the general scaling law $n_{d}=\tau^{-\frac{m\nu}{z\nu+1}}$, with $m=2$, and $\nu=z=1$, where $\nu$ is the critical exponent of the correlation length, and $z$ is the dynamic critical exponent. Moreover, in the $\tau>>1$ limit for a fixed $\tau$ the defect density is found to be maximum at $\alpha=0$ and $\pi/2$, and at $\pi/4$ exhibits a local maximum.

• 2.

In the limit of very slow quench the defect correlation function is found to scale as $\tau^{-1}e^{-A/\tau}$ with the quench time $\tau$. Furthermore, the correlation function is found to be spatially symmetric in the $x-y$-plane (or, equivalently, in the $\mathbf{a}_{1}-\mathbf{a}_{2}$ plane), and asymmetric in the other two planes for $J_{x}=J_{y}$. However, in general, the defect correlation is an-isotropic for $J_{x}\neq J_{y}$.

• 3.

We have evaluated the residual entropy in the final state after quenching as a function of the system-parameter. The dependence of the entropy on the quench rate is qualitatively the same as that of the correlation function.

In the following, we discuss all the above mentioned findings one by one.

The scaling behavior of the defect density obtained in the 3D Kitaev model considered hare is a direct consequence of the way the system enters into the non-adiabatic regime. The non-adiabatic condition, governed by $|\Delta_{1k}|^{2}=0$, where $\Delta_{1k}$ is the relevant coupling between the levels determining the defect density, does not ensure the spectrum of the system to be gapless. The fact that the defect density satisfies the general scaling law $n_{d}=\tau^{-\frac{m\nu}{z\nu+1}}$, proposed in Ref. Polkovnikov (2005); Mondal et al. (2008), indicates a deeper connection between the gapless condition for the energy spectrum and the non-adiabatic condition. For the present model, the parameter $\Delta_{1k}$ appears within the condition for the spectrum to be gapless, but the zeros of $\Delta_{1k}$ do not ensure the gaplessness of the spectrum. Incidentally, the zeros of $\Delta_{1k}$ and the spectrum both constitute their own set of $d-m=1$ dimensional critical lines which ensure the general scaling law being satisfied. This is unlike the 2D Kitaev model where the gapless condition itself determines the non-adiabaticity Mondal et al. (2008). We further speculate that the departure of the scaling law from the above mentioned behavior may be possible in a multi-level L-Z problem, if the relevant non-adiabatic condition exhibits a different dimensional dependence than the condition for the energy spectrum to be gapless. Finding such a possibility or lack of it, we believe, is a definitive direction of theoretical investigation in the future.

The dependence of the defect density on $\alpha$ originates form the nature of the effective-DOS corresponding to the function $\Delta_{1k}$ in the gapless phase as shown in Fig. 5. Physically the effective-DOS quantifies the number of $\mathbf{k}$-states satisfying the non-adiabatic condition. When the quenching process takes the system across the gapless phase the number of $\mathbf{k}$-states satisfying the non-adiabatic condition is the maximum at $\alpha=0$ and $\pi/2$. This leads to a larger defect formation at these values of $\alpha$. On the other hand, the number of $\mathbf{k}$-states satisfying the non-adiabatic condition at $\alpha=\pi/4$ is locally maximum, making the defect density locally maximum as well.

Moreover, when $J_{x}=J_{y}$ the lattice retains its full 3D structure and the nearest-neighbor dimers are placed along the $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ directions. From Fig. 1 it is easy to recognize that when $J_{x}=J_{y}$ a dimer finds an equal number of nearest-neighbor dimers in both $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ directions. This makes the defect correlation in the $\mathbf{a}_{1}-\mathbf{a}_{2}$ plane to become symmetric with respect to the $|\mathbf{a}_{1}|=|\mathbf{a}_{2}|$ line. For other planes such a condition is not satisfied. It turns out that when either $J_{x}$ or $J_{y}$ vanishes, the lattice becomes a set of decoupled chains (see Fig. 10) and only those correlations survive for which a dimer can find another dimer in the same chain.

Starting from a trivial product state of ferromagnetic dimers the quenching process takes the system to a final state whose density matrix can be approximated as a suitable reduced density matrix exhibiting a finite value of the von-Neumann entropy. Most certainly, the final state contains the ground state of the Hamiltonian at $t=\infty$, and the other excited states. The post quench (at $t=\infty$) ground state of the Hamiltonian exhibits a structure where anti-ferromagnetic dimers are sitting on each $z$-links. In this regard, it is worthwhile to point out that our definition of the defect density quantifies the ferromagnetic component of the final state. The maxima of the correlation functions and the entropy are reached around same values of $J\tau$. This indicates that the quenching process gives rise to the generation of a non-local correlation which is further manifested in the finite values entropy.

There exists a class of 3D materials, such as $\beta-$ and $\gamma-\text{Li}_{2}\text{IrO}_{3}$, containing Kitaev interactions as the effective low energy description Winter et al. (2017). These materials offer a possibility to experimentally realize the quench dynamics in the 3D Kitaev model studied here. However, in these materials apart from the Kitaev interactions other interactions such as, Heisenberg interaction, asymmetric interaction of the type $S^{x}S^{y}$, or even Dzyaloshinskii-Moriya interaction are present which destroy the exact solvability of the system. We expect that the external pressure or magnetic field may become useful in realizing a similar quench dynamics given the fact that these parameters have been useful to realize the much-coveted spin-liquid phase in the Kitaev materials Takayama et al. (2015); Breznay et al. (2017); Veiga et al. (2017). In particular, it has been pointed out that $x-x,~{}~{}y-y$ or $z-z$ type interactions can be changed by applying pressure Veiga et al. (2017). Therefore, one can, in principle, achieve the quench protocol we have considered, viz., the variation of the strength of $z-z$ interaction. We also believe that an alternative way to realize our quench study, as well as quench studies in other multi-level L-Z problem, would be in the optical lattices as envisaged in a recent experiment Keesling et al. (2019).

Acknowledgements

S.M. acknowledges the earlier collaboration with Anirban Dutta with whom the project actually started but could not be continued. S.S. acknowledges the partial financial support provided through a Kreitman postdoctoral fellowship and the Israel Science Foundation (ISF) Grant No. 1360/17, and Sayak Ray for critically reading the manuscript.

Appendix A Ground states of the single and the two-particle sector

We have performed a unitary transformation to diagonalize the upper-left and lower-right $2\times 2$ blocks of (3). As a result, we get the following relationship between the old $\Psi_{k}$ fermions (defined in (3)) and new $\eta_{i,k}$ fermions:

[TABLE]

where $u_{k}$ is the $2\times 2$ unitary matrix used in the above-mentioned diagonalization, and ‘$k$’ has been defined only over the HBZ. The explicit expression for $u_{k}$ is given by,

[TABLE]

The initial state for the four-level problem, which we denote by $|\mathcal{I}\rangle$, is obtained by filling only the $\eta_{0k}$ of the vacuum of $\eta_{ik}$ ($|\tilde{\emptyset}\rangle$). Similarly, the ground state $|\mathcal{G}\rangle$ for the six-level problem is obtained by filling $\eta_{0k}$ and $\eta_{1k}$. We then find the following relations,

[TABLE]

In the above, $|\emptyset\rangle$ is the original vacuum of the fermions representing a zero-fermion state, and we used the fact that $|\emptyset\rangle$ is already a vacuum for $\eta_{2k}$ and $\eta_{3k}$.

All four singly-occupied states corresponding to the four-level problem are given by, $|0\rangle=\eta_{0k}^{\dagger}|\tilde{\emptyset}\rangle$, $|1\rangle=\eta_{1k}^{\dagger}|\tilde{\emptyset}\rangle$ , $|2\rangle=\eta_{2k}^{\dagger}|\tilde{\emptyset}\rangle$, and $|3\rangle=\eta_{3k}^{\dagger}|\tilde{\emptyset}\rangle$. With these notations, we have defined $\Phi_{4,k}$ in (9). Similarly, for the six-level problem, we enumerate all the six doubly occupied states in terms of the single particle states and their respective energies as follows:

[TABLE]

Within this notation, we have defined, in the Schrodinger equation governing the time evolution of the system, $\Phi^{\dagger}_{6,k}=(|\tilde{0}\rangle,|\tilde{1}\rangle,|\tilde{2}\rangle,|\tilde{3}\rangle,|\tilde{4}\rangle,|\tilde{5}\rangle)$.

Appendix B Transition amplitudes: four-level problem

To calculate the transition amplitudes, we first consider all possible semi-classical paths corresponding to Fig. 2(a). As mentioned in section II, the ground state at $t=-\infty$ is $|0\rangle$. First we note that there are four crossing points, viz.,

[TABLE]

We apply the LZ formula by using the above coupling constants corresponding to the respective crossing points. We note that whenever a state crosses another state with coupling constant $g_{k}$ (as well as $-g_{k}$), the probability amplitude that the system remains in the same state is given by $u_{g}=\exp\left(-\frac{\pi|g_{k}|^{2}(J\tau)}{4J^{2}}\right)$. Therefore, the probability amplitude that the system makes a transition is $i\xi\sqrt{1-(u_{g})^{2}}$, where $\xi=1$ when the sign of the coupling constant is positive and $\xi=-1$, when negative. These expressions for the crossing points with coupling $\gamma_{k}$ are same except $|g_{k}|$ shall be replaced by $|\gamma_{k}|$. In the following, let us write down the transition amplitudes $S_{ij}$ for all the semi-classically allowed paths. In writing such paths we need to consider the effect of the dynamic phase $e^{i\phi_{d}^{ij}(\nu)}$ in a given path $\nu$ which is given by,

[TABLE]

where $\epsilon_{\nu k}(t)$ are the diabatic energy levels participating in the trajectory. This dynamic phase is nothing but the area under the trajectory $\nu$ under consideration. All the transition amplitudes corresponding to the transitions from the state $|0\rangle$ to the rest of the states are given by,

[TABLE]

In the above equation the dynamic phases in the respective trajectories get canceled when the square of the modulus of the transition amplitudes are taken. Therefore we have dropped the term $e^{i\phi_{d}^{0j}(\nu)}$ corresponding to the dynamic phase. Similarly, all the transition amplitudes corresponding to the transitions from $|1\rangle$ to the rest of the states can be written down. Here, while considering all the semi-classical trajectories, we need to consider the effects of the dynamic phases on the transition amplitude in the corresponding trajectories,

[TABLE]

where in the above equation $e^{i\phi_{d}^{ij}(1)}$ and $e^{i\phi_{d}^{ij}(2)}$ are the dynamic phases corresponding to the trajectories (1) and (2), respectively, and the phase difference is given by, $\Delta\phi_{d}(2,1)=\phi_{d}^{ij}(2)-\phi_{d}^{ij}(1)=4J_{z}T^{2}+|\Delta_{2k}|T$, which is the area in the $\epsilon-t$ plane bounded by both the trajectories. $\pm T$ are the points on the $t$-axis where the diabatic levels cross each other and the time axis too. It is easy to recognize from Fig. 13 that the phase difference is the same for all the four transition amplitudes where the dynamic phase plays its role. Likewise, all the transition amplitudes corresponding to the transitions from $|2\rangle$ are given by,

[TABLE]

Lastly, all transition amplitudes corresponding to transitions from the state $|3\rangle$ are given by,

[TABLE]

Therefore, the transition probability matrix is given by,

[TABLE]

where the matrix elements of the above matrix is given by, $P_{ij}=S_{ij}^{*}S_{ij}$, and

[TABLE]

The quantity $P^{(4)}_{00}$ used in (18) is given by $P^{(4)}_{00}=(u_{g}u_{\gamma})^{2}$.

Appendix C Transition amplitudes: six-level problem

The six-level matrix representation for the 3D Kitaev Model is given by,

[TABLE]

with $J_{3}=Jt/\tau$, where the degeneracy removal term $\pm\epsilon$ has been used for the horizontal levels at zero energy. This is perhaps the only way to apply ICA in a system exhibiting more than one degenerate horizontal energy levels passing through zero energy. However, to take the limit $\epsilon\rightarrow 0$ is mandatory at the end of the calculation, and in this limit $T_{1}$ also vanishes. The couplings at the crossings are given by,

[TABLE]

It is worthwhile to note that there is no direct coupling between the states $|0\rangle$ and $|5\rangle$. At each crossing point, the crossing involves a horizontal level (with zero slope) and a state with definite non-zero slope $\pm 4J/\tau$. Therefore, the denominator in the corresponding two-level L-Z transition amplitude is always $4J/\tau$ for all the available independent two-level L-Z transitions. There are only two types of the transition amplitudes for diabatic passage through any crossing point, viz., $u_{g}=\exp{\left(-\frac{\pi|g|^{2}J\tau}{4J^{2}}\right)}$, and $u_{\gamma}=\exp{\left(-\frac{\pi|\gamma|^{2}J\tau}{4J^{2}}\right)}$ whereas, the transition amplitude for the adiabatic passage through the crossing points are $i\xi\sqrt{1-u_{g}^{2}}$, and $i\xi\sqrt{1-u_{\gamma}^{2}}$, with $\xi=\pm 1$ as mentioned earlier.

Let us start by finding all the amplitudes corresponding to the transitions from state $|0\rangle$ to all other states, viz.,

[TABLE]

where $\phi_{0i}$’s are the dynamic phases associated with the paths. In the above,

[TABLE]

In the above equations, the dynamic phases associated with the paths are given by,

[TABLE]

[TABLE]

It is worthwhile to note that once again (as was the case for four-state model) the dynamic phase difference is equal to the area under the curve. The transition probabilities are given by, $P_{00}=|S_{00}|^{2}=(u_{g}u_{\gamma})^{4}$, this is the quantity we have defined as defect probability in the equation (18). Similarly, $P_{01}=|S_{01}|^{2}=P_{04}=|S_{04}|^{2}=(u_{g}u_{\gamma})^{2}(1-u_{\gamma}^{2})$, (note that they are equally probable) $P_{03}=|S_{03}|^{2}=(u_{\gamma})^{4}u_{g}^{2}(1-u_{g}^{2})$, and

[TABLE]

Although we have defined $P_{00}$ only as the defect probability, the other probabilities except $P_{05}$ are still representing probabilities of the system going into different adiabatic states other than the adiabatic ground state $|5\rangle$ at $t=\infty$. As a whole, the total probability must be one, i.e., $(P_{00}+P_{01}+P_{02}+P_{03}+P_{04}+P_{05})=1$. Furthermore, it is easy to determine the relation between the time $T_{2}$ and the other parameters of the system to be $T_{2}=\frac{|\Delta_{2k}|}{2J^{2}}(J\tau)$, from the fact that the slope of the diabatic energy curve, $\frac{2|\Delta_{2k}|}{T_{2}}=\frac{4J}{\tau}$. In the limit of $J\tau>>1$, we expect the dynamic phase difference to be large enough and can be considered as a random phase.

In the above, analysis we have considered all the transitions from the state $|\tilde{0}\rangle$ since it is the ground state of our 3D-Kitaev model at $t=-\infty$. One can follow the methodology explained above to get the entire transition probability matrix as has been done for 4-state representation of the model however, this is unnecessary for our purpose.

Appendix D Fitting of numerically obtained defect density

In the following we have plotted the $\ln(n_{d})$ corresponding to the six-level representation as a function of $\ln(J\tau)$ for different values of $\alpha$. We have found that when $\alpha=0$ and $\pi/2$ the exponent of the defect density is ‘$1/2$’ which is as expected since in these two limits the 3D Kitaev model become disconnected Kitaev chains exhibiting a scaling of $\tau^{-1/2}$. For other values of $\alpha$ the defect density scales as $\tau^{-1}$ according to ICA and as shown in the Fig. 15 (a) the numerically obtained scaling is $\tau^{-0.9}$ which matches with the ICA value (within the numerical accuracy). We have checked that the four-level problem exhibit exactly the same quantitative behavior with respect to the scaling of the defect density. However, we have not provided the corresponding plots to avoid redundancy.

Appendix E Derivation of the correlation function expression corresponding to the equation (28)

The correlation function as defined in (27) can be expressed in terms of the original $\psi(\mathbf{R})$ and $\phi(\mathbf{R})$ fermions in real space as follows,

[TABLE]

Now we use the Fourier transformation $\psi(\mathbf{R})=\frac{1}{\sqrt{N}}\int d^{k}e^{i\mathbf{k}\cdot\mathbf{R}}\psi_{k}$. We notice that contributions from the 1st and the last term are zero, and contributions from the second and the third are conjugate to each other. With the above expressions we use the definitions in (35) to re-express correlation function in terms of $\alpha_{ik}$ and using orthogonality of states we arrive at (28).