Latest revision as of 13:25, 4 June 2024

Abstract

For open quantum systems, a short-time evolution is usually well described by the effective non-Hermitian Hamiltonians, while long-time dynamics requires the Lindblad master equation, in which the Liouvillian superoperators characterize the time evolution. In this paper, we constructed an open system by adding suitable gain and loss operators to the Chern insulator to investigate the time evolution of quantum states at long times by numerical simulations. Finally, we also propose a topolectrical circuits to realize the dissipative system for experimental observation. It is found that the opening and closing of the Liouvillian gap leads to different damping behaviours of the system and that the presence of non-Hermitian skin effects leads to a phenomenon of chiral damping with sharp wavefronts. Our study deepens the understanding of quantum dynamics of dissipative system.

Keywords: Open quantum system, chiral damping, topolectrical circuits

1. Introduction

With the laboratory advances in modulating dissipation and quantum coherence,the theory of open and nonequilibrium systems has received renewed attention [1,2]. Non-Hermitian Hamiltonians have been used to describe a large number of non-conservative systems, such as classical waves with gain and loss [3-8], solids with finite quasiparticles lifetimes [9-11], and open quantum systems [12-14]. The unique features of non-Hermitian systems have been recognized in a variety of physical settings, in particular the non-Hermitian skin effect (NHSE) [15,16], where the eigenstates of the system are exponentially localized on the boundary. In recent years, the impact of NHSE has been extensively studied [17-29].

NHSE was also found in open quantum systems [30]. For open quantum systems, the non-Hermitian effective Hamiltonian describes the time evolution of the wavefunction under post-selection conditions, while the time evolution of the density matrix (without post-selection) is driven by the Liouvillian superoperator in the master equation [2,31-33]. It has been found that the Liouvillian superoperator can also exhibit non-Hermitian skin effects and that such effects can significantly affect the dynamical behaviour of the system at long times [30,34-43]. In a large class of open quantum systems, the quantum state in the long time limit converges to the steady state by algebraic damping under periodic boundary conditions and exponential damping under open boundary conditions [30].

In recent years, it has been discovered that topolectrical circuits can be used as platform to simulate the lattice systems, thus enabling the study of topological states in topolectrical circuits and gradually developing the field of topological circuitry [44-46]. Some of the early experiments and theories were extensively studied in Hermitian systems [45,47]. Since the phenomena of non-Hermitian systems are more rich than that of Hermitian systems, increasing attentions are contributed into the non-Hermitian physics, and some interesting phenomena have also been realized by topolectrical circuits [48-52].

Previous studies on open quantum dynamics and topolectrical circuits have mainly focused on one-dimensional non-Hermitian models, and relatively few studies on higher-dimensional non-Hermitian models. In this paper, we consider a two-dimensional open quantum system based on Chern insulators. Following the method developed in Song et al. [30], we study the dynamics of this system in terms of the damping matrix derived from the Liouvillian superoperator, and give a model of topolectrical circuit realization of the damping matrix based on Kirchhoff’s theory. It is found that due to the NHSE of the damping matrix, the long-time dynamics of the system under open boundary conditions is significantly different from that under periodic boundary conditions.

Our paper is organized as follows: in section 2, we briefly review the general framework on how to convert Liouvillian operators with linear jumps to non-Hermitian damping matrix. In sections 3 and 4, we compute and numerically simulate the long-time evolution of the model. In section 5, we give the circuit model of the non-Hermitian damping matrix . Finally, we conclude in section 6.

2. General formalism of damping matrix

An open quantum system undergoing Markovian damping satisfies the Lindblad master equation

where ${\displaystyle \rho }$ is the density matrix of the system, ${\displaystyle H}$ is the Hamiltonian that represents unitary evolution of the system, and ${\displaystyle L_{\mu }}$ are Lindblad dissipation operators describing the quantum jumps induced by the coupling to the environment. The above equation can be abbreviated as ${\displaystyle {\frac {d\rho }{dt}}=L\rho }$, where ${\displaystyle L}$ is called the Liouvillian superoperator. By regarding the density matrix ${\displaystyle \rho }$ as a vector that consists of matrix elements ${\displaystyle {\rho }_{i,j}}$, ${\displaystyle L}$ is represented as a matrix whose elements are given by [53]

These representations enable one to treat the Lindblad equation as a linear equation. In other words, the dynamics of the system can be understood in terms of the eigenvalue problem of the Liouvillian matrix: ${\displaystyle L{\rho }^{\left(i\right)}={\lambda }_{i}{\rho }^{\left(i\right)}.}$ The Hamiltonian and dissipators can be expressed in terms of 2n Majorana fermions [54]

where ${\displaystyle {\gamma }_{i}}$ are Majorana fermions satisfying ${\displaystyle \lbrace {\gamma }_{i},{\gamma }_{j}\rbrace =2{\delta }_{ij}}$. The matrix ${\displaystyle H^{M}}$ is chosen to be an antisymmetric matrix, ${\displaystyle {\left(H^{M}\right)}^{T}=-H^{M}}$. Defining ${\displaystyle M_{ij}={\sum }_{\mu }l_{\mu ,i}^{\ast }l_{\mu ,j}^{}}$, ${\displaystyle M_{{}_{ij}}^{M}={\sum }_{\mu }{\left(l_{\mu ,i}^{M}\right)}^{\ast }l_{\mu ,j}^{M}}$, we have ${\displaystyle M_{{}_{}}^{M}={\frac {1}{4}}M\otimes (1+{\sigma }_{y})}$. Under the third quantization [54,55], the Liouvillian superoperator is expressed as a quadratic form of the 2n complex fermions (4n Majorana fermions)

where ${\displaystyle Z=H^{M}+iRe{\left(M^{M}\right)}^{T},Y=2Im{\left(M^{M}\right)}^{T}}$, and ${\displaystyle c=(c_{1},c_{2},...,c_{2n})}$ are third quantized complex fermions. Through the above expression, we can obtain the Liouvillian eigenspectrum [54,55]

with ${\displaystyle v_{i}\in \lbrace 0,1\rbrace }$, where ${\displaystyle \left\{E_{i}\right\}}$ is the eigenspectrum of ${\displaystyle 4iZ}$. Here ${\displaystyle \lambda }$ contains valuable information of the full density-matrix dynamics, and it can be easily obtained from the damping matrix ${\displaystyle X}$ with ${\displaystyle X_{ij}=ih_{ji}-{\sum }_{\mu }l_{\mu ,j}^{\ast }l_{\mu ,i}^{}}$ [36]. Rewriting ${\displaystyle M}$ as ${\displaystyle M=M_{r}+iM_{i}}$, where ${\displaystyle M_{r},M_{i}}$ are real matrices, we have ${\displaystyle M^{M}={\frac {1}{4}}(M_{r}+iM_{i})\otimes (1+{\sigma }_{y})}$. ${\displaystyle Z}$ can be further written as

Therefore,

The eigenvalue of ${\displaystyle 4iZ}$ are the union of the eigenvalues of ${\displaystyle X}$ and ${\displaystyle X^{\ast }}$, which gives the Liouvillian eigenspectrum.

Then we outline the general form of the Lindblad damping matrix in open quantum systems [30]. We consider tight-binding models whose Hamiltonian can generally be written as ${\displaystyle H=\sum {}_{ij}h_{ij}c_{i}^{\dagger }c_{j}}$, where ${\displaystyle c_{i}^{\dagger },c_{i}}$ are the creation and annihilation operators on lattice site ${\displaystyle i}$, and ${\displaystyle h_{ij}=h_{ij}^{\ast }}$ is the hopping amplitude between the lattice points of the system (${\displaystyle i\not =j}$) or onsite potential (${\displaystyle i=j}$). It is convenient to define the single-particle correlation function ${\displaystyle {\Delta }_{ij}(t)=Tr[c_{i}^{\dagger }c_{j}\rho (t)]}$ to observe the time evolution of the density matrix. Each cell is coupled to the environment through the gain jump operator ${\displaystyle L_{\mu }^{g}=\sum {}_{i}D_{\mu i}^{g}c_{i}^{\dagger }}$ and loss jump operator ${\displaystyle L_{\mu }^{l}=\sum {}_{i}D_{\mu j}^{l}c_{i}}$. Substituting the Lindblad quantum master equation into the time evolution of the single-particle correlation function, we can obtained

where ${\displaystyle X=ih^{T}-(M_{l}^{T}+M_{g})}$ is the damping matrix with ${\displaystyle {\left(M_{g}\right)}_{i}{}_{j}=\sum {}_{\mu }D_{\mu i}^{g\ast }D_{\mu j}^{g}}$ and ${\displaystyle {\left(M_{l}\right)}_{i}{}_{j}=\sum {}_{\mu }D_{\mu i}^{l\ast }D_{\mu j}^{l}}$. The steady state correlation ${\displaystyle {\Delta }_{s}=\Delta (\infty )}$, to which the long-time evolution of any initial state converges, is determined by ${\displaystyle d{\Delta }_{s}/dt=0}$ or ${\displaystyle X{\Delta }_{s}+{\Delta }_{s}X^{\dagger }+2M_{g}=0}$. Focusing on the deviation towards the steady state ${\displaystyle {\tilde {\Delta }}(t)=\Delta (t)-{\Delta }_{s}}$, whose time evolution is ${\displaystyle d{\tilde {\Delta }}(t)/dt=X{\tilde {\Delta }}(t)+{\tilde {\Delta }}(t)X^{\dagger }}$, we can integrate it with Eq. (1) to obtain

Therefore, the dynamical behaviour of the system can be characterized by the damping matrix.

3. Model

In this paper, we consider the Chern insulator model with the Hamiltonian in momentum space as

where ${\displaystyle {\epsilon }_{k}=m+t_{x}\cos k_{x}+t_{y}\cos k_{y}}$. Let each unit cell contain a single loss and gain dissipator,

where ${\displaystyle x}$ denotes the lattice site, ${\displaystyle A,B}$ refer to the sublattice. The Fourier transformation of ${\displaystyle X}$ is ${\displaystyle X(k)=ih^{T}(-k)-M_{l}^{T}(-k)-M_{g}(k)}$. The gain and loss dissipators are intra-cell, so these ${\displaystyle M(k)}$ matrices are independent of ${\displaystyle k}$, ${\displaystyle M_{l}(k)={\frac {\sqrt {2}}{2}}\lambda +{\frac {1}{2}}{\sigma }_{x}-}$${\displaystyle {\frac {1}{2}}{\sigma }_{y},M_{g}(k)={\frac {\sqrt {2}}{2}}\lambda +}$${\displaystyle {\frac {1}{2}}{\sigma }_{x}+{\frac {1}{2}}{\sigma }_{y}}$. Then, the damping matrix in momentum space is

It can be written in the form of left and right eigenvectors,

where ${\displaystyle X^{\dagger }\vert u_{Ln}\rangle ={\lambda }_{n}^{\ast }\vert u_{Ln}\rangle ,X\vert u_{Rn}\rangle =}$${\displaystyle {\lambda }_{n}\vert u_{Rn}\rangle }$. It is worth noting that our ${\displaystyle M_{l}}$ and ${\displaystyle M_{g}}$ satisfy ${\displaystyle M_{l}^{T}+M_{g}=2M_{g}}$, guaranteeing that ${\displaystyle {\Delta }_{S}={\frac {1}{2}}I_{2L\times 2L}}$ is a steady state solution, where ${\displaystyle L=N_{x}\times N_{y}}$, ${\textstyle L}$ is the system size, and ${\displaystyle N_{x},N_{y}}$ are the size in ${\textstyle x,y}$ direction, respectively. We assume that the initial state of the system is the completely filled state, i.e., ${\displaystyle \Delta (0)}$ is an identity matrix. Therefore, Eq.(9) can be re-expressed as

According to the dissipative property, ${\displaystyle Re\lbrace {\lambda }_{n}\rbrace \leq 0}$ always holds. The Liouvillian gap ${\displaystyle \Lambda =\min[2Re(-{\lambda }_{n})]}$ plays a decisive role in long-time dynamics. The opening gap (${\displaystyle \Lambda \not =0}$) implies an exponential rate of convergence to the steady state, while the closing gap (${\displaystyle \Lambda =0}$) implies algebraic convergence [34].

4. Chiral damping

For simplicity, the parameters of our model are taken as ${\displaystyle l_{x}=l_{y}=1}$, ${\displaystyle t_{x}=t_{y}=-1}$. We first study the dynamical behaviour under the periodic boundary conditions. Diagonalizing ${\displaystyle X(k)}$, we obtain the energy spectrum as shown in Figure 1. It is found that the Liouvillian gap vanishes at ${\displaystyle m=1.5}$, while the gap opens at ${\displaystyle m=2.5}$. So we expect the damping rate to be algebraic and exponential in each case, respectively.

Figure 1. Eigenvalues of the damping matrix X. Blue: periodic boundary; Red: open boundary . The Liouvillian gap under periodic boundary condition vanishes for (a) and (b), while it is nonzero for (c) and (d). Under open boundary condition, the Liouvillian gap is nonzero in all four cases. This significant difference between open and periodic boundary comes from the NHSE of ${\textstyle X}$. (a) ${\displaystyle \lambda =0.1,m=1.5}$. (b) ${\displaystyle \lambda =0.5,m=1.5}$. (c) ${\displaystyle \lambda =0.1,m=2.5}$. (d) ${\displaystyle \lambda =0.5,}$ ${\displaystyle m=2.5}$

To verify this, we define the site-averaged fermion number deviation from the steady state ${\displaystyle R(t)={\sqrt {{\frac {1}{N_{X}N_{Y}}}{\sum _{x}R_{x}(t)}^{2}}}={\sqrt {{\frac {1}{N_{X}N_{Y}}}{\sum _{x}\left({\frac {n_{x}(t)-n_{x}(t-{\delta }_{t})}{{\delta }_{t}}}\right)}^{2}}}}$, where ${\displaystyle R_{x}(t)=n_{x}(t)-n_{x}(\infty )}$, and ${\displaystyle n_{x}(t)={\Delta }_{xA,xA}(t)+{\Delta }_{xB,xB}(t)}$. The numerical results are shown in Figure 2. As anticipated, it is observed that the damping of ${\textstyle R(t)}$ is algebraic for cases black and red lines with ${\textstyle m=1.5}$, while exponential for blue and green lines with ${\textstyle m=2.5}$ under the periodic boundary condition.

Figure 2. Damping of site-averaged fermion number towards the steady state under periodic boundary condition with size ${\displaystyle L=30\times 30}$. ${\textstyle m=1.5}$ (black and red) exhibits a slow algebraic damping, while ${\textstyle m=2.5}$ (blue and green) is an exponential damping. The initial state is the completely filled state

Next we turn to the open boundary conditions. Since the damping matrix ${\textstyle X}$ has NHSE, its energy spectrum is no longer that of the periodic boundary conditions. At this point all the energy spectrums have a non-zero energy gap (red part of Figure 1), therefore, we expect an exponential long-time damping of ${\displaystyle {\tilde {\Delta }}(t)}$. The numerical simulation in Figure 3 confirms this exponential behaviour with ${\textstyle R(t)}$ having a period of algebraic damping before entering into the exponential damping. The time of the algebraic damping increases with the size ${\textstyle L}$ (Figure 3(a)). To better understand this feature, we plot the damping in several unit cells in the same x dimension (${\textstyle ix=1}$), as shown in Figure 3(b). It can be seen that the left end (${\textstyle iy=1}$) enters the exponential damping immediately, and the other sites enter the exponential damping in turn according to their different distances to the left end.Due to a process of algebraic damping that occurs before entering the exponential stage,there is a "damping wavefront" from left (${\displaystyle ix=1,iy=1}$) to right (${\displaystyle ix=1,iy=N_{y}}$). This phenomenon is known as "chiral damping".

Figure 3. (a) Site-averaged particle number damping under periodic boundary conditions (solid line) and open boundary conditions (dashed line) for several sizes ${\textstyle L}$. The long-time damping of ${\textstyle R(t)}$ follows a power law under periodic boundary condition, while the damping follows an exponential law after an initial power law stage under open boundary condition. (b) Particle number damping on several sites. The system size is ${\textstyle 30\times 30}$, and the left end (${\displaystyle ix=1,iy=1}$) enters the exponential phase from the beginning, followed by the other sites in turn. For (a) and (b), the initial state is completely filled state, ${\displaystyle m=1.5,\lambda =0.1}$

The phenomenon of chiral damping can be observed more intuitively as shown in Fig. 4(a) where the colour shades indicate the value of ${\textstyle R(t)}$. Under the periodic boundary condition, the time evolution follows a slow power law while under the open boundary condition, a wavefront moving to the upper right is observed. This can be intuitively linked to the phenomenon that all eigenstates of ${\textstyle X}$ are localized in the upper right corner, which arises from the non-Hermitian skin effect of the damping matrix ${\textstyle X}$. If the matrix ${\textstyle X}$ does not have NHSE under the open boundary condition, the fermion number of the system should have a similar behaviour of damping under different boundary conditions. Therefore, the non-Hermitian skin effect plays an important role in open quantum systems and significantly affects the dynamical behaviour of open quantum systems.

Figure 4. Evolution of ${\textstyle R(t)}$ at each lattice site under open boundary conditions (a) and periodic boundary conditions (b)

5. Experiment realized

Next we give the scheme of topolectrical circuits to simulate the damping matrix. Based on the similarity between the Kirchhoff equation and the Schrödinger equation, it is possible to simulate the Hamiltonian of the system using different circuit components, and the different parameters in the Hamiltonian can be adjusted independently by various components. The circuit Laplacian corresponding to the Hamiltonian can be written as

where ${\textstyle W}$ and ${\textstyle D}$ are diagonal matrices containing the total conductance from each node to the ground and to the rest of the circuit, respectively. ${\textstyle C}$ is the adjacency matrix of conductances [44].

Figure 5. Structure of topolectrical circuit to realize damping matrix under periodic boundary conditions. (a) Connection relations between the nodes. The blue solid line box containing two “sublattice” nodes A (red) and B (blue) simulates a unit cell of ${\textstyle X}$. The black (grey) solid line indicates the coupling between nodes in the ${\textstyle x(y)}$-direction. (b) Circuit element structure is detailed for the green dashed framed rectangle in (a). (c) Internal circuit diagram of the INIC element, consisting of an operational amplifier and impedances ${\displaystyle Z_{\pm }}$. The impedance ${\displaystyle Z}$ is the target element, and different conductance in different directions of ${\displaystyle V_{i,r}}$ can be achieved by connecting the INIC in series. ${\displaystyle Z_{\pm }}$ satisfies ${\displaystyle Z_{+}=Z_{-}}$. (d) Grounding module of the nodes. The resistances ${\displaystyle R_{A,}R_{B}}$ and capacitance ${\textstyle C}$ are used to simulate the onsite potential, and inductance L allows the Laplacian eigenvalue spectrum to be shifted uniformly as desired

Considering the periodic boundary conditions first, the topolectrical circuit for realizing the damping matrix X is illustrated in Figure 5. Figure 5 depicts the schematic diagram of the overall circuit structure, which gives the connection relationship between the nodes. Figure 5 shows the detailed circuit component of the unit which is the green dashed box in Figure 5(a). The blue box in Figure 5(a) represents a unit cell in the system, and the two nodes inside it correspond to sublattices A (red) and B (blue). The circuit connections in the x and y directions are distinguished by black and gray. From Figure 5(b) we can obtain the matrices C and D in Eq. (9), so that

with ${\displaystyle C_{1}=-\lambda }$, ${\displaystyle C_{2}=l_{x}/2}$, ${\displaystyle {\frac {1}{R_{0}}}={\frac {1}{R_{1}}}-2({\frac {1}{R_{3}}}+{\frac {1}{R_{4}}})}$, ${\displaystyle R_{1}=-1/\left(\omega \lambda \right)}$, ${\displaystyle R_{2}=-R_{3}=-2/t_{x}}$, ${\displaystyle R_{4}=-R_{7}=-2/t_{y}}$, ${\displaystyle R_{5}=-R_{6}=-2/l_{y}.}$

Comparing it with the damping matrix, we need to add grounding elements to match the onsite potential. The grounding elements of nodes A and B are shown in Figure 5(d), where the resistors ${\displaystyle R_{A,}R_{B}}$ and capacitors ${\displaystyle C}$ simulate the lattice potential, and ${\displaystyle R_{A,}R_{B}}$ satisfies ${\displaystyle R_{A}=-R_{B}}$. So the diagonal matrix ${\textstyle W}$ is