Dynamics of charge imbalance resolved negativity after a global quench in free scalar field theory

In this paper, we consider the time evolution of charge imbalance resolved negativity after a global quench in the 1+1 dimensional complex Klein-Gordon theory. We focus on two types of global quenches which are called boundary state quench and mass quench respectively. We first study the boundary state quench where the post-quench dynamic is governed by a massless Hamiltonian. In this case, the temporal evolution of charged imbalance resolved negativity can be obtained first by evaluating the correlators of the fluxed twist field in the upper half plane and then applying Fourier transformation. We numerical test our analytical formulas in the underlying lattice model. We also study the mass quench in the complex harmonic chain where the system is evolved according to a massive Hamiltonian after the quench. We argue that our results can be understood in the framework of quasi-particle picture.


Introduction
In the past two decades, the investigation of various entanglement measures has greatly sharpened our understanding of quantum many-body system, quantum field theory and quantum gravity. In condensed matter theory, entanglement is a powerful tool to characterize different phases of matter [1][2][3]. In the AdS/CFT correspondence, the Ryu-Takayanagi formula [4,5] firstly opens the route of understanding spacetime from entanglement and this idea turns out to have a key role in the black hole information loss paradox [6][7][8]. Entanglement is also an important concept in the studies of equilibration and thermalization of isolated quantum systems. Among all these progress, entanglement entropy is the most successful entanglement measure to characterize the bipartite entanglement of a pure state. When the system is prepared in a pure state |ψ , the reduced density matrix (RDM) of subsystem A is defined by tracing out its complement B, ρ A = tr B |ψ ψ|. From the moments of ρ A , i.e. Trρ n A , one can obtain the Von Neumann entropy through the replica trick [9] S ≡ −Tr(ρ A log ρ A ) = lim n→1 S n (1.1) where S n is the Rényi entropies S n = 1 1 − n log Trρ n A . (1.2) Now suppose we are interested in the entanglement between two subsystems A 1 and A 2 , which are not necessarily complementary to each other. In this situation, ρ A 1 ∪A 2 is general a mixed state, and Von Neumann entropy is no longer a good measure of entanglement. Among different proposals, a computable measure of mixed state entanglement, entanglement negativity (or logarithmic negativity equivalently) turns out to be very useful [10][11][12]. The definition is where ||O|| = Tr √ O † O denotes the trace norm of the operator O and ρ T 2 A is the partial transpose of RDM ρ A with respect to degree of freedom of subsystem A 2 . Let |e (1) i and |e (2) j be two arbitrary bases of the Hilbert spaces associated to the degree of freedom on A 1 and A 2 respectively. The partial transpose (with respect to the second space) of ρ A is defined as i e (2) l | ρ A |e (1) k e (2) j .
In this paper, we will consider the time evolution of the charge imbalance resolved negativity in a special kind of non-equilibrium state known as global quantum quenches. A global quantum quench describes a process in which the sudden change of the Hamiltonian H 0 → H at a given time that we set as t = 0, then the initial ground state |ψ 0 of the pre-quench Hamiltonian H 0 evolves according to post-quench Hamiltonian H. Thus, we have |ψ(t) = e −iHt |ψ 0 . The quench dynamics of various measures of entanglement have been extensively studied in the literature [45][46][47][48][49][50][51][52][53][54][55]. In this paper, we will focus on two types of global quenches in 1+1 dimensional complex Klein-Gordon field theory. The first type is the boundary state quench in which after the quench the mass of the scalar field is zero, so the time evolution is governed by a conformal field theory (CFT). The second type we will consider is a mass quench in the underlying lattice model (complex harmonic chain), where the post-quench dynamic is governed by a massive Hamiltonian.
The remaining part of this manuscript is organized as follows. In section 2, we briefly review some basic facts about charge imbalance resolution of entanglement negativity. In section 3, we discuss how to calculate the dynamics of charged Rényi negativity in the boundary state quench setup using CFT techniques. In section 4, we will consider the quench dynamics of charged (logarithmic) Rényi negativity between two finite intervals in the boundary state quench protocol. In section 5, by applying Fourier transformation, we can obtain the charge imbalance resolved negativity from the results derived in the previous section. In section 6, we check our analytical predictions against numerical computation in the complex harmonic chain. In section 7, we further discuss the mass quench in the underlying lattice model numerically and using quasi-particle interpretation to predict the quench dynamics of charged logarithmic negativity, finding perfect agreement. Finally, we conclude in section 8 and discuss some possible interesting extensions of this work.

Charge imbalance resolution of negativity
In this section, we will briefly review the decomposition of entanglement negativity under a global internal symmetry.
Let us first review some basic facts about the symmetry resolution of the entanglement entropy. We assume that our system exhibit a global U (1) symmetry generated by a local charge Q. If [ρ, Q] = 0 (this can be achieved if ρ only acts non-trivially on the eigenspace of Q), then we have [ρ A , Q A ] = 0, which implies that ρ A admits charge decomposition according to eigenvalues q of local charge Q A where P q is the projection operator which projects the space to the eigenspace corresponding to eigenvalue q. The symmetry resolved Rényi entropies are defined as It's convenient to first introduce the charged moments of ρ A , Then it's sufficient to compute its Fourier transform to obtain the Rényi entropies of the sector with charge q as Finally, the symmetry resolved entanglement entropy can be obtained by taking the replica limit S(q) = lim n→1 S n (q). Now we briefly review the symmetry decomposition of entanglement negativity under the U (1) charge Q. Since the charge Q is local, we can write Q A = Q 1 + Q 2 , where Q 1 and Q 2 are the charges corresponding to sub subsystem A 1 and A 2 respectively. From the relation [ρ A , Q A ] = 0, performing a partial transposition with respect to the second region A 2 of subsystem A, we obtain where we have introduced the charge imbalance operator Q A and we will denote its eigenvalues as q to make a distinction with the eigenvalues of Q A . Then ρ T 2 A has a block matrix form, each block was characterized by different eigenvalues q of the imbalance operator Q A . If we write . (2.7) Then we have where p(q) = Tr(P q ρ T 2 A ) is the probability of finding q as the outcome of measurement of Q A . The charge imbalance resolved negativity is defined as The total negativity is given by the sum of charge imbalance resolved negativity weighted by the corresponding probability N = q p(q)N (q). (2.10) It's useful to define the charge imbalance resolved Rényi negativity Then the charge imbalance entanglement negativity can be obtained by taking the limit The projection operator P q has the following integral representation It's convenient to first introduce the charged Rényi negativity and it's Fourier transformation then the charge imbalance resolved Rényi negativity are related by Fourier transform It's also useful to introduce the charged Rényi logarithmic negativity E n (µ), defined as The charged logarithmic negativity is defined by taking the following replica limit Then the charge imbalance resolved negativity can be obtained as In section 5, we will use eq. (2.19) and eq. (2.16) to evaluate the charge imbalance resolved negativity.

Boundary state quench
In this paper, we will consider the 1+1 dimensional complex free scalar field theory with the Euclidean action given by This action exhibit a U (1) symmetry, i.e. the phase transformation of the field φ → e iθ φ, φ † → e −iθ φ † leaves the action invariant. The Hamiltonian of this theory is with π is the canonical momentum for φ.
We can rewrite this theory in terms of two real scalar fields φ (1) and φ (2) with φ = 1 √ 2 (φ (1) + iφ (2) ). The Hamiltonian in terms of these variables becomes The Hamiltonian can be diagonalized in terms of particle and anti-particle mode operators a(p), a † (p) and b(p), b † (p) with the commutation relation [a(p), a † (q)] = [b(p), b † (q)] = 2πδ(p − q) and all other commutators vanishing. We have with e(p) = m 2 + p 2 , while the conserved charge corresponding to the global U (1) symmetry is In the second equal sign of above equation, the conserved charge are expressed as integral of local density in real space. Thus its value in a given subsystem A is the same integral restricted to A, The lattice version of the complex Klein-Gordon field theory is the complex harmonic chain which is equivalent to two decoupled real harmonic chains. The Hamiltonian of the real harmonic chain made by L sites reads where periodic boundary conditions φ L ≡ φ 0 , π L ≡ π 0 are imposed and variables π j and φ j satisfy The lattice version of the complex scalar field theory is the sum of two of the above harmonic chain. In terms of the variables φ (1) , π (1) and φ (2) , π (2) , the Hamiltonian is The Hamiltonian eq. (3.8) can be diagonalized by introducing the creation and annihilation oper- In terms of these operators, the Hamiltonian eq. (3.8) is diagonal The conserved charge is local and can also be written in the position space and for a given subsystem A reads

The path integral approach to boundary state quench
In this section, we will focus on the case in which after the quench the mass of the scalar field is zero, so time evolution is governed by a CFT. Let's first briefly review the imaginary time formalism for the global quenches where the post-quench dynamic is govern by a CFT. This quench setup is called a boundary state quench as the correlation measures after quench can be described by boundary conformal field theory (BCFT), where the gapped initial state are represent by a spacetime boundary.
The expectation value of a product of equal-time local operators in time dependent state |ψ(t) is where two factors e −τ 0 H have been introduced to make the path integral representation of this expectation value absolutely convergent and Z = ψ 0 | e −2τ 0 H |ψ 0 is the normalization factor. The above prescription is equivalent to assuming that the initial state has the form |ψ 0 ∝ e −τ 0 H |B , where |B is some conformal boundary state. One could interpret τ 0 as being proportional to the correlation length of the initial state. Thus, the predictions made by this approach are expected to be valid only in the spacetime scaling limit, The above correlator admits a path integral represent where L is the Euclidean Lagrangian corresponding to the post-quench Hamiltonian H. We need to identify τ 1 = 0 and τ 2 = 2τ 0 . We should compute the path integral considering τ real and only at the end of the computation to analytically continue it to τ = τ 0 + it.

Evolution of charged moments of RDM
In this section, we briefly review the strategy of computing the charged Rényi negativity in twodimensional CFT (see [40] for more details). In a generic two-dimensional QFT, we can view the charged moments Z n (µ) (defined in eq. (2.3)) as the partition function on the Riemann surface R n,N pierced by an Aharonov-Bohm flux, such that the total phase accumulated by the field upon going through the entire surface is µ. The presence of the flux corresponds to impose an additional twist on the boundary of subsystem A. This additional twist fuses with the replica (ordinary) twist field at the endpoints of subsystem A and can be implemented by two local fields T n,µ andT n,µ named as fluxed twist fields and fluxed anti-twist fields. These fluxed twist fields take into account not only the internal permutational symmetry among the replicas but also the presence of the flux. The partition function on the fluxed Riemann surface is thus proportional to the 2N -point function of these fluxed twist operators.
We consider the subsystem A consists of two disjoint intervals on the real axis, The charged moments of RDM of vacuum state ρ A , i.e. trρ n A e iµQ A is equivalent to a four-point function of the fluxed twist fields where w = u + iτ with τ ∈ R, τ ∈ (0, 2π) is the coordinate of the strip. The fluxed twist fields T n,µ and fluxed anti-twist fieldT n,µ are primary operators and they have the same dimension [23,27] The charged moments of the partially transposed RDM can be obtained from the correlator above with the fluxed twist field T n,µ andT n,µ at the endpoints of A 2 exchanged while the remaining ones keep the same, giving By conformal map z(w) = e πw/2τ 0 , the four point function on the strip are mapped to the upper half plane (UHP), then we can write By global conformal symmetry, four-point functions of primary fields on the upper half plane depend on 6 cross ratios η i,j and have the following structure depend on the full operator content of the theory and usually is very complicated. The constants c n,µ are also non-universal and are known for some specific theories.
While for the charged Rényi negativity, we have Fortunately, in the spacetime scaling regime, we have η i,j → 0, 1 or ∞, and the non-universal functions F n ({η j,k }), G n ({η j,k }) are just constant. In the following section we will just omit these non-universal functions and just focus on the universal part.
There is one important comment we must address. To compute the four-point function of fluxed twist fields on the strip, we apply a conformal mapping from the strip to the UHP. This is somewhat similar to the case of computation of the negativity at finite temperature, where it is well-known that if the partial transposition involves an infinite part of an infinite system at finite temperature, using the conformal map from the cylinder to the complex plane is naively wrong [56]. This is not the case if one, for example, is interested in the negativity between two (adjacent or disjoint) finite intervals [57,58]. In our situation, we are concentrated on the charged Rényi negativity between two finite intervals. Thus it's indeed free of these troubles.

Evolution of charged logarithmic negativity in boundary state quench
In this section, we will consider the dynamics of the charged logarithmic negativity between two intervals after a global quench whose post-quench evolution is governed by a massless Hamiltonian.

Bipartite system
It's convenient to first study the case in which A = A 1 ∪ A 2 is the entire system. In this case, ρ A corresponds to a pure state. As explained in previous section, in this case, the temporal evolution of the charged Rényi negativity is governed by T 2 n,µ (w 1 )T 2 n,µ (w 2 ) on the strip. The strip two-point function can be computed from the one in the UHP which has the standard form where ∆ (2) n,µ = ∆ n,2µ , odd n 2∆ n 2 ,µ , even n (4.2) After mapping the correlator in the UHP to the one in the strip, we find In the spacetime scaling limit regime t Then it's straightforward to derive the charged Rényi logarithmic negativity In particular, for n = 1, the result is The charged logarithmic negativity is obtained by taking replica limit n e → 1 in E ne (µ) (4.7)
The temporal evolution of the charged logarithmic negativity is obtained by taking replica limit n e → 1 in E ne (µ). The result is (4.11) From eq. (4.9), we have where h 3 (µ) = − µ 2 2π 2 .
The charged logarithmic negativity is easily obtained as . As before, we also report the expression of E 1 (µ) here since it will be useful in the computation of charge imbalance resolved negativity. The result is (4.16) We stress that eq. (4.9) and eq. (4.14) are only valid for CFT in which there is a perfect linear dispersion. However this is not the case for the underlying lattice model, where the excitation has non-linear dispersion. We will discuss how to adapt eq. (4.9) and eq. (4.14) to describe the dynamics of charged Rényi (logarithmic) negativity in the mass quench protocol of the complex harmonic chain.

Charge imbalance resolved negativity
In this section, using the results obtained in the last section, we will compute the dynamics of the charge imbalance resolved negativity between two intervals after a global quench to a conformal Hamiltonian. According to the discussion in section 2, our strategy is to first compute Z T 2 (q) and p(q) (cf. eq. (2.19) and eq. (2.16)) from the charged moments of partially transposed RDM, and then the charge imbalance resolved negativity is given by eq. (2.19).

Bipartite system
In this time region, we have After Fourier transformation, we find where Erfi(x) is the imaginary error function A very similar expression exists for p(q) which we omit here. In the spacetime scaling limit, we have much more concise results Z T 2 (q) = 2e πt 2τ 0 t/τ 0 πq 2 + 4πt 2 /τ 2 0 (5.4) and the probability distribution is given by Thus the charge imbalance resolved negativity is obtained from eq. (2.19) In this time region, we can simply make the replacement t → l 1 /2 to obtain corresponding quantities. In particular, the charge imbalance resolved negativity is given by

Two adjacent intervals
Without loss of generality, we will assume l 1 < l 2 , the other case can be worked out similarly.

Two disjoint intervals
The computation of charge imbalance resolved negativity for two disjoint intervals is similar to previous subsection, here we just report the final results (the space time scaling limit has been taken). In this circumstance, the different size relation between l 1 , l 2 and d may lead to different behavior. However, for simplicity, we will only consider the case d < l 1 < l 2 .
At this early time stage, we have E(µ) = E 1 (µ), therefore N (q) = 0. (5.25) d/2 < t < l 1 /2 d/2 < t < l 1 /2 d/2 < t < l 1 /2 In this time period, we have (5.26) Then the charge imbalance resolved negativity is obtained straightforward as In this time region, the following results are obtained (in the spacetime scaling limit) (5.28) The charge imbalance resolved negativity is given by In this time period, we have (5.30) The charge imbalance resolved negativity is In this late time region, we have E(µ) = E 1 (µ) again, and charge imbalance resolved negativity is zero N (q) = 0. (5.32)

Total negativity
Having obtained the quench dynamics of the charge imbalance resolved negativity, it's an easy task to check whether we can recover the known results of quench dynamics of entanglement negativity. Using the formula  Figure 1: Numerical data of charged Rényi logarithmic negativity E 2 (µ) as a function of t/l for µ = π 8 , π 2 , π in the complex harmonic chain. The full lines are the CFT predictions (cf. eq. (4.9) and eq. (4.14)). Left panel: Adjacent interval with L = 3000, l 1 = 200, l 2 = 300. Right panel: Disjoint interval with L = 3000, l 1 = 200, l 2 = 300, d = 100. We have used the best fitted value of τ 0 . and the reconstruction formula of negativity eq. (2.10), it's easy to check that indeed in each time period, we can obtain the total entanglement negativity from the charge imbalance resolved one, i.e. in the spacetime scaling limit, we have the following relation In this way, we have recovered the known quench dynamics [50] of the total logarithmic negativity from the charge imbalance resolved ones, providing evidence of the correctness of our results.

Numerical test
In this section, we test our analytical predictions against numerical computation in the complex harmonic chain which is the lattice version of our complex Klein-Gordon field theory. We will use the correlation matrix techniques to obtain the charged Rényi (logarithmic) negativity.
In section 3, we have already reported the Hamiltonian and the U (1) conserved charge of the complex harmonic chain (cf. eq. (3.9) and eq. (3.10)). From the dispersion relation, it's easy to see that the Hamiltonian has zero modes for k = 0 and m = 0, and the group velocities are obtained from the dispersion relation as The maximum velocity v max ≡ max p v p determine s the spreading of entanglement and correlations.
In the global quench protocol, for t < 0, the system is prepared in the ground state of the Hamiltonian At t = 0, the system evolves unitarily according to a new Hamiltonian H CHC (m) with a different value m, namely |ψ(t) = e −iH CHC (m)t |ψ 0 . The dynamics of symmetry resolved negativity after a global quench described in the previous section corresponds to set the new parameter m = 0. In this section, we will focus on the global quench to a massless Hamiltonian (i.e. m = 0) to test our CFT predictions. In the next section we will discuss the case m = 0.  Figure 2: Numerical data of charged logarithmic negativity E(µ) as a function of t/l for µ = π 6 , π 3 , π in the complex harmonic chain. The full lines are the CFT predictions (cf. eq. (4.10) and eq. (4.15)). Left panel: Adjacent interval with L = 2000, l 1 = 200, l 2 = 250. Right panel: Disjoint interval with L = 5000, l 1 = 250, l 2 = 400, d = 100. We have used the best fitted value of τ 0 .
Since the Hamiltonian of a complex harmonic chain is just two copies of a real harmonic chain's Hamiltonian, we could focus on the real harmonic chain first, and take the double copies into account only at the end of the calculations. In the correlation matrix method of computing entanglement measures, we first need to know the following two-point correlators of the real scalars The explicit form of these correlators is given by where X k (t) = 1 e k e k e 0,k cos 2 (e k t) + e 0,k e k sin 2 (e k t) , P k (t) = e k e k e 0,k sin 2 (e k t) + e 0,k e k cos 2 (e k t) , M k (t) = e k e 0,k − e 0,k e k sin(e k t) cos(e k t). (6.6) From the above expressions, it's easy to find that the contribution from the zero mode m = 0 and k = 0 is finite The evolution of charged Rényi (logarithmic) negativity of any subsystem A containing l lattice sites can be computed from these time-dependent correlators. Firstly, we should consider the correlation matrices X A (t),P A (t) and M A (t) obtained by restricting the the indices of the corresponding correlation matrices to the sites belonging to A. Given X A ,P A and M A , the covariance matrix Γ A and the symplectic matrix J A associated to the subsystem A are where 0 l and I l are l × l zero matrix and identity matrix respectively. Then find the eigenvalues of the 2l × 2l matrix iJ A Γ A (t) which we denoted it by {±σ 1 (t), · · · , ±σ l (t)}. It's also convenient to introduce the Fock space basis |n ≡ ⊗ l j=1 |n j , defined by products of eigenstates of the number operator in the subsystem A, the reduced density matrix of A (in a real harmonic chain) can be written as In the Fock basis {|n }, Q T 2 2 = Q 2 and the operator Q A = Q 1 −Q T 2 2 = Q 1 −Q 2 becomes exactly the charge imbalance operator. For the complex harmonic chain, the charged Rényi negativity factorises as (6.10) In bosonic system, the net effect of partial transposition with respect to A 2 is changing the sign of the momenta corresponding to A 2 . Thus the momenta correlators in the partial transposed density matrix can be obtained from P A by simply change the sign of the momenta that in A 2 , i.e. (6.11) If we denote the eigenvalues of iJ A Γ A (t) T 2 by {±τ 1 (t), ±τ 2 (t), · · · , ±τ l (t)}, then the charged Rényi logarithmic negativity is given by 1 6.12) and the charged logarithmic negativity is The numerical data of the dynamics of the charged Rényi (logarithmic) negativity are shown in fig. 1 and fig. 2, in which the CFT predictions are drawn with the full line for comparison.

Global mass quench of the harmonic chain
In this section, we consider a global quantum quench in which the complex harmonic chain is initially prepared in the ground state with mass m 0 and at time t = 0 the mass is quenched to a different (non-zero) value m = m 0 .  Figure 3: Charged Rényi logarithmic negativity E 2 (µ) as a function of t/l for different µ after the mass quench from m 0 = 1 to m = 2 in the complex harmonic chain. The full lines are the quasi-particle predictions (cf. eq. (7.10) and eq. (7.13)). Left panel: Adjacent intervals with L = 2500, l 1 = 200, l 2 = 300 and for µ = π 8 , π 3 , π 2 . Right panel: Disjoint intervals with L = 2500, l 1 = 150, l 2 = 200, d = 50 and for µ = π 10 , π 5 , π 2 . As shown in the figure, the agreement is extremely excellent.

Quasi-particle picture for entanglement entropy
In free scalar theory, since the particle numbers in each momentum mode are conserved quantities, the time evolution leads to local and quasi-local observables converging to their average values in the Generalized Gibbs Ensemble (GGE) instead of ordinary Gibbs ensemble.
The quasi-particle picture of the time evolution of entanglement after a global quantum quench has been proposed in [45], The basic idea is that one can view the initial state as a source of quasi-particle excitations. Pairs of particles emitted from the same point in space are highly entangled whereas particles produced from points far apart are incoherent. We assume that the pairs of quasi-particle are created uniformly with opposite momenta (p, −p). After the production, these quasi-particles travel ballistically with velocity v p = −v −p . The entanglement entropy and Rényi entropy of subsystem A are proportional to the pairs of entangled quasi-particle shared with its complement at a given time t.
In free models, we have S n (t) = dp 2π s GGE (p) is momentum space density of the Rényi entropies in the GGE thermodynamic state [59]. Recall that the well-known result of the CFT prediction of the quench to massless dynamics of Rényi entropies S n (t) = 1 n − 1 π∆ n,0 2τ 0 min(2t, l).
We can formally obtained the formula eq. (7.1) from the CFT prediction eq. (7.2) by replacing t → |v p |t, − π∆n,µ 2τ 0 → s 7.2 Quasi-particle picture for charged logarithmic negativity The complex harmonic chain is a free models and the stationary behavior of local and quasi-local observables is described by the GGE  Figure 4: Charged logarithmic negativity E(µ) as a function of t/l for different µ after the mass quench from m 0 = 1 to m = 2 in the complex harmonic chain. The full lines are the quasi-particle predictions (cf. eq. (7.11) and eq. (7.14)). Left panel: Adjacent intervals with L = 2500, l 1 = 200, l 2 = 300 and for µ = π 8 , π 4 , π. Right panel: Disjoint intervals with L = 2500, l 1 = 200, l 2 = 300, d = 150 and for µ = π 8 , π 4 , π. As shown in the figure, the agreement is perfect.
where λ k are Lagrange mulitpliers and Z is the normalisation constant such that Trρ GGE = 1. We have The mode occupation numbern (7.7) We have find the density of the logarithmic charged moment in momentum space n,µ (k) = −2 log |(1 + n k ) n − e iµ n n k |. Then we have all the ingredients to give conjectures about the dynamics of charged Rényi (logarithmic) negativity after a global mass quench in our free lattice model.

Conclusion
In this paper, we discuss the temporal evolution of charge imbalance resolved negativity after two types of global quenches. Firstly, we have considered the boundary state quench in which the post-quench dynamics is governed by a conformal Hamiltonian. Evaluating the correlators of the fluxed twist fields in the UHP, we obtained the quench dynamics of charged Rényi negativity. Then the charge resolved negativity is obtained by Fourier transformation. The total negativity can be reconstructed from these charge resolved ones. We also discussed the mass quench in the underlying lattice model and made conjectures based on the quasi-particle picture.
It would be very interesting to investigate the evolution of other charge resolved entanglement measures in the same quenches discussed in this paper, or studying other quench setups such as local joining quench or local operator quenches in field theories and in holographic theories. One may find the breakdown of quasi-particle picture in certain circumstance [60,61]. We will turn to these problems in the future work.