On the effect of exceptional points in the Liouvillian dynamics of a 1 D quantum Lorentz gas

On the effect of exceptional points in the Liouvillian dynamics of a 1D quantum Lorentz gas Kazunari Hashimoto1,∗, Kazuki Kanki2, Savannah Garmon2, Satoshi Tanaka2, and Tomio Petrosky3,4 1Graduate School of Interdisciplinary Research, University of Yamanashi, Kofu 400-8511, Japan 2Department of Physical Science, Osaka Prefecture University, Sakai 599-8531, Japan 3Institute of Industrial Science, The University of Tokyo, Tokyo 153-8505, Japan 4Center for Studies in Complex Quantum Systems, The University of Texas at Austin, TX 78712, USA ∗E-mail: kuzu.uncertainworld@gmail.com, hashimotok@yamanashi.ac.jp


Introduction
Particularly in light of studies of PT-symmetric systems initiated by C. Bender and colleagues [1,2], in recent decades the importance of non-Hermitian Hamiltonian operators has been recognized in the description of a variety of physical systems, especially certain optical systems [3][4][5][6][7].Another important arena where non-Hermitian operators play a role is in the description of irreversible processes in open quantum systems [8][9][10][11][12][13].For example, when we describe the decay process of an unstable particle, the effective Hamiltonian takes a non-Hermitian form for which the imaginary part of the eigenvalue represents the decay rate of the particle [14].The non-hermiticity of the effective Hamiltonian is associated with the environmental component of the open quantum system.Similarly, in statistical mechanics, the effective Liouvillian takes a non-Hermitian form in the thermodynamic limit, in which both intensive and extensive variables exist and the imaginary part of its eigenvalue represents the transport coefficients such as the diffusion coefficient [15][16][17].
It has been observed in many recent studies that these dynamical features of non-Hermitian operators can experience sudden changes at so-called exceptional points (EPs).For example, exceptional points are associated with the sudden appearance of complex eigenvalues in systems described by a non-Hermitian Hamiltonian [18].Their physical significance has been studied theoretically [7,[19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36] and experimentally [37][38][39][40][41]. Characteristic features appear due to an intricate structure PTEP 2016, 053A02 K. Hashimoto et al. of intersecting Riemann sheets around an EP [20].For example, a nontrivial geometrical phase is obtained when encircling the EP in the parameter space.This has been confirmed experimentally in a microwave cavity [38] and in a chaotic optical microcavity [41].Another characteristic feature is nondiagonalizability of the Hamiltonian at an EP.This leads to the Jordan block structure [29].As a result of this structure, the usual exponential time behavior of the survival amplitude of an unstable state is modified near the EP as texp[−γ t] [42].Therefore, the time evolution of the survival probability follows t2 exp[−2γ t] [29].Such a time behavior at the EP has been confirmed in a microwave cavity [39].
In our recent investigations, we have found that EPs also appear in the spectrum of the Liouvillian for various systems both in classical and quantum mechanics.In our previous paper [42], we showed that the spectrum of the Liouvillian for a 1D quantum Lorentz gas contains second-order EPs (EP2s) on the real axis of the wavenumber space.The same structure also appears in the Liouvillian spectrum of a 2D classical Lorentz gas [43] and in a 1D polaron system (Y.Sato et al., unpublished work).Although the physical effects of EPs in the Hamiltonian dynamics have been studied extensively, there are fewer works on their role in the Liouvillian dynamics.
The purpose of this paper is to study how the EPs affect the Liouvillian dynamics.More specifically, here we discuss the effect of the EPs on the time evolution of the density matrix by analyzing the time evolution of the Wigner distribution function in terms of the complex spectral analysis of the Liouville-von Neumann equation [16].The system that we study is the 1D quantum Lorentz gas [42][43][44], which is one of the simplest systems having EPs in the spectrum of the Liouvillian.In a previous paper [42], we showed that the spatial propagation is described by the telegraph equation due to the EP2s in the spectrum of the Liouvillian.In this paper, we study the details of the time evolution.
Here we compare the spatial time development of the Wigner distribution function with two different initial conditions on either side of the exceptional point.As a result, we shall show that the EPs are a threshold of two different shifting motions in space; namely, that one distribution shows a diffusive motion, while the other gives a ballistic motion.
The paper is organized as follows: In Sect.2, we introduce the system and the complex spectral representation of the Liouvillian and solve the associated eigenvalue problem.The results presented in this section have been detailed in a previous work [42], but we rehash the essential formal background here for completeness.In Sect.3, we analyze the time evolution of the Wigner distribution function relying on the spectrum presented in Sect.2, and discuss how the EPs affect the spatial time development.In Sect.4, we give concluding remarks.In the appendix, we briefly summarize some properties of the algebra in the Liouville space.

Complex spectral representation of the Liouvillian for the 1D quantum Lorentz gas
In this section, we introduce the system and the formalism of the complex spectral representation of the Liouvillian [16], summarizing the results of our previous paper [42].

System
The Lorentz gas consists of one light-mass particle (the test particle) with mass m, which interacts with N heavy particles with mass M. The Hamiltonian of the system is given by PTEP 2016, 053A02 K. Hashimoto et al. with where g is a dimensionless coupling constant, V q = V |q| is a Fourier component of the interaction potential, ≡ L/2π and q ≡ n q with q ≡ 1/ and n = 0, ±1, ±2, . ... In this paper, we consider the weak-coupling case (g 1).We also assume that V q is a continuous function at q = 0 in the limit q → 0, and that it also satisfies the condition for this limit, in order to avoid a singular transport process characteristic of the 1D system (see Eq. ( 21)).We shall also take the limit m/M → 0, in which the system is called the perfect Lorentz gas [45].
In this paper, we consider the thermodynamic limit for the heavy particles L → ∞, N → ∞ with c ≡ N /L = finite, where c is the concentration of heavy particles.In this limit we have q → 0, thus both the wavenumber and the momentum become continuous variables.At the appropriate stage, we shall replace summation with integration and the Kronecker delta δ Kr with the Dirac δ-function as −1 q → dq and δ Kr P − P → δ P − P with ≡ / , respectively.The time evolution of the density matrix for the total N + 1-particle system ρ(t) obeys the Liouville-von Neumann equation, where L H is the Liouvillian that is defined by the commutation relation with the Hamiltonian as , where L 0 is the unperturbed Liouvillian associated with H 0 and gL V is the interaction Liouvillian associated with gV .In order to discuss the "coordinate" and "momentum" dependence of the distribution of the quantum particles in parallel with classical mechanics, we introduce the Wigner distribution function, which is a quantum analog of the phase space distribution function in classical mechanics [16,46], defined by ρ W X, X j , P, P j , t ≡ 1 where {X j }, {P j }, and {k j } each represent a set of variables for the N heavy particles.The Fourier component is defined in terms of the momentum state of N + 1 particles | p, { p j } , which is an eigenstate of H 0 , as where the "wavenumbers" and the "momenta" in the Wigner representation are defined by Here the volume factor N +1 comes from the normalization in the Wigner representation (see (A8) in Appendix A).In this paper, we focus on the time evolution of the reduced density matrix for the test particle defined by where Tr hev is the partial trace over the N heavy particles.We assume that the initial condition of the system is given in the form where ρ eq hev is the Maxwell distribution of the heavy particles with temperature T , and k B is the Boltzmann constant.In the thermodynamic limit, the time evolution of the density matrix associated with the heavy particles is negligible since its deviation from ρ eq hev is proportional to 1/L.

Complex eigenvalue problem of the Liouvillian
Here we apply the formalism of the complex spectral representation of the Liouvillian [16] to the system.For this purpose, it is convenient to introduce the Liouville space representation of the operators, which is a function space spanned by operators in wave function space [16].A brief summary of the Liouville space representation is given in Appendix A.
The eigenvalue problem of the Liouvillian is represented as where the double bra-and ket-vectors stand for vectors in the Liouville space, the indices α and ν specify an eigenvalue (ν is associated with the spatial correlations (see Ref. [16])), and |F (ν) α and F(ν) α | are right-and left-eigenstates, respectively.We use the Wigner basis |k, k j ; P, P j in (66) to represent the operators as well as vectors in the Liouville space.Using the Wigner basis, the Fourier component of the Wigner distribution function ( 6) is represented in a simple form as ρ k,{k j } P, P j , t = k, k j ; P, P j |ρ(t) . ( We apply the Brillouin-Wigner-Feshbach formalism [16,17] with projection operators and Q (k) ≡ 1 − P (k) to solve the eigenvalue problem (11).For the projection operator P (k) , we have P (k) L 0 P (k) |k, {0}; P, {P j } = (k P/m) |k, {0}; P, P j and g P (k) L V P (k) = 0 because of V 0 = 0, which results from condition (3).By applying the projection operators P (k) and Q (k) on the first equation in (11), we obtain a set of equations for the P (k) and Q (k) components.By solving for the P (k) component, we find where (k) (z) is the effective Liouvillian, which is also called the collision operator in nonequilibrium statistical mechanics [16,17].The effective Liouvillian can be expanded in a power series of the coupling constant as (k) In this paper, we consider the weak-coupling case, where we take into account only the second-order term in g in the effective Liouvillian, i.e. (15).
In the following analysis, we focus our attention on the situation where the wavenumber k is much smaller than the inverse of the interaction range between the test particle and an arbitrary heavy particle denoted by d −1 , i.e., As shown below, the effective Liouvillian reduces to the collision operator in the phenomenological Boltzmann equation in this approximation.In this situation, the typical value of q appearing in V q in the third term in ( 1) is much larger than k, i.e., |k| |q|, and thus we can approximate the effective Liouvillian as where +i0 means that the effective Liouvillian 2 (z) is evaluated on the real axis of the complex z plane as approached from the upper half-plane to ensure that the time evolution is future-oriented t > 0 [16].
We focus on the dynamics of the test particle; thus, we trace out the variables for the heavy particles.We also take the limit corresponding to the perfect Lorentz gas m/M → 0. Thus we have a reduced effective Liouvillian for the test particle under the condition (16): In the thermodynamic limit, a matrix element of this operator in the Wigner representation is given by where |k; P is a reduced state for the test particle and the operator ∂ q/2 P is a displacement operator defined by ∂ q/2 P ≡ e 2 q ∂ ∂ P − e − 2 q ∂ ∂ P , where exp [a∂/∂ P] acts on a function of P as exp [a∂/∂ P] f (P) = f (P + a).Furthermore, we have assumed that the wavenumber k is negligible in the denominator in (19) as compared with q.The operator in (19) has the form ψ (k) = k P/m + ψ (0) , which is the form of the phenomenological Boltzmann operator (see, e.g., Ref. [47]).Note that the expression in (19) does not depend on the temperature of the heavy particles.This is because there is no energy transfer between the test particle and heavy particles in the limit [44].
Performing the q integration in (19), a matrix element of the effective Liouvillian ψ (k) is expressed as k; P|ψ (k)  where γ P → 0 for P → 0 (see (3)).Hence, it has nonvanishing matrix elements only between the states |k; P and |k; −P .Physically, this is because only forward and backward scattering is allowed in this 1D system.Therefore, in terms of this basis, the effective Liouvillian is represented by the following 2 × 2 non-Hermitian matrix: In terms of ψ (k) , the time evolution equation for the reduced density matrix for the test particle is given by with p(k) ≡ ∞ −∞ d P|k; P k; P|.This is equivalent to the Boltzmann equation for the perfect quantum Lorentz gas [43,48], for which the first term in the square brackets in ( 19) is called the flow term and the second term is called Boltzmann's collision term.

EPs in the spectrum of the effective Liouvillian
Let us summarize the solution of the eigenvalue problem of the effective Liouvillian (22), which was studied in Ref. [43].Its spectrum contains exceptional points in the wavenumber space, as discussed in Ref. [42].
The eigenvalues of the effective Liouvillian are given by z (k) where the wavenumber is equal to the inverse of the mean-free-length l P of the test particle with momentum P. The corresponding eigenvectors are |φ (k) Note that the eigenvectors have not yet been normalized in order to emphasize the fact that they cannot be normalized at the exceptional points k = ±k P (see (33)).In the above expressions, we have chosen the branch of the square root function as 6/15  In Fig. 1, we show the k-dependence of the real and imaginary parts of the eigenvalues for P = 0.In the figures, the dashed lines and the dot-dashed lines represent the eigenvalues z (k)B +;P and z (k)B −;P , respectively.The solid lines along the x-axis for |k/k P | < 1 in Fig. 1(A) or for |k/k P | > 1 in Fig. 1(B), respectively, are included to emphasize that these two curves connect.
The eigenstates have second-order exceptional points (EP2s) at for which the two eigenvalues and corresponding eigenvectors coalesce.Explicitly, at k = k P , we have z and, for k = −k P , Since there is only one linearly independent eigenvector at each exceptional point, the effective Liouvillian ( 22) can no longer be diagonalized.Instead, the effective Liouvillian exhibits a Jordan block structure at these points (see, e.g., Ref. [49] as well as Eq.(81) in Ref. [42]).
For k = k P , the eigenvectors are normalized in the usual way as with 7/15 They satisfy the following bi-orthonormality and bi-completeness relations:

Time evolution of the Wigner distribution function and the exceptional point
In this section, we discuss the relation of the exceptional points in the spectrum of the Liouvillian to transport processes in the system.The results presented in this section are summarized as follows: We found a shifting motion of the peak of the distribution function in space in addition to spreading as a result of a diffusion-type process.However, the mechanism of the shifting motion changes dramatically on either side of the EPs k = ±k P as follows: (1) For |k| ≤ k P , a shifting motion arises as a result of the difference between the two purely imaginary eigenvalues in this region (see Fig. 1(A)).( 2) For |k| > k P , a shifting motion is caused by the nonzero real part of the eigenvalue that leads to wave propagation with the initial velocity P/m.
To see this, we shall analyze the time evolution of the Wigner distribution function for the test particle, with By multiplying k; P| by the formal solution of the Boltzmann equation (23), i.e., we have the following expression for the Fourier transformation of the Wigner distribution function f k (P, t): The first two terms correspond to the time evolution due to the forward scattering and the third term corresponds to the time evolution due to the backward scattering.Note that the expression in ( 39) is regular at the EPs where z (k) −;P .Indeed, by taking the limit k → k P or k → −k P , we have f ±k P (P, t) = e − g 2 γ P 2 t f ±k P (P, 0) ∓ i k P P m te − g 2 γ P 2 t f ±k P (P, 0) which is identical with the expression derived in Eq. (97) in Ref. [42] in terms of the Jordan block representation.In the second and third terms in (40), we have critical damping factors t exp − g 2 γ P /2 t .Since the eigenvalues z (k) +;P and z (k) −;P change their values from pure imaginary values for |k| < k P to complex values for |k| > k P , the time evolution of the Fourier component (39) changes from an over-damping to a damped oscillation as we pass through the EPs, while critical damping appears directly at the EP itself.
In order to clearly see the relation of the spectral structure to the time evolution of the Wigner distribution function, we consider the following function as an initial condition: where χ k b (k) is a step function defined as and h(P) is a momentum distribution function that is normalized as To extract the essence of the mechanism of the shifting motion, we here assume that h(P) ≥ 0 for P ≥ 0, and h(P) = 0 for P < 0, i.e., the initial distribution is composed of particles with positive momentum.We use units in which l P = 1/k P = 1 and τ P = 1/ g 2 γ P = 1 in the numerical calculations (see, e.g., Fig. 2).With these units, the eigenvalues and the eigenvectors are independent of the value of P (see (25) and ( 27)).Therefore, the dynamics described by the Boltzmann equation is universal in the sense that the time evolution equation in these units is the same regardless of the actual value of P.

Time evolution for |k| ≤ k P
Let us first consider the situation where the initial distribution is composed of the Fourier components with k in the region |k| ≤ k P .We thus take In this case, the time evolution of the Wigner distribution function for P > 0 is expressed by where  In Fig. 2 we present the time evolution according to Eq. ( 46) in X space with a specific value of momentum P > 0, which is implicitly given as a function of l P and τ P . 1 In the figure the solid line is the initial distribution and the dashed lines are the distribution at later times.As shown in the figure, the peak of the distribution function initially shifts rightward toward X = l P in the early stage of the evolution.Afterwards, the peak no longer shifts but instead the height of the peak decreases and the width increases consistent with the diffusion equation.This can be seen by the fact that ( 46) can be approximated after the above-mentioned first stage as which is a solution of the diffusion equation with a diffusion coefficient D P given by We note that, for t → ∞, the momentum distribution function is stationary as f 0 (P, t) = f 0 (−P, t) = 1 2 f 0 (P, 0).
This implies that the momentum relaxation has been completed, and this is the physical interpretation of the fact that the peak of the distribution no longer moves.

Mechanism of the shift for the component with |k| ≤ k P
In the above, we have shown that the peak of the distribution function in X space shifts in spite of the fact that there is no real component of the spectrum for |k| ≤ k P .We now show that this is a result of the fact that the components of the distribution function corresponding to the z

Fig. 1 .
Fig. 1.Eigenvalues of the Boltzmann collision operator (25) are drawn as functions of k. (A) is the real part and (B) is the imaginary part.In each figure, the dashed lines represent an eigenvalue with α = + and the dot-dashed lines represent an eigenvalue with α = −.The solid lines are included to emphasize that these two lines are overlapping.The gray lines in (A) represent eigenvalues of the Liouvillian for a free light-mass particle y = ±(1/2)(k/k P ).

Fig. 2 .
Fig. 2. Time evolution of the Wigner distribution function for P > 0. The solid line represents the initial distribution, given by Eq. (41) with k b = k P .The initial distribution evolves as indicated by the arrows, to the distributions represented by the dotted lines.
eigenvalues(25) have asymmetry in X space, and their decay rates are different as Im z