Lattice Lindblad simulation

We perform the real-time lattice simulation of an open quantum system, which is based on the Schwinger-Keldysh path integral representation of the Lindblad formalism. Although the real-time simulation generally suffers from the sign problem, there exist a few exceptional cases. We focus on a sign-problem-free system of a non-relativistic spinless fermion and analyze time evolution under driving and dissipation.


Introduction
Real-time simulation is one of the ultimate goals of lattice gauge theory. Since the real-time path integral has a complex-valued probability, the Monte Carlo sampling does not work due to the sign problem. Theoretical physicists have proposed many candidates for the solution, e.g., the complex Langevin method [1][2][3], the Lefschetz thimble method [4][5][6], and the tensor network approach [7][8][9]. Algorithms on quantum computers are also being discussed [10][11][12][13][14][15][16]. Although these proposals are hopeful, their applications are currently limited to systems with small numbers of degrees of freedom. Large-scale simulation of non-equilibrium quantum phenomena still remains an open issue.
The Lindblad master equation describes the time evolution of open quantum theory [17,18]. A system gains or loses internal quantities, such as energy-momentum or particle number, through the interaction with an environment. In general, the path integral representation of the Lindblad equation is not positive definite. In very special cases, however, it does not encounter the sign problem. The known examples are purely dissipative systems [19][20][21][22] and a purely fermionic system [23,24]. The fermionic system is of particular interest because it can contain both of unitary and dissipative dynamics. If one can find the system where the fermion determinant is positive definite as a consequence of miracle cancellation, its non-equilibrium property is accessible by the real-time lattice simulation.
In this paper, we study the Lindblad evolution of another fermionic model, which is different from the models in earlier works [19][20][21][22][23][24]. The lattice formulation is based on the Schwinger-Keldysh path integral representation, which is reviewed in Sec. 2. As explicitly shown in Sec. 3, the model has a significant property; the fermion determinant is exactly quenched and the sign problem is absent. Owing to this property, the real-time lattice simulation is possible. We show simulation results in Sec. 4 and discuss theoretical aspects of the model in Sec. 5.

Path integral formalism of the Lindblad evolution
Let us briefly summarize the path integral representation of the Lindblad formalism. A comprehensive review can be found in the literature [25]. The Lindblad equation was originally formulated in terms of a density matrix in a doubled Hilbert space. In the path integral representation, it is translated to two branches of the real-time axis: one in the forward direction and the other in the backward direction. They draw a closed counter, which is called the Schwinger-Keldysh path, in the time plane. The counter is discretized with the lattice spacing δ. The fermion field ψ contains three pieces on the discretized time path in Fig. 1. Fig. 1 Closed path in the real and imaginary time plane. The red, blue, and white circles stand for ψ + , ψ − , and ψ 0 , respectively. The anti-periodic boundary condition is imposed at the ends of the imaginary-time direction.
The Schwinger-Keldysh path integral is given by For the sake of brevity, we consider a non-relativistic one-component fermion on the threedimensional spatial lattice. The continuous-time limit δ → 0 is assumed to be taken eventually, while the continuous-space limit is not considered here. The time evolution at t > 0 is governed by the real-time action The Hamiltonian H ± determines unitary dynamics of ψ ± . The second line represents dissipations caused by interactions with an environment. The index k labels the species of the dissipations, γ k are 2/8 their strengths, and L k± are the so-called jump operators. The initial condition at t = 0 is controlled by the imaginary-time action with the same Hamiltonian as H ± but of ψ 0 . This means that the initial state is a thermal state at the temperature T = 1/(N τ δ).
The expectation value of an operator is given by the average of the two branches, where O + and O − are the corresponding operators on the upper and lower branches. On the other hand, the difference is identically zero for any physical operator [25]. The average and difference are called the "classical" and "quantum" parts of the observable, respectively.

Model without the sign problem
For the unitary part, we consider the tight-binding model on a three-dimensional cubic lattice, where w is the hopping parameter and e k is the unit lattice vector in the x k direction (k = 1, 2, 3). The external field A k (x) couples to the electric current. We choose the electric current operators as the jump operators When the dissipation is isotropic, γ 1 = γ 2 = γ 3 ≡ γ, the dissipation term is written as the completed square form Here we used the Hermiticity j k± = j * k± . After the Hubbard-Stratonovich transformation, the dissipation term is rewritten by the auxiliary field B k (x) as and 3/8 Thus, the time evolution is unitary, but it is subject to the Gaussian dephasing noise. Since the fermion action is in the bilinear form ψ * Dψ, we can integrate out the fermion field, and then get When we choose the backward difference for the time derivative, where h 0 and h n (n = 1, · · · , N t ) are the matrix representation of Eqs. (8) and (13), respectively. The matrix structure in the spatial directions is omitted in Eq. (15). The fermion determinant is given by When N τ is an even number, this is positive definite because The Monte Carlo sampling is free from the sign problem.
The fermion determinant has a significant property. We can evaluate the determinant (16) as This is a B k -independent constant in the continuous-time limit δ → 0. (Note that the auxiliary field action in Eq. (14) is O(δ).) This property can be more clearly understood in another discretization manner. Following the discretization in the time-ordered form [26], we can deform Eq. (15) as The two discretizations (15) and (19) are equivalent up to the higher order of δ, which is irrelevant for the continuous-time limit. In this form, the determinant is This is completely independent of B k . We do not need to compute the fermion determinant in the Monte Carlo simulation. This looks like the famous "quenched approximation". In this model, however, the quenching is not approximate but exact. This property comes from the specialty of the dissipation term. Through the Hubbard-Stratonovich transformation, the difference j k+ − j k− in the dissipation term is replaced by B k . As mentioned in Sec. 2, the quantum part j q is zero so B k must be zero. If the determinant had a linear term proportional to B k , B k would be nonzero. Therefore, the leading B k -dependent term is not O(δB k ) but O(δ 2 B 2 k ) as long as the determinant is a polynomial of the dimensionless combination δB k . This argument holds true for any bilinear Hermitian jump operator written as the completed square form (10).

Simulation results
We show several results obtained by the simulation. As mentioned above, the fermion determinant can be quenched and the Monte Carlo weight is the Gaussian distribution. This property greatly reduces simulation cost. The simulation is just a white noise sampling and any update algorithm is not necessary. The hopping parameter is fixed at w = 0.1/δ. All the dimensional quantities in this section are scaled in the unit of w. Spatial lattice volume is N 3 s = 8 3 and imaginary-time length is N τ = 4. Since the path integral is the grand canonical ensemble constructed from coherent states, a particle number is defined as an expectation value. We numerically checked that the expectation value of the particle number densitȳ is conserved in the real-time evolution. The time-splitting form of Eq. (21) originates from a chemical potential on the lattice [27]. Im P(t) wt γ / w = 0 γ / w = 1 Fig. 2 The real-time propagator. The real-time length is N t = 32.

5/8
A fundamental quantity of interest is the fermion propagator in the real-time direction. The propagator on the positive branch is plotted in Fig. 2. The external field is zero, A k (x) = 0, here. Without the dissipation γ = 0, real-time evolution is oscillatory. Since the initial state is not an eigenstate of the Hamiltonian, the propagator is given by a superposition of many oscillations. With the dissipation γ = 0, the propagator damps in time. Quantum coherence is lost under the Gaussian dephasing noise, i.e., quantum fluctuation of the electric current. The path integral is a functional integral over all possible classical paths. Even though the total path integral has space-time translation symmetry, each path can break the symmetry. Individual particles exchange energies and momenta with the environment, and gradually lose their original information. This picture can be seen from the occupation number The momentum integral of the occupation number gives the particle number density,n c = p n c (t, p), so the integral is conserved while the distribution is not. As shown in Fig. 3, the occupation number changes from the initial thermal distribution to a uniform one as time goes by. The original information of the initial state is lost by the dissipation. We can also study the time evolution of physical quantities. We set the external field such that the constant electric field E drives the electric current in the x 1 direction. The strength of the electric field is set E/w 2 = 1. The total energy per volume 6/8 and the electric current in the x 1 direction are plotted in Fig. 4. The factor −i is multiplied to make ε q and j q real. As the electric field is constantly applied, the electric current j c is induced and linearly increased. The system energy ε c is also increased. The energy gain is purely due to the electric field because only the external field explicitly breaks time translation symmetry and other terms, including the dissipation term, do not violate energy conservation. The quantum parts ε q and j q are always zero within the margin of error as expected.

Discussion
Because of the special property of the fermion determinant, the dissipation term is equivalent to the Gaussian quantum fluctuation. The resultant theory is very close to the stochastic Schrödinger equation. Particles are stochastically scattered and result in random trajectories. In this model, the stochastic noise couples to the current, i.e., momentum, so the momentum distribution changes in time while the total momentum is conserved. It is known that such an environment-induced diffusion leads to an efficient quantum transport [28,29]. Our formulation serves as an efficient numerical method for its large-scale simulation. The real-time propagator is a damping function. Since non-local correlation functions are given by the products of the propagator, all the non-local-time correlators will be damping functions, too. On the other hand, local-time observables of conserved quantities, such as ε c (t) and j c (t), are not damped because the dissipation term does not couple to these observables. The dissipation couples only to the non-observable j q . Non-equilibrium steady states cannot be realized by the competition between the dissipation and the driving in this model. In general, bilinear and Hermitian jump operators can realize only trivial steady states, where the density matrix is unity [30,31]. We need to discover some other models without the sign problem to study the formation of non-trivial steady states.
As mentioned in Sec. 3, our method works for any bilinear Hermitian jump operators, although we chose the current operator in this model. We can generally examine the evolution of an open 7/8 free-fermion theory under the dephasing process that randomizes the phase of the wave function, and the resultant diffusive dynamics such as the charge transport by real-time lattice simulation.