Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                

Optical pumping through the Liouvillian skin effect

De-Huan Cai Hefei National Laboratory, University of Science and Technology of China, Hefei 230088, China    Wei Yi wyiz@ustc.edu.cn CAS Key Laboratory of Quantum Information, University of Science and Technology of China, Hefei 230026, China Anhui Province Key Laboratory of Quantum Network, University of Science and Technology of China, Hefei, 230026, China CAS Center For Excellence in Quantum Information and Quantum Physics, Hefei 230026, China Hefei National Laboratory, University of Science and Technology of China, Hefei 230088, China    Chen-Xiao Dong cxdong@hfnl.cn Hefei National Laboratory, University of Science and Technology of China, Hefei 230088, China
Abstract

The Liouvillian skin effect describes the boundary affinity of Liouvillian eignemodes that originates from the intrinsic non-Hermiticity of the Liouvillian superoperators. Dynamically, it manifests as directional flow in the transient dynamics, and the accumulation of population near open boundaries at long times. Intriguingly, similar dynamic phenomena exist in the well-known process of optical pumping, where the system is driven into a desired state (or a dark-state subspace) through the interplay of dissipation and optical drive. In this work, we show that typical optical pumping processes can indeed be understood in terms of the Liouvillian skin effect. By studying the Liouvillian spectra under different boundary conditions, we reveal that the Liouvillian spectra of the driven-dissipative pumping process sensitively depend on the boundary conditions in the state space, a signature that lies at the origin of the Liouvillian skin effect. Such a connection provides insights and practical means for designing efficient optical-pumping schemes through engineering Liouvillian gaps under the open-boundary condition. Based on these understandings, we show that the efficiency of a typical side-band cooling scheme for trapped ions can be dramatically enhanced by introducing counterintuitive dissipative channels. Our results provide a useful perspective for optical pumping, with interesting implications for state preparation and cooling.

pacs:
67.85.Lm, 03.75.Ss, 05.30.Fk

I Introduction

Optical pumping is a fundamentally important technique in the study of atomic, molecular, and optical physics [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. Originally developed to achieve population inversion necessary for lasing [2, 3], it has become the standard practice to cyclically pump atoms to a given quantum state [1, 2, 3, 15, 16], often with a well-defined magnetic quantum number. More generally, through the ingenious design of optical drive and dissipation, a quantum open system can be driven into a desired steady state (or a desired dark-state subspace) at long times [17, 18, 19, 20, 21, 22, 23]. Such general optical pumping processes are widely used for state preparation and cooling [21, 24, 25, 26, 27, 28], and offer promising paradigms for quantum simulation with atoms [17, 18, 19, 20]. Phenomenologically, a typical optical pumping process manifests two salient features, the directional flow of population in the state space, and the long-time population accumulation in the final steady state, which, given its dark-state nature, can be considered as a boundary in the state space. Intriguingly, these features also manifest in systems with the non-Hermitian skin effect, a phenomenon that has attracted extensive interest in recent years [29, 32, 33, 35, 36, 37, 38, 30, 31, 34, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55].

The non-Hermitian skin effect describes the accumulation of eigenstates near the boundaries of certain non-Hermitian systems [29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42]. It derives from the instability of the eigenvalue problems of non-Hermitian matrices to boundary perturbations, and has profound impact on the band and spectral topologies [39, 40, 42], as well as the bulk dynamics [43, 44, 45, 46, 47, 36, 37]. Experimentally, the non-Hermitian skin effect and its various manifestations have been observed in classical systems with gain and/or loss [48, 49, 50, 51, 52], and in the conditional dynamics of quantum open systems subject to post selection [53, 54, 55]. But the non-Hermitian skin effect also arises in the full-fledged quantum dynamics governed by the Lindblad master equation, wherein the Liouvillian superoperator can be mapped to a non-Hermitian matrix in an enlarged Hilbert space. Alternatively, under the master equation, the single-particle correlation evolves according to a non-Hermitian damping matrix [46, 56]. The corresponding non-Hermitian skin effect in quantum open systems, dubbed the Liouvillian skin effect [57], hosts chiral damping and directional bulk flow in the transient dynamics, as well as various boundary-sensitive long-time behaviors, such as the time scale at which the steady state is approached, and the boundary affinity of steady-state population [23, 57, 58, 59, 60, 61]. While the Liouvillian skin effect has yet to be explicitly demonstrated in experiments, the resemblance of its dynamic consequences with those of optical pumping strongly suggests an intimate, if not direct, connection between them.

In this work, we show that typical optical pumping processes can indeed be understood in terms of the Liouvillian skin effect of the underlying quantum master equation. As illustrated in Fig. 1(a), we focus on a generic optical pumping setup, where a series of otherwise independent quantum-state sectors (labeled by l𝑙litalic_l) are connected by directional dissipation. The quantum states within each sector are coupled by coherent optical fields, and may subject to additional incoherent dissipative processes in between. A discrete translational symmetry in l𝑙litalic_l is possible, but not necessary. Typical examples of such a general setup include the simplest optical pumping process in a three-level system, and the side-band cooling in trapped ions. In these examples, an open boundary condition (OBC) is naturally present, with the final state of the pumping process forming an open boundary. However, for the sake of discussion, a formal periodic boundary condition (PBC) can also be enforced by connecting the left-most and right-most sectors [as illustrated in Fig. 1(a)]. We take a typical side-band cooling configuration as an example, and study the Liouvillian spectra of the system. We find that the eigenspectra sensitively depend on the boundary conditions, a signature that lies at the origin of the Liouvillian (or non-Hermitian) skin effect. The existence of the Liouvillian skin effect is further confirmed by the directional bulk flow and the accumulation of the steady-state population at the open boundary, both of which are also natural consequences of the side-band cooling (or optical pumping) setup.

Such a connection provides insights on the further design of efficient optical pumping schemes. Specifically, since the time for the system to reach the steady state is determined by the Liouvillian gap, the efficiency of the optical pumping process can be enhanced by engineering larger Liouvillian gaps. Through analytic and numerical analyses, we identify the condition to maximize the Liouvillian gap of our system, which is surprisingly achieved by introducing dissipative processes that are opposite in direction to the bulk flow.

Our work is organized as follows. In Sec. II, we introduce the model that we consider, demonstrate the dynamic signatures of the Liouvillian skin effect, and discuss its connection with optical pumping. In Sec. III, we discuss the origin of the Liouvillian skin effect through analytic and numerical characterization of the Liouvillian spectrum. In Sec. IV, we show how the efficiency of the optical pumping process can be enhanced by optimizing the Liouvillian gap. We summarize in Sec. V.

Refer to caption
Figure 1: (a) Schematic illustration of a general optical pumping process, where a series of subsystems are connected by directional dissipation ending in a steady state. Both coherent and incoherent couplings exist between quantum states within each subsystem. (b) A typical example of the general scheme in (a), where |gket𝑔|g\rangle| italic_g ⟩ and |eket𝑒|e\rangle| italic_e ⟩ are the electronic ground and excited states, and |lket𝑙|l\rangle| italic_l ⟩ are the states in the ground- and excited-state manifolds, with energy offsets ωlsubscript𝜔𝑙\omega_{l}italic_ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT between adjacent states in the same manifold. The purple and red dashed arrows indicate the spontaneous decay from the excited state |e,lket𝑒𝑙|e,l\rangle| italic_e , italic_l ⟩ to the ground states |g,l+1ket𝑔𝑙1|g,l+1\rangle| italic_g , italic_l + 1 ⟩ and |g,lket𝑔𝑙|g,l\rangle| italic_g , italic_l ⟩, respectively. The solid arrows with different colors indicate the resonant optical couplings between states in the ground-state manifold and the excited states. The effective Rabi coupling rates ΩlsubscriptΩ𝑙\Omega_{l}roman_Ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT and the decay rates γ0,1subscript𝛾01\gamma_{0,1}italic_γ start_POSTSUBSCRIPT 0 , 1 end_POSTSUBSCRIPT are also illustrated.

II Liouvillian skin effect in optical pumping

Refer to caption
Figure 2: Liouvillian spectrum and density-matrix dynamics. (a) The red square dots and blue dots denote the eigenvalues of Liouvillian superoperator \mathcal{L}caligraphic_L in OBC and PBC, respectively. (b) and (c) are the eigenmodes of zero eigenvalue for Liouvillian superoperator \mathcal{L}caligraphic_L under OBC and PBC, respectively. (d) The Liouvillian gap as a function of the size system, for different boundary conditions. (e) and (f) are the time evolution of distribution for eigenmodes under OBC and PBC respectively, with the initial state |g,15=1ket𝑔151|g,15\rangle=1| italic_g , 15 ⟩ = 1. The green solid line indicates the time evolution of the average value for the energy state’s index, i.e., n=nnρnndelimited-⟨⟩𝑛subscript𝑛𝑛subscript𝜌𝑛𝑛\langle n\rangle=\sum_{n}n\rho_{nn}⟨ italic_n ⟩ = ∑ start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT italic_n italic_ρ start_POSTSUBSCRIPT italic_n italic_n end_POSTSUBSCRIPT. The dimension of the Hilbert space of the system is N=60𝑁60N=60italic_N = 60 except for (d). Other parameters: Ω=0.25Ω0.25\Omega=0.25roman_Ω = 0.25,γ1subscript𝛾1\gamma_{1}italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT=1.

As illustrated in Fig. 1(b), we consider a concrete example of the general optical pumping process, where external light fields couple transitions from the ground to the excited states, and, aided by dissipative processes, eventually pump the system to a given steady state. Specifically, a set of ground states with energy intervals {ωl}subscript𝜔𝑙\{\omega_{l}\}{ italic_ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT }(l=1,2𝑙12l=1,2...italic_l = 1 , 2 …) are labeled as {|g,l=|n=2l1}ket𝑔𝑙ket𝑛2𝑙1\{|g,l\rangle=|n=2l-1\rangle\}{ | italic_g , italic_l ⟩ = | italic_n = 2 italic_l - 1 ⟩ }, and the corresponding excited states are labeled as {|e,l=|n=2l}ket𝑒𝑙ket𝑛2𝑙\{|e,l\rangle=|n=2l\rangle\}{ | italic_e , italic_l ⟩ = | italic_n = 2 italic_l ⟩ }. The Rabi frequencies of the coherent optical couplings are {Ωl}subscriptΩ𝑙\{\Omega_{l}\}{ roman_Ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT }, and γ0subscript𝛾0\gamma_{0}italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT and γ1subscript𝛾1\gamma_{1}italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT are the decay rates from an excited state to different states in the ground-state manifold. Physically, l𝑙litalic_l can label magnetic quantum numbers in the ground- and excited-state hyperfine manifolds [15], in which case the scheme in Fig. 1(b) corresponds to a typical optical pumping for state preparation. Alternatively, l𝑙litalic_l can label phonon side bands in trapped ions, in which case Fig. 1(b) depicts side-band cooling [21, 24, 25, 26, 27, 28]. Regardless of the physical correspondence, the time evolution of the density matrix under the couplings of Fig. 1(b) is determined by the Lindblad master equation (we take =1Planck-constant-over-2-pi1\hbar=1roman_ℏ = 1[62, 63]

dρdt=i[H,ρ]+l,p(2Ll,pρLl,p{Ll,p,Ll,p}ρ)(ρ).d𝜌d𝑡𝑖𝐻𝜌subscript𝑙𝑝2subscript𝐿𝑙𝑝𝜌subscriptsuperscript𝐿𝑙𝑝subscriptsuperscript𝐿𝑙𝑝subscript𝐿𝑙𝑝𝜌𝜌\displaystyle\frac{\mathrm{d}\rho}{\mathrm{d}t}=-i[H,\rho]+\sum_{l,p}(2L_{l,p}% \rho L^{{\dagger}}_{l,p}-\{L^{{\dagger}}_{l,p},L_{l,p}\}\rho)\equiv\mathcal{L}% (\rho).divide start_ARG roman_d italic_ρ end_ARG start_ARG roman_d italic_t end_ARG = - italic_i [ italic_H , italic_ρ ] + ∑ start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT ( 2 italic_L start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT italic_ρ italic_L start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT - { italic_L start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT , italic_L start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT } italic_ρ ) ≡ caligraphic_L ( italic_ρ ) . (1)

Here the coherent Hamiltonian reads

H=𝐻absent\displaystyle H=italic_H = l[(j=1lωj)(|g,l+1g,l+1|+|e,l+1e,l+1|)]subscript𝑙delimited-[]superscriptsubscript𝑗1𝑙subscript𝜔𝑗ket𝑔𝑙1bra𝑔𝑙1ket𝑒𝑙1bra𝑒𝑙1\displaystyle\sum_{l}\Big{[}(\sum_{j=1}^{l}\omega_{j})(|g,l+1\rangle\langle g,% l+1|+|e,l+1\rangle\langle e,l+1|)\Big{]}∑ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT [ ( ∑ start_POSTSUBSCRIPT italic_j = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_l end_POSTSUPERSCRIPT italic_ω start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ) ( | italic_g , italic_l + 1 ⟩ ⟨ italic_g , italic_l + 1 | + | italic_e , italic_l + 1 ⟩ ⟨ italic_e , italic_l + 1 | ) ]
+lΩl(|e,lg,l+1|+H.c.),\displaystyle+\sum_{l}\Omega_{l}(|e,l\rangle\langle g,l+1|+\mathrm{H.c.}),+ ∑ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT roman_Ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT ( | italic_e , italic_l ⟩ ⟨ italic_g , italic_l + 1 | + roman_H . roman_c . ) , (2)

and the quantum jump operators are

Ll,p=0=γ0|g,l+1e,l|,Ll,p=1=γ1|g,le,l|.formulae-sequencesubscript𝐿𝑙𝑝0subscript𝛾0ket𝑔𝑙1bra𝑒𝑙subscript𝐿𝑙𝑝1subscript𝛾1ket𝑔𝑙bra𝑒𝑙\displaystyle L_{l,p=0}=\sqrt{\gamma_{0}}|g,l+1\rangle\langle e,l|,\quad L_{l,% p=1}=\sqrt{\gamma_{1}}|g,l\rangle\langle e,l|.italic_L start_POSTSUBSCRIPT italic_l , italic_p = 0 end_POSTSUBSCRIPT = square-root start_ARG italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG | italic_g , italic_l + 1 ⟩ ⟨ italic_e , italic_l | , italic_L start_POSTSUBSCRIPT italic_l , italic_p = 1 end_POSTSUBSCRIPT = square-root start_ARG italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG | italic_g , italic_l ⟩ ⟨ italic_e , italic_l | . (3)

We denote the Hilbert-space dimension of the system as N𝑁Nitalic_N, with nmax=2lmax=Nsubscript𝑛𝑚𝑎𝑥2subscript𝑙𝑚𝑎𝑥𝑁n_{max}=2l_{max}=Nitalic_n start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT = 2 italic_l start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT = italic_N. Then the right and left eigenmodes of the Liouvillian superoperator \mathcal{L}caligraphic_L, defined in an N2superscript𝑁2N^{2}italic_N start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT-dimensional extended Hilbert space, are given by

(ρμR)=λμρμR,(ρμL)=λμρμL,formulae-sequencesuperscriptsubscript𝜌𝜇𝑅subscript𝜆𝜇superscriptsubscript𝜌𝜇𝑅superscriptsuperscriptsubscript𝜌𝜇𝐿superscriptsubscript𝜆𝜇superscriptsubscript𝜌𝜇𝐿\displaystyle\mathcal{L}(\rho_{\mu}^{R})=\lambda_{\mu}\rho_{\mu}^{R},\quad% \quad\mathcal{L}^{{\dagger}}(\rho_{\mu}^{L})=\lambda_{\mu}^{\ast}\rho_{\mu}^{L},caligraphic_L ( italic_ρ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT ) = italic_λ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT , caligraphic_L start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ( italic_ρ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_L end_POSTSUPERSCRIPT ) = italic_λ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT italic_ρ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_L end_POSTSUPERSCRIPT , (4)

with μ=1,2,3,,N2𝜇123superscript𝑁2\mu=1,2,3,...,N^{2}italic_μ = 1 , 2 , 3 , … , italic_N start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. The right and left eigenmodes are normalized as ρμR|ρμR=ρμL|ρμL=1inner-productsubscriptsuperscript𝜌𝑅𝜇subscriptsuperscript𝜌𝑅𝜇inner-productsubscriptsuperscript𝜌𝐿𝜇subscriptsuperscript𝜌𝐿𝜇1\sqrt{\langle\rho^{R}_{\mu}|\rho^{R}_{\mu}\rangle}=\sqrt{\langle\rho^{L}_{\mu}% |\rho^{L}_{\mu}\rangle}=1square-root start_ARG ⟨ italic_ρ start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT | italic_ρ start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT ⟩ end_ARG = square-root start_ARG ⟨ italic_ρ start_POSTSUPERSCRIPT italic_L end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT | italic_ρ start_POSTSUPERSCRIPT italic_L end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT ⟩ end_ARG = 1, and are orthogonal to each other (ρμL|ρνR=0inner-productsubscriptsuperscript𝜌𝐿𝜇subscriptsuperscript𝜌𝑅𝜈0\sqrt{\langle\rho^{L}_{\mu}|\rho^{R}_{\nu}\rangle}=0square-root start_ARG ⟨ italic_ρ start_POSTSUPERSCRIPT italic_L end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT | italic_ρ start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT ⟩ end_ARG = 0) when their eigenvalues are different (λμλνsubscript𝜆𝜇subscript𝜆𝜈\lambda_{\mu}\neq\lambda_{\nu}italic_λ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT ≠ italic_λ start_POSTSUBSCRIPT italic_ν end_POSTSUBSCRIPT). In particular, the eigenmodes of \mathcal{L}caligraphic_L with vanishing eigenvalues are the steady states of the system, with (ρss)=0subscript𝜌𝑠𝑠0\mathcal{L}(\rho_{ss})=0caligraphic_L ( italic_ρ start_POSTSUBSCRIPT italic_s italic_s end_POSTSUBSCRIPT ) = 0.

It follows that the density matrix of the initial state can be expanded as

ρini=μ=1N2cμρμR,subscript𝜌inisuperscriptsubscript𝜇1superscript𝑁2subscript𝑐𝜇superscriptsubscript𝜌𝜇𝑅\displaystyle\rho_{\mathrm{ini}}=\sum_{\mu=1}^{N^{2}}c_{\mu}\rho_{\mu}^{R},italic_ρ start_POSTSUBSCRIPT roman_ini end_POSTSUBSCRIPT = ∑ start_POSTSUBSCRIPT italic_μ = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT italic_c start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT , (5)

where cμ=ρμL|ρini/ρμL|ρμRsubscript𝑐𝜇inner-productsubscriptsuperscript𝜌𝐿𝜇subscript𝜌iniinner-productsubscriptsuperscript𝜌𝐿𝜇subscriptsuperscript𝜌𝑅𝜇c_{\mu}=\langle\rho^{L}_{\mu}|\rho_{\mathrm{ini}}\rangle/\langle\rho^{L}_{\mu}% |\rho^{R}_{\mu}\rangleitalic_c start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT = ⟨ italic_ρ start_POSTSUPERSCRIPT italic_L end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT | italic_ρ start_POSTSUBSCRIPT roman_ini end_POSTSUBSCRIPT ⟩ / ⟨ italic_ρ start_POSTSUPERSCRIPT italic_L end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT | italic_ρ start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT ⟩ according to the completeness condition μ|ρμRρμL|/ρμL|ρμR=1subscript𝜇ketsubscriptsuperscript𝜌𝑅𝜇brasubscriptsuperscript𝜌𝐿𝜇inner-productsubscriptsuperscript𝜌𝐿𝜇subscriptsuperscript𝜌𝑅𝜇1\sum_{\mu}|\rho^{R}_{\mu}\rangle\langle\rho^{L}_{\mu}|/\langle\rho^{L}_{\mu}|% \rho^{R}_{\mu}\rangle=1∑ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT | italic_ρ start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT ⟩ ⟨ italic_ρ start_POSTSUPERSCRIPT italic_L end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT | / ⟨ italic_ρ start_POSTSUPERSCRIPT italic_L end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT | italic_ρ start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT ⟩ = 1. Thus, the time evolution of the density matrix can be written as

ρ(t)=μ=1N2cμeλμtρμR.𝜌𝑡superscriptsubscript𝜇1superscript𝑁2subscript𝑐𝜇superscript𝑒subscript𝜆𝜇𝑡superscriptsubscript𝜌𝜇𝑅\displaystyle\rho(t)=\sum_{\mu=1}^{N^{2}}c_{\mu}e^{\lambda_{\mu}t}\rho_{\mu}^{% R}.italic_ρ ( italic_t ) = ∑ start_POSTSUBSCRIPT italic_μ = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT italic_c start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_e start_POSTSUPERSCRIPT italic_λ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_t end_POSTSUPERSCRIPT italic_ρ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT . (6)

Note that the real parts of the eigenvalues of the excited eigenmodes (those that are not steady states) must be negative to ensure that their contributions in Eq. 6 would be exponentially small after a long enough time evolution, as the system approaches the steady states. Here we set ρss=ρμ=1Rsubscript𝜌𝑠𝑠superscriptsubscript𝜌𝜇1𝑅\rho_{ss}=\rho_{\mu=1}^{R}italic_ρ start_POSTSUBSCRIPT italic_s italic_s end_POSTSUBSCRIPT = italic_ρ start_POSTSUBSCRIPT italic_μ = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT, and assume that all eigenvalues are indexed in descending order according to their real parts: 0=λ1>Re[λ2]Re[λ3]Re[λN2]0subscript𝜆1Redelimited-[]subscript𝜆2Redelimited-[]subscript𝜆3Redelimited-[]subscript𝜆superscript𝑁20=\lambda_{1}>\mathrm{Re}[\lambda_{2}]\geq\mathrm{Re}[\lambda_{3}]...\geq% \mathrm{Re}[\lambda_{N^{2}}]0 = italic_λ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT > roman_Re [ italic_λ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ] ≥ roman_Re [ italic_λ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT ] … ≥ roman_Re [ italic_λ start_POSTSUBSCRIPT italic_N start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_POSTSUBSCRIPT ]. Equation 6 can then be rewritten as

ρ(t)=ρss+μ=2N2cμeλμtρμR.𝜌𝑡subscript𝜌𝑠𝑠superscriptsubscript𝜇2superscript𝑁2subscript𝑐𝜇superscript𝑒subscript𝜆𝜇𝑡superscriptsubscript𝜌𝜇𝑅\displaystyle\rho(t)=\rho_{ss}+\sum_{\mu=2}^{N^{2}}c_{\mu}e^{\lambda_{\mu}t}% \rho_{\mu}^{R}.italic_ρ ( italic_t ) = italic_ρ start_POSTSUBSCRIPT italic_s italic_s end_POSTSUBSCRIPT + ∑ start_POSTSUBSCRIPT italic_μ = 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT italic_c start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_e start_POSTSUPERSCRIPT italic_λ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT italic_t end_POSTSUPERSCRIPT italic_ρ start_POSTSUBSCRIPT italic_μ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT . (7)

Importantly, the Liouvillian gap is defined as Δ=|Re[λ2]|ΔRedelimited-[]subscript𝜆2\Delta=|\mathrm{Re}[\lambda_{2}]|roman_Δ = | roman_Re [ italic_λ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ] |, which describes the asymptotic decay rate of the system toward the steady states at long times [64].

We first consider the simple case with ωl=0subscript𝜔𝑙0\omega_{l}=0italic_ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = 0, Ωl=ΩsubscriptΩ𝑙Ω\Omega_{l}=\Omegaroman_Ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = roman_Ω, and γ0=0subscript𝛾00\gamma_{0}=0italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0. It follows that Hamiltonian (2) is simplified to H=lΩ(|e,lg,l+1|+H.c.)H=\sum_{l}\Omega(|e,l\rangle\langle g,l+1|+\mathrm{H.c.})italic_H = ∑ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT roman_Ω ( | italic_e , italic_l ⟩ ⟨ italic_g , italic_l + 1 | + roman_H . roman_c . ), and only a single quantum jump process exists for each pair of ground and excited states, given by Ll,1subscript𝐿𝑙1L_{l,1}italic_L start_POSTSUBSCRIPT italic_l , 1 end_POSTSUBSCRIPT. In Hamiltonian (2), states with the smallest and largest n𝑛nitalic_n indices are not coupled. This corresponds to an OBC in the state space. By contrast, one may consider an artificial PBC, where all states are cyclically coupled. Such a PBC is achieved by adding the term (Ω|e,lmaxg,1|+H.c.)(\Omega|e,l_{\mathrm{max}}\rangle\langle g,1|+\mathrm{H.c.})( roman_Ω | italic_e , italic_l start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ⟩ ⟨ italic_g , 1 | + roman_H . roman_c . ) to Eq. (2), where lmaxsubscript𝑙maxl_{\mathrm{max}}italic_l start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT is the maximum l𝑙litalic_l. Although the PBC is unphysical, it offers insights to the setup as we detail below. Alternatively, one may consider the state label n𝑛nitalic_n as lattice sites along a synthetic dimension. Different boundary conditions in the synthetic dimension then directly correspond to boundary conditions in the state space. With these understandings, we now study the Liouvillian spectrum and dynamic evolution of the master equation, under different boundary conditions.

As depicted in Fig. 2(a), the eigenvalues of the Liouvillian superoperator \mathcal{L}caligraphic_L under the PBC form a closed loop on the complex plane, enclosing those under the OBC. This is reminiscent of the spectral topology of non-Hermitian Hamiltonians with the skin effect, and is an outstanding signature for the Liouvillian skin effect. In either case, the drastic difference in the eigenspectrum under different boundary conditions originates from the instability of non-Hermitian matrices to boundary perturbations. Fig.  2(b) shows the density-matrix elements ρnmsubscript𝜌𝑛𝑚\rho_{nm}italic_ρ start_POSTSUBSCRIPT italic_n italic_m end_POSTSUBSCRIPT of the steady state under the OBC. Here the density-matrix element is defined as ρnm=n|ρμ=1R|msubscript𝜌𝑛𝑚quantum-operator-product𝑛superscriptsubscript𝜌𝜇1𝑅𝑚\rho_{nm}=\langle n|\rho_{\mu=1}^{R}|m\rangleitalic_ρ start_POSTSUBSCRIPT italic_n italic_m end_POSTSUBSCRIPT = ⟨ italic_n | italic_ρ start_POSTSUBSCRIPT italic_μ = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_R end_POSTSUPERSCRIPT | italic_m ⟩. The steady state is indeed localized in |g,l=1ket𝑔𝑙1|g,l=1\rangle| italic_g , italic_l = 1 ⟩, corresponding to an open boundary. The corresponding steady state under the PBC is shown in Fig. 2(c), where uniform distributions in l𝑙litalic_l are observed for both the ground and excited states. A closer look reveals that, in the steady state under the PBC, the majority of the population is in the ground state.

Another drastic distinction between the Liouvillian spectrum under OBC and PBC is the Liouvillian gap. As shown in Fig. 2(d), the Liouvillian gap ΔΔ\Deltaroman_Δ tends to zero as the size of the system increases under the PBC. By contrast, the gap is independent of the system size under the OBC. A finite Liouvillian gap implies that the density matrix in Eq. (7) converges exponentially fast to the steady state at long times. Whereas a vanishing Liouvillian gap implies an algebraic convergence, such that the relaxation time diverges for Δ0Δ0\Delta\rightarrow 0roman_Δ → 0 [65].

Taking the size of the system as N=60𝑁60N=60italic_N = 60 in Fig. 2(e) and (f), we evolve the system according to Eq. (1), while setting the initial state to |g,15=|n=29ket𝑔15ket𝑛29|g,15\rangle=|n=29\rangle| italic_g , 15 ⟩ = | italic_n = 29 ⟩. Under the OBC, the occupation rapidly flows toward the boundary and eventually evolves to the steady state as shown in Fig. 2(b). This is the dynamic signature of the Liouvillian skin effect. In the context of optical pumping, such a directional flow is the underlying mechanism for state preparation and cooling. For instance, in trapped ions, the index l𝑙litalic_l corresponds to the phonon modes. The coherent optical drives are implemented by side-band couplings, and the directional flow toward l=0𝑙0l=0italic_l = 0 corresponds to cooling of the external ion motion. The timescale or efficiency of the cooling process is then determined by the Liouvillian gap under the OBC. Under the PBC, since the Liouvillian gap is much smaller, the time it takes to relax to the steady state is much longer, and diverges in the thermodynamic limit.

Refer to caption
Figure 3: The Liouvillian spectrum under the OBC (red square) and PBC (blue dots), respectively. (a) We take Ωl=lΩ1subscriptΩ𝑙𝑙subscriptΩ1\Omega_{l}=\sqrt{l}\Omega_{1}roman_Ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = square-root start_ARG italic_l end_ARG roman_Ω start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT, ω=0𝜔0\omega=0italic_ω = 0. (b) We take Ωl=lΩ1subscriptΩ𝑙𝑙subscriptΩ1\Omega_{l}=\sqrt{l}\Omega_{1}roman_Ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = square-root start_ARG italic_l end_ARG roman_Ω start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT, ω=0.2𝜔0.2\omega=0.2italic_ω = 0.2. Other parameters are Ω1=0.25subscriptΩ10.25\Omega_{1}=0.25roman_Ω start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = 0.25, γ1=1subscript𝛾11\gamma_{1}=1italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = 1, γ0=0subscript𝛾00\gamma_{0}=0italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0.

More generally, we have ωl0subscript𝜔𝑙0\omega_{l}\neq 0italic_ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT ≠ 0, and state-dependent ΩlsubscriptΩ𝑙\Omega_{l}roman_Ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT (but still with γ0=0subscript𝛾00\gamma_{0}=0italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0). The Liouvillian spectrum under the PBC no longer encloses the one under OBC, but they remain different, as shown in Fig. 3. The Liouvillian skin effect persists, and the steady state under the OBC remains the same as that in Fig. 2(b). The long-time evolution of the system generates a directional flow, similar to the results shown in Fig. 2(e), and the relaxation time depends on the Liouvillian gap. In the next section, we will illustrate the structure of the Liouvillian spectrum and the origin of the Liouvillian skin effect for the general case through analytic methods.

Refer to caption
Figure 4: The matrix structure of the Liouvillian superoperator. (a) shows the block-diagonlized structure of the Liouvillian operator; (b) shows the block upper-triangular structure of 0subscript0\mathcal{L}_{0}caligraphic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT with OBC; (c) shows the block-circulant structure of 0subscript0\mathcal{L}_{0}caligraphic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT with PBC. Here blocks with different colors represent different matrix elements, and the orange blocks indicate the matrix elements from the recycling terms l,p2Ll,pρLl,psubscript𝑙𝑝2subscript𝐿𝑙𝑝𝜌subscriptsuperscript𝐿𝑙𝑝\sum_{l,p}2L_{l,p}\rho L^{{\dagger}}_{l,p}∑ start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT 2 italic_L start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT italic_ρ italic_L start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT. We fix the system size to N=8𝑁8N=8italic_N = 8.

III Analytic study of the Liouvillian

Refer to caption
Figure 5: Liouvillian gap and the long-time damping dynamic. (a) The long-time damping follows an exponential law after an initial power law stage, where the system size is N=50𝑁50N=50italic_N = 50. (b) The Liouvillian gaps are independent of the dimension of the system in OBC, where different Rabi frequencies and energy intervals result in different gaps and diffrent dampings in (a). (c) The Liouvillian gap as a function of energy interval ω𝜔\omegaitalic_ω and Rabi frequency ΩΩ\Omegaroman_Ω with N=50𝑁50N=50italic_N = 50. We perform numerical calculations with OBC, the decay rates are γ0=0subscript𝛾00\gamma_{0}=0italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0 and γ1=1subscript𝛾11\gamma_{1}=1italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = 1.

In this section, we analytically solve the spectrum of the Liouvillian superoperator to elucidate the origin of the Liouvillian skin effect described in the previous section.

First, we rearrange the Lindblad equation (1) into

dρdt=i[Heff,ρ]+l,p2Ll,pρLl,p,d𝜌d𝑡𝑖subscript𝐻eff𝜌subscript𝑙𝑝2subscript𝐿𝑙𝑝𝜌subscriptsuperscript𝐿𝑙𝑝\frac{\mathrm{d}\rho}{\mathrm{d}t}=-i[H_{\text{eff}},\rho]+\sum_{l,p}2L_{l,p}% \rho L^{{\dagger}}_{l,p},divide start_ARG roman_d italic_ρ end_ARG start_ARG roman_d italic_t end_ARG = - italic_i [ italic_H start_POSTSUBSCRIPT eff end_POSTSUBSCRIPT , italic_ρ ] + ∑ start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT 2 italic_L start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT italic_ρ italic_L start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT , (8)

where the effective non-Hermitian Hamilton is

Heff=Hil,pLl,pLl,p.subscript𝐻eff𝐻𝑖subscript𝑙𝑝subscriptsuperscript𝐿𝑙𝑝subscript𝐿𝑙𝑝H_{\mathrm{eff}}=H-i\sum_{l,p}L^{{\dagger}}_{l,p}L_{l,p}.italic_H start_POSTSUBSCRIPT roman_eff end_POSTSUBSCRIPT = italic_H - italic_i ∑ start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT italic_L start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT italic_L start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT . (9)

We observe that the effective non-Hermitian Hamiltonian for the setup in Fig. 1 is block-diagonal. This is because both the coherent Hamiltonian H𝐻Hitalic_H and the terms il,pLl,pLl,p𝑖subscript𝑙𝑝subscriptsuperscript𝐿𝑙𝑝subscript𝐿𝑙𝑝-i\sum_{l,p}L^{{\dagger}}_{l,p}L_{l,p}- italic_i ∑ start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT italic_L start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT italic_L start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT are block-diagonal with respect to the subsystems shown in Fig. 1(a). We hence denote

Heff=H1H2Hm,subscript𝐻effdirect-sumsubscript𝐻1subscript𝐻2subscript𝐻𝑚H_{\text{eff}}=H_{1}\oplus H_{2}\oplus\cdots\oplus H_{m},italic_H start_POSTSUBSCRIPT eff end_POSTSUBSCRIPT = italic_H start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ⊕ italic_H start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ⊕ ⋯ ⊕ italic_H start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT , (10)

where m𝑚mitalic_m represents the number of subsystems and each Hisubscript𝐻𝑖H_{i}italic_H start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT represents an individual subsystem with the Hilbert-space dimension nisubscript𝑛𝑖n_{i}italic_n start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT, with j=1mnj=Nsuperscriptsubscript𝑗1𝑚subscript𝑛𝑗𝑁\sum_{j=1}^{m}n_{j}=N∑ start_POSTSUBSCRIPT italic_j = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT italic_n start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT = italic_N.

For our model in Fig. 1(b), we find that the effective Hamilton is composed of two single-level systems with on-site energies 00 and l=1N/21ωliγ1superscriptsubscript𝑙1𝑁21subscript𝜔𝑙𝑖subscript𝛾1\sum_{l=1}^{N/2-1}\omega_{l}-i\gamma_{1}∑ start_POSTSUBSCRIPT italic_l = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N / 2 - 1 end_POSTSUPERSCRIPT italic_ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT - italic_i italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT, and a series of (N/21)𝑁21(N/2-1)( italic_N / 2 - 1 ) two-level subsystems each described by the Hamiltonian

Hj=(l=1jωl)|g,j+1g,j+1|+(l=1j1ωliγ1)|e,je,j|+Ωj(|e,jg,j+1|+|g,j+1e,j|),subscript𝐻𝑗superscriptsubscript𝑙1𝑗subscript𝜔𝑙ket𝑔𝑗1quantum-operator-product𝑔𝑗1superscriptsubscript𝑙1𝑗1subscript𝜔𝑙𝑖subscript𝛾1𝑒𝑗bra𝑒𝑗subscriptΩ𝑗ket𝑒𝑗bra𝑔𝑗1ket𝑔𝑗1bra𝑒𝑗\begin{split}H_{j}=&(\sum_{l=1}^{j}\omega_{l})|g,j+1\rangle\langle g,j+1|\\ +&(\sum_{l=1}^{j-1}\omega_{l}-i\gamma_{1})|e,j\rangle\langle e,j|\\ +&\Omega_{j}(|e,j\rangle\langle g,j+1|+|g,j+1\rangle\langle e,j|),\end{split}start_ROW start_CELL italic_H start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT = end_CELL start_CELL ( ∑ start_POSTSUBSCRIPT italic_l = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_j end_POSTSUPERSCRIPT italic_ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT ) | italic_g , italic_j + 1 ⟩ ⟨ italic_g , italic_j + 1 | end_CELL end_ROW start_ROW start_CELL + end_CELL start_CELL ( ∑ start_POSTSUBSCRIPT italic_l = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_j - 1 end_POSTSUPERSCRIPT italic_ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT - italic_i italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ) | italic_e , italic_j ⟩ ⟨ italic_e , italic_j | end_CELL end_ROW start_ROW start_CELL + end_CELL start_CELL roman_Ω start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ( | italic_e , italic_j ⟩ ⟨ italic_g , italic_j + 1 | + | italic_g , italic_j + 1 ⟩ ⟨ italic_e , italic_j | ) , end_CELL end_ROW (11)

where j=1,2,,N/212𝑁21,2,\cdots,N/21 , 2 , ⋯ , italic_N / 2.

Under the PBC, we observe that all N/2𝑁2N/2italic_N / 2 subsystems in the effective Hamilton are two-level systems, given by the Hamiltonian in Eq. (11), but with

HN/2=(l=1N/21ωliγ1)|e,N/2e,N/2|+ΩN/2(|e,N/2g,1|+|g,1e,N/2|).subscript𝐻𝑁2superscriptsubscript𝑙1𝑁21subscript𝜔𝑙𝑖subscript𝛾1ket𝑒𝑁2bra𝑒𝑁2subscriptΩ𝑁2ket𝑒𝑁2bra𝑔1ket𝑔1bra𝑒𝑁2\begin{split}H_{N/2}=&(\sum_{l=1}^{N/2-1}\omega_{l}-i\gamma_{1})|e,N/2\rangle% \langle e,N/2|\\ +&\Omega_{N/2}(|e,N/2\rangle\langle g,1|+|g,1\rangle\langle e,N/2|).\end{split}start_ROW start_CELL italic_H start_POSTSUBSCRIPT italic_N / 2 end_POSTSUBSCRIPT = end_CELL start_CELL ( ∑ start_POSTSUBSCRIPT italic_l = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N / 2 - 1 end_POSTSUPERSCRIPT italic_ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT - italic_i italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ) | italic_e , italic_N / 2 ⟩ ⟨ italic_e , italic_N / 2 | end_CELL end_ROW start_ROW start_CELL + end_CELL start_CELL roman_Ω start_POSTSUBSCRIPT italic_N / 2 end_POSTSUBSCRIPT ( | italic_e , italic_N / 2 ⟩ ⟨ italic_g , 1 | + | italic_g , 1 ⟩ ⟨ italic_e , italic_N / 2 | ) . end_CELL end_ROW (12)

Additionally, we observe that the contribution from the recycling terms l,p2Ll,pρLl,psubscript𝑙𝑝2subscript𝐿𝑙𝑝𝜌subscriptsuperscript𝐿𝑙𝑝\sum_{l,p}2L_{l,p}\rho L^{{\dagger}}_{l,p}∑ start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT 2 italic_L start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT italic_ρ italic_L start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT exists either between two adjacent subsystems, or within an individual subsystem (defined as 0subscript0\mathcal{L}_{0}caligraphic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT below). Hence the overall Liouvillian superoperator is also block-diagonal in its matrix form, as illustrated in Fig 4(a). The large block with intra-block recycling-term contribution is given by the Liouvillian

0=ij=1m(HjnjnjHj)+l,p2Ll,pρLl,p,subscript0𝑖subscriptsuperscript𝑚𝑗1tensor-productsubscript𝐻𝑗subscriptsubscript𝑛𝑗tensor-productsubscriptsubscript𝑛𝑗subscript𝐻𝑗subscript𝑙𝑝2subscript𝐿𝑙𝑝𝜌subscriptsuperscript𝐿𝑙𝑝\mathcal{L}_{0}=-i\sum^{m}_{j=1}(H_{j}\otimes\mathcal{I}_{n_{j}}-\mathcal{I}_{% n_{j}}\otimes H_{j})+\sum_{l,p}2L_{l,p}\rho L^{{\dagger}}_{l,p},caligraphic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = - italic_i ∑ start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_j = 1 end_POSTSUBSCRIPT ( italic_H start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ⊗ caligraphic_I start_POSTSUBSCRIPT italic_n start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_POSTSUBSCRIPT - caligraphic_I start_POSTSUBSCRIPT italic_n start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_POSTSUBSCRIPT ⊗ italic_H start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ) + ∑ start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT 2 italic_L start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT italic_ρ italic_L start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_l , italic_p end_POSTSUBSCRIPT , (13)

with the dimension j=1mnj2subscriptsuperscript𝑚𝑗1superscriptsubscript𝑛𝑗2\sum^{m}_{j=1}n_{j}^{2}∑ start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_j = 1 end_POSTSUBSCRIPT italic_n start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. Other blocks are given by

lj=i(HjnlnjHl),subscript𝑙𝑗𝑖tensor-productsubscript𝐻𝑗subscriptsubscript𝑛𝑙tensor-productsubscriptsubscript𝑛𝑗subscript𝐻𝑙\mathcal{L}_{lj}=-i(H_{j}\otimes\mathcal{I}_{n_{l}}-\mathcal{I}_{n_{j}}\otimes H% _{l}),caligraphic_L start_POSTSUBSCRIPT italic_l italic_j end_POSTSUBSCRIPT = - italic_i ( italic_H start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ⊗ caligraphic_I start_POSTSUBSCRIPT italic_n start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_POSTSUBSCRIPT - caligraphic_I start_POSTSUBSCRIPT italic_n start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_POSTSUBSCRIPT ⊗ italic_H start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT ) , (14)

with dimensions nlnjsubscript𝑛𝑙subscript𝑛𝑗n_{l}n_{j}italic_n start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT italic_n start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT, where l,j=1,2,,N/21formulae-sequence𝑙𝑗12𝑁21l,j=1,2,\cdots,N/2-1italic_l , italic_j = 1 , 2 , ⋯ , italic_N / 2 - 1 and lj𝑙𝑗l\neq jitalic_l ≠ italic_j. Here nsubscript𝑛\mathcal{I}_{n}caligraphic_I start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT is the identity matrix with dimension n𝑛nitalic_n.

Due to the block-diagonal structure of the Liouvillian, its eigenspectrum is analytically solvable by diagonalizing 0subscript0\mathcal{L}_{0}caligraphic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT and ljsubscript𝑙𝑗\mathcal{L}_{lj}caligraphic_L start_POSTSUBSCRIPT italic_l italic_j end_POSTSUBSCRIPT, respectively. Specifically, in our model, since the dimensions of Hjsubscript𝐻𝑗H_{j}italic_H start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT are less than or equal to 2, the dimension of any given ljsubscript𝑙𝑗\mathcal{L}_{lj}caligraphic_L start_POSTSUBSCRIPT italic_l italic_j end_POSTSUBSCRIPT is less than or equal to 4. And 0subscript0\mathcal{L}_{0}caligraphic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT is a special matrix that is easy to diagonalize.

We first study the case with OBC. In this case, we observe that the dissipation between two adjacent subsystems is directional, which makes 0subscript0\mathcal{L}_{0}caligraphic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT a block upper-triangular matrix, as illustrated in Fig 4(b). The eigenspectrum of 0subscript0\mathcal{L}_{0}caligraphic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT is then the union of the spectra of the diagonal blocks. In the presence of translational symmetry with ωl=0subscript𝜔𝑙0\omega_{l}=0italic_ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = 0 and Ωl=ΩsubscriptΩ𝑙Ω\Omega_{l}=\Omegaroman_Ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = roman_Ω, diagonal blocks of 0subscript0\mathcal{L}_{0}caligraphic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT are invariant with increasing system size. Consequently, the Liouvillian gap remains constant as the system size changes, consistent with discussions in the previous section. Furthermore, due to the block upper-triangular structure of 0subscript0\mathcal{L}_{0}caligraphic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, some eigenvectors from 0subscript0\mathcal{L}_{0}caligraphic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT are localized within the subsystems near the boundaries of the entire state space. As we detail in the Appendix, such a localization persists even as the translational symmetry is broken (for general values of ωlsubscript𝜔𝑙\omega_{l}italic_ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT and/or ΩlsubscriptΩ𝑙\Omega_{l}roman_Ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT).

Under the PBC, when ωl=0subscript𝜔𝑙0\omega_{l}=0italic_ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = 0 and Ωl=ΩsubscriptΩ𝑙Ω\Omega_{l}=\Omegaroman_Ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = roman_Ω, due to the translational symmetry, 0subscript0\mathcal{L}_{0}caligraphic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT forms a block-circulant matrix, illustrated in Fig 4(c). We can thus visualize it as a four-band non-Hermitian one-dimensional lattice model with PBC. The eigenspectrum is analytically solvable, and we find that the Liouvillian gap approaches zero when the system size tends to infinity. (More details are shown in the Appendix).

Therefore, the Liouvillian skin effect observed in the previous section mathematically originates from the difference in 0subscript0\mathcal{L}_{0}caligraphic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT under different boundary conditions. Physically, the Liouvillian skin effect in our system arises from the divisibility of the effective Hamilton and the non-reciprocal recycling terms. This phenomenon is analogous to the non-Hermitian skin effect observed in non-Hermitian lattice models. Finally, we remark that our discussions here can be generalized to generic optical pumping setups illustrated in Fig. 1(a).

IV Designing efficient pumping scheme

In this section, we show that the pumping scheme in Fig. 1(b) can be optimized based on the understandings above. Here we set ωl=ωsubscript𝜔𝑙𝜔\omega_{l}=\omegaitalic_ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = italic_ω and Ωl=ΩsubscriptΩ𝑙Ω\Omega_{l}=\Omegaroman_Ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = roman_Ω to simplify discussions, but our results qualitatively hold for schemes without the translational symmetry. The latter can be important for side-band cooling in trapped ions in the Lamb-Dicke regime, where the coupling strength between different side bands scale as n𝑛\sqrt{n}square-root start_ARG italic_n end_ARG [22].

In our system, any initial state evolves towards a steady state. To quantify the damping dynamics, we calculate the particle-number deviation from that of the steady state, defined as n~(t)=Tr[ρ(t)ρ(t)]~𝑛𝑡Trdelimited-[]𝜌𝑡𝜌𝑡\tilde{n}(t)=\text{Tr}[\rho(t)-\rho(t\to\infty)]over~ start_ARG italic_n end_ARG ( italic_t ) = Tr [ italic_ρ ( italic_t ) - italic_ρ ( italic_t → ∞ ) ]. As shown in Fig. 5(a), the damping of n~(t)~𝑛𝑡\tilde{n}(t)over~ start_ARG italic_n end_ARG ( italic_t ) depends on the initial state and the Liouvillian gap. With the same initial states, the damping dynamic accelerates when the Liouvillian gap increases.

Next, we explore the relationship between the system parameters and the Liouvillian gap. As discussed earlier, when γ0=0subscript𝛾00\gamma_{0}=0italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0, the spectrum of our system is independent of the system size, as illustrated in Fig. 5(b). Generally, in the experiments, the energy offset ω𝜔\omegaitalic_ω is smaller than the Rabi frequency ΩΩ\Omegaroman_Ω. Specifically, when ω=0𝜔0\omega=0italic_ω = 0, the Liouvillian gap increases with Ω/γ1Ωsubscript𝛾1\Omega/\gamma_{1}roman_Ω / italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT when Ω/γ1<1/4Ωsubscript𝛾114\Omega/\gamma_{1}<1/4roman_Ω / italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT < 1 / 4, reaching a maximum of γ1/4subscript𝛾14\gamma_{1}/4italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / 4 when Ω/γ1>1/4Ωsubscript𝛾114\Omega/\gamma_{1}>1/4roman_Ω / italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT > 1 / 4, as illustrated in Fig. 5(c). Subsequently, if ω0𝜔0\omega\neq 0italic_ω ≠ 0, the Liouvillian gap consistently decreases with increasing ω𝜔\omegaitalic_ω. As a result, the maximum possible Liouvillian gap is γ1/4subscript𝛾14\gamma_{1}/4italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / 4 when γ0=0subscript𝛾00\gamma_{0}=0italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0, which is consistent with previous studies [21]. In the following, we aim to further increase the Liouvillian gap by introducing new decay channels.

Refer to caption
Figure 6: Liouvillian gap for γ00subscript𝛾00\gamma_{0}\neq 0italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ≠ 0. (a) and (b) show the Liouvillian gap with increasing γ0subscript𝛾0\gamma_{0}italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT for different ΩΩ\Omegaroman_Ω, for ω=0𝜔0\omega=0italic_ω = 0 and 0.30.30.30.3, respectively. (c) The change rate of the Liouvillian gap at γ0subscript𝛾0\gamma_{0}italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT (dΔ/dγ0|γ0=0evaluated-atdΔdsubscript𝛾0subscript𝛾00\left.\mathrm{d}\Delta/\mathrm{d}\gamma_{0}\right|_{\gamma_{0}=0}roman_d roman_Δ / roman_d italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT | start_POSTSUBSCRIPT italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0 end_POSTSUBSCRIPT). (d) and (e) show the maximum Liouvillian gap and the optimal decay rate γ0,maxsubscript𝛾0max\gamma_{0,\text{max}}italic_γ start_POSTSUBSCRIPT 0 , max end_POSTSUBSCRIPT for the maximum Liouvillian gap, under different ω𝜔\omegaitalic_ω and ΩΩ\Omegaroman_Ω. The magenta, red, green and blue lines (include solid and dashed lines) in (d) and (e) correspond to the parameters of ω=0𝜔0\omega=0italic_ω = 0, ω=0.3𝜔0.3\omega=0.3italic_ω = 0.3, ω=1𝜔1\omega=1italic_ω = 1, and ω=2𝜔2\omega=2italic_ω = 2, respectively. In (d), dashed lines represent the Liouvillian gap when γ0=0subscript𝛾00\gamma_{0}=0italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0, solid lines represent the maximum Liouvillian gap for γ00subscript𝛾00\gamma_{0}\neq 0italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ≠ 0. For all the plots, we set γ1=1subscript𝛾11\gamma_{1}=1italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = 1.

We first introduce an additional decay term given by the jump operator

Ll,p=2=γ2|g,le,l+1|,subscript𝐿𝑙𝑝2subscript𝛾2ket𝑔𝑙bra𝑒𝑙1L_{l,p=2}=\sqrt{\gamma_{2}}|g,l\rangle\langle e,l+1|,italic_L start_POSTSUBSCRIPT italic_l , italic_p = 2 end_POSTSUBSCRIPT = square-root start_ARG italic_γ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_ARG | italic_g , italic_l ⟩ ⟨ italic_e , italic_l + 1 | , (15)

which enhances the dissipation in the direction of the steady state. While such a term does not change the discussion on the Liouvillian superoperator under OBC, it contributes to an increased decay rate within the subsystem, effectively transforming γ1subscript𝛾1\gamma_{1}italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT to γ1+γ2subscript𝛾1subscript𝛾2\gamma_{1}+\gamma_{2}italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_γ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT in Eq. (11).Consequently, the maximum Liouvillian gap becomes γ1/4+γ2/4subscript𝛾14subscript𝛾24\gamma_{1}/4+\gamma_{2}/4italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / 4 + italic_γ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT / 4, and the pumping efficiency is enhanced. Likewise, we can introduce longer-distance decay terms to similar effects.

Alternatively, we consider the decay term Ll,0subscript𝐿𝑙0L_{l,0}italic_L start_POSTSUBSCRIPT italic_l , 0 end_POSTSUBSCRIPT, leading to transitions within each subsystem. As illustrated in Fig. 1(b), the direction of the dissipation is opposite to that of the directional flow toward the steady state. From numerical calculations, we identify two distinct behaviors of the Liouvillian gap when varying γ0subscript𝛾0\gamma_{0}italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT. First, the Liouvillian gap monotonically decreases to 0 with increasing γ0subscript𝛾0\gamma_{0}italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, shown as dashed lines in Fig. 6(a)(b). Second, the Liouvillian gap increases to a maximum value before decreasing to 0, shown as solid lines in Fig. 6(a)(b). Here we use dΔ/dγ0|γ0=0evaluated-atdΔdsubscript𝛾0subscript𝛾00\left.\mathrm{d}\Delta/\mathrm{d}\gamma_{0}\right|_{\gamma_{0}=0}roman_d roman_Δ / roman_d italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT | start_POSTSUBSCRIPT italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0 end_POSTSUBSCRIPT to differentiate the parameter regimes for these different behaviors, as shown in Fig. 6(c). When dΔ/dγ0|γ0=0<0evaluated-atdΔdsubscript𝛾0subscript𝛾000\left.\mathrm{d}\Delta/\mathrm{d}\gamma_{0}\right|_{\gamma_{0}=0}<0roman_d roman_Δ / roman_d italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT | start_POSTSUBSCRIPT italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0 end_POSTSUBSCRIPT < 0, the Liouvillian gap monotonically decreases with increasing γ0subscript𝛾0\gamma_{0}italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT; otherwise, the Liouvillian gap increases to a maximum value before decreasing to 0, resulting in a larger Liouvillian gap for appropriate values of γ0subscript𝛾0\gamma_{0}italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT compared to the case where γ0=0subscript𝛾00\gamma_{0}=0italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0. We then numerically calculate the maximum Liouvillian gap for different ω𝜔\omegaitalic_ω and ΩΩ\Omegaroman_Ω, as shown in Fig. 6(d)(e). In general, the maximum Liouvillian gap increases with larger ΩΩ\Omegaroman_Ω and smaller ω𝜔\omegaitalic_ω. The optimal decay rate γ0,maxsubscript𝛾0max\gamma_{0,\text{max}}italic_γ start_POSTSUBSCRIPT 0 , max end_POSTSUBSCRIPT for achieving the maximum Liouvillian gap shows intricate behavior in conjunction with other parameters. Introducing the decay term Ll,0subscript𝐿𝑙0L_{l,0}italic_L start_POSTSUBSCRIPT italic_l , 0 end_POSTSUBSCRIPT yields a potential maximum Liouvillian gap of γ1/2subscript𝛾12\gamma_{1}/2italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / 2, achievable under the parameters ΩΩ\Omega\rightarrow\inftyroman_Ω → ∞, ω=0𝜔0\omega=0italic_ω = 0, and γ1=γ0subscript𝛾1subscript𝛾0\gamma_{1}=\gamma_{0}italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT.

V Summary

To summarize, we show that typical optical pumping processes can be understood from the perspective of the Liouvillian skin effect. We confirm this understanding through the Liouvillian eigenspectrum and open-system dynamics for a concrete optical pumping setup involving coherent optical drives and directional dissipation. We further illustrate that such an understanding provides means to optimize the pumping efficiency. Our results are helpful for state preparation and cooling in quantum simulation and computation where optical pumping is inevitable.

Acknowledgements.
This work is supported by the National Natural Science Foundation of China (Grant No. 12374479), and by the Innovation Program for Quantum Science and Technology (Grant Nos. 2021ZD0301200, 2021ZD0301904).

Appendix A The block upper-triangular matrix

Here we discuss the eigen problem of a block upper-triangular matrix M𝑀Mitalic_M with

M=(A1,1A1,2A1,3A1,m0A2,2A2,3A2,m00A3,3000Am,m),𝑀matrixsubscript𝐴11subscript𝐴12subscript𝐴13subscript𝐴1𝑚0subscript𝐴22subscript𝐴23subscript𝐴2𝑚00subscript𝐴33000subscript𝐴𝑚𝑚M=\begin{pmatrix}A_{1,1}&A_{1,2}&A_{1,3}&\cdots&A_{1,m}\\ 0&A_{2,2}&A_{2,3}&\cdots&A_{2,m}\\ 0&0&A_{3,3}&\cdots&\vdots\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&0&\cdots&A_{m,m}\\ \end{pmatrix},italic_M = ( start_ARG start_ROW start_CELL italic_A start_POSTSUBSCRIPT 1 , 1 end_POSTSUBSCRIPT end_CELL start_CELL italic_A start_POSTSUBSCRIPT 1 , 2 end_POSTSUBSCRIPT end_CELL start_CELL italic_A start_POSTSUBSCRIPT 1 , 3 end_POSTSUBSCRIPT end_CELL start_CELL ⋯ end_CELL start_CELL italic_A start_POSTSUBSCRIPT 1 , italic_m end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL italic_A start_POSTSUBSCRIPT 2 , 2 end_POSTSUBSCRIPT end_CELL start_CELL italic_A start_POSTSUBSCRIPT 2 , 3 end_POSTSUBSCRIPT end_CELL start_CELL ⋯ end_CELL start_CELL italic_A start_POSTSUBSCRIPT 2 , italic_m end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL italic_A start_POSTSUBSCRIPT 3 , 3 end_POSTSUBSCRIPT end_CELL start_CELL ⋯ end_CELL start_CELL ⋮ end_CELL end_ROW start_ROW start_CELL ⋮ end_CELL start_CELL ⋮ end_CELL start_CELL ⋮ end_CELL start_CELL ⋱ end_CELL start_CELL ⋮ end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL ⋯ end_CELL start_CELL italic_A start_POSTSUBSCRIPT italic_m , italic_m end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ) , (16)

where Ai,jsubscript𝐴𝑖𝑗A_{i,j}italic_A start_POSTSUBSCRIPT italic_i , italic_j end_POSTSUBSCRIPT are matrices with dimensions ni×njsubscript𝑛𝑖subscript𝑛𝑗n_{i}\times n_{j}italic_n start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT × italic_n start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT, respectively.

We observe that the eigenvalues of the block upper-triangular matrix M𝑀Mitalic_M is the union of the spectra of the diagonal blocks Ai,isubscript𝐴𝑖𝑖A_{i,i}italic_A start_POSTSUBSCRIPT italic_i , italic_i end_POSTSUBSCRIPT. In the following, we prove it by induction.

First, the conclusion obviously holds for m=1𝑚1m=1italic_m = 1. Then, assuming the statement is valid for m=l𝑚𝑙m=litalic_m = italic_l, we will show below that it also holds for m=l+1𝑚𝑙1m=l+1italic_m = italic_l + 1. To simplify discussions, we set

Mlsubscript𝑀𝑙\displaystyle M_{l}italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT =(A1,1A1,2A1,3A1,l0A2,2A2,3A2,l00A3,3000Al,l),absentmatrixsubscript𝐴11subscript𝐴12subscript𝐴13subscript𝐴1𝑙0subscript𝐴22subscript𝐴23subscript𝐴2𝑙00subscript𝐴33000subscript𝐴𝑙𝑙\displaystyle=\begin{pmatrix}A_{1,1}&A_{1,2}&A_{1,3}&\cdots&A_{1,l}\\ 0&A_{2,2}&A_{2,3}&\cdots&A_{2,l}\\ 0&0&A_{3,3}&\cdots&\vdots\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&0&\cdots&A_{l,l}\\ \end{pmatrix},= ( start_ARG start_ROW start_CELL italic_A start_POSTSUBSCRIPT 1 , 1 end_POSTSUBSCRIPT end_CELL start_CELL italic_A start_POSTSUBSCRIPT 1 , 2 end_POSTSUBSCRIPT end_CELL start_CELL italic_A start_POSTSUBSCRIPT 1 , 3 end_POSTSUBSCRIPT end_CELL start_CELL ⋯ end_CELL start_CELL italic_A start_POSTSUBSCRIPT 1 , italic_l end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL italic_A start_POSTSUBSCRIPT 2 , 2 end_POSTSUBSCRIPT end_CELL start_CELL italic_A start_POSTSUBSCRIPT 2 , 3 end_POSTSUBSCRIPT end_CELL start_CELL ⋯ end_CELL start_CELL italic_A start_POSTSUBSCRIPT 2 , italic_l end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL italic_A start_POSTSUBSCRIPT 3 , 3 end_POSTSUBSCRIPT end_CELL start_CELL ⋯ end_CELL start_CELL ⋮ end_CELL end_ROW start_ROW start_CELL ⋮ end_CELL start_CELL ⋮ end_CELL start_CELL ⋮ end_CELL start_CELL ⋱ end_CELL start_CELL ⋮ end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL ⋯ end_CELL start_CELL italic_A start_POSTSUBSCRIPT italic_l , italic_l end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ) , (17)
Clsubscript𝐶𝑙\displaystyle C_{l}italic_C start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT =(A1,l+1A2,l+1Al,l+1).absentmatrixsubscript𝐴1𝑙1subscript𝐴2𝑙1subscript𝐴𝑙𝑙1\displaystyle=\begin{pmatrix}A_{1,l+1}\\ A_{2,l+1}\\ \vdots\\ A_{l,l+1}\\ \end{pmatrix}.= ( start_ARG start_ROW start_CELL italic_A start_POSTSUBSCRIPT 1 , italic_l + 1 end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL italic_A start_POSTSUBSCRIPT 2 , italic_l + 1 end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL ⋮ end_CELL end_ROW start_ROW start_CELL italic_A start_POSTSUBSCRIPT italic_l , italic_l + 1 end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ) . (18)

For m=l+1𝑚𝑙1m=l+1italic_m = italic_l + 1, we have

Ml+1=(MlCl0Al+1,l+1).subscript𝑀𝑙1matrixsubscript𝑀𝑙subscript𝐶𝑙0subscript𝐴𝑙1𝑙1\displaystyle M_{l+1}=\begin{pmatrix}M_{l}&C_{l}\\ 0&A_{l+1,l+1}\end{pmatrix}.italic_M start_POSTSUBSCRIPT italic_l + 1 end_POSTSUBSCRIPT = ( start_ARG start_ROW start_CELL italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_CELL start_CELL italic_C start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL italic_A start_POSTSUBSCRIPT italic_l + 1 , italic_l + 1 end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ) . (19)

Any eigenvalue of Mlsubscript𝑀𝑙M_{l}italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT is thus also an eigenvalue of Ml+1subscript𝑀𝑙1M_{l+1}italic_M start_POSTSUBSCRIPT italic_l + 1 end_POSTSUBSCRIPT. Specifically, for any given eigenvalue α𝛼\alphaitalic_α and the corresponding eigenvector |ψket𝜓\left|{\psi}\right>| italic_ψ ⟩ of Mlsubscript𝑀𝑙M_{l}italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT, we have

(MlCl0Al+1,l+1)(|ψ0)=(Ml|ψ0)=α(|ψ0).matrixsubscript𝑀𝑙subscript𝐶𝑙0subscript𝐴𝑙1𝑙1matrixket𝜓0matrixsubscript𝑀𝑙ket𝜓0𝛼matrixket𝜓0\begin{pmatrix}M_{l}&C_{l}\\ 0&A_{l+1,l+1}\end{pmatrix}\begin{pmatrix}\left|{\psi}\right>\\ 0\end{pmatrix}=\begin{pmatrix}M_{l}\left|{\psi}\right>\\ 0\end{pmatrix}=\alpha\begin{pmatrix}\left|{\psi}\right>\\ 0\end{pmatrix}.( start_ARG start_ROW start_CELL italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_CELL start_CELL italic_C start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL italic_A start_POSTSUBSCRIPT italic_l + 1 , italic_l + 1 end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ) ( start_ARG start_ROW start_CELL | italic_ψ ⟩ end_CELL end_ROW start_ROW start_CELL 0 end_CELL end_ROW end_ARG ) = ( start_ARG start_ROW start_CELL italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT | italic_ψ ⟩ end_CELL end_ROW start_ROW start_CELL 0 end_CELL end_ROW end_ARG ) = italic_α ( start_ARG start_ROW start_CELL | italic_ψ ⟩ end_CELL end_ROW start_ROW start_CELL 0 end_CELL end_ROW end_ARG ) . (20)

We then show that any eigenvalue of Al+1,l+1subscript𝐴𝑙1𝑙1A_{l+1,l+1}italic_A start_POSTSUBSCRIPT italic_l + 1 , italic_l + 1 end_POSTSUBSCRIPT is also an eigenvalue of Ml+1subscript𝑀𝑙1M_{l+1}italic_M start_POSTSUBSCRIPT italic_l + 1 end_POSTSUBSCRIPT. For that purpose, we focus on a given eigenvalue al+1superscript𝑎𝑙1a^{l+1}italic_a start_POSTSUPERSCRIPT italic_l + 1 end_POSTSUPERSCRIPT and the corresponding eigenstate |ψl+1ketsuperscript𝜓𝑙1\left|{\psi^{l+1}}\right>| italic_ψ start_POSTSUPERSCRIPT italic_l + 1 end_POSTSUPERSCRIPT ⟩ of Al+1,l+1subscript𝐴𝑙1𝑙1A_{l+1,l+1}italic_A start_POSTSUBSCRIPT italic_l + 1 , italic_l + 1 end_POSTSUBSCRIPT. If al+1superscript𝑎𝑙1a^{l+1}italic_a start_POSTSUPERSCRIPT italic_l + 1 end_POSTSUPERSCRIPT is also an eigenvalue of Mlsubscript𝑀𝑙M_{l}italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT, with the corresponding eigenstate |ϕ0ketsubscriptitalic-ϕ0\left|{\phi_{0}}\right>| italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ⟩, we have

(MlCl0Al+1,l+1)(|ϕ00)=(Ml|ϕ00)=al+1(|ϕ00).matrixsubscript𝑀𝑙subscript𝐶𝑙0subscript𝐴𝑙1𝑙1matrixketsubscriptitalic-ϕ00matrixsubscript𝑀𝑙ketsubscriptitalic-ϕ00superscript𝑎𝑙1matrixketsubscriptitalic-ϕ00\begin{split}\begin{pmatrix}M_{l}&C_{l}\\ 0&A_{l+1,l+1}\end{pmatrix}\begin{pmatrix}\left|{\phi_{0}}\right>\\ 0\end{pmatrix}&=\begin{pmatrix}M_{l}\left|{\phi_{0}}\right>\\ 0\end{pmatrix}\\ &=a^{l+1}\begin{pmatrix}\left|{\phi_{0}}\right>\\ 0\end{pmatrix}\end{split}.start_ROW start_CELL ( start_ARG start_ROW start_CELL italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_CELL start_CELL italic_C start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL italic_A start_POSTSUBSCRIPT italic_l + 1 , italic_l + 1 end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ) ( start_ARG start_ROW start_CELL | italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ⟩ end_CELL end_ROW start_ROW start_CELL 0 end_CELL end_ROW end_ARG ) end_CELL start_CELL = ( start_ARG start_ROW start_CELL italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT | italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ⟩ end_CELL end_ROW start_ROW start_CELL 0 end_CELL end_ROW end_ARG ) end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL = italic_a start_POSTSUPERSCRIPT italic_l + 1 end_POSTSUPERSCRIPT ( start_ARG start_ROW start_CELL | italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ⟩ end_CELL end_ROW start_ROW start_CELL 0 end_CELL end_ROW end_ARG ) end_CELL end_ROW . (21)

Otherwise, we have

(MlCl0Al+1,l+1)(|ϕ|ψl+1)=(Ml|ϕ+Cl|ψl+1Al+1,l+1|ψl+1)=(Ml|ϕ+Cl|ψl+1al+1|ψl+1),matrixsubscript𝑀𝑙subscript𝐶𝑙0subscript𝐴𝑙1𝑙1matrixketitalic-ϕketsuperscript𝜓𝑙1matrixsubscript𝑀𝑙ketitalic-ϕsubscript𝐶𝑙ketsuperscript𝜓𝑙1subscript𝐴𝑙1𝑙1ketsuperscript𝜓𝑙1matrixsubscript𝑀𝑙ketitalic-ϕsubscript𝐶𝑙ketsuperscript𝜓𝑙1superscript𝑎𝑙1ketsuperscript𝜓𝑙1\begin{split}\begin{pmatrix}M_{l}&C_{l}\\ 0&A_{l+1,l+1}\end{pmatrix}\begin{pmatrix}\left|{\phi}\right>\\ \left|{\psi^{l+1}}\right>\end{pmatrix}&=\begin{pmatrix}M_{l}\left|{\phi}\right% >+C_{l}\left|{\psi^{l+1}}\right>\\ A_{l+1,l+1}\left|{\psi^{l+1}}\right>\end{pmatrix}\\ &=\begin{pmatrix}M_{l}\left|{\phi}\right>+C_{l}\left|{\psi^{l+1}}\right>\\ a^{l+1}\left|{\psi^{l+1}}\right>\end{pmatrix}\end{split},start_ROW start_CELL ( start_ARG start_ROW start_CELL italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_CELL start_CELL italic_C start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL italic_A start_POSTSUBSCRIPT italic_l + 1 , italic_l + 1 end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ) ( start_ARG start_ROW start_CELL | italic_ϕ ⟩ end_CELL end_ROW start_ROW start_CELL | italic_ψ start_POSTSUPERSCRIPT italic_l + 1 end_POSTSUPERSCRIPT ⟩ end_CELL end_ROW end_ARG ) end_CELL start_CELL = ( start_ARG start_ROW start_CELL italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT | italic_ϕ ⟩ + italic_C start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT | italic_ψ start_POSTSUPERSCRIPT italic_l + 1 end_POSTSUPERSCRIPT ⟩ end_CELL end_ROW start_ROW start_CELL italic_A start_POSTSUBSCRIPT italic_l + 1 , italic_l + 1 end_POSTSUBSCRIPT | italic_ψ start_POSTSUPERSCRIPT italic_l + 1 end_POSTSUPERSCRIPT ⟩ end_CELL end_ROW end_ARG ) end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL = ( start_ARG start_ROW start_CELL italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT | italic_ϕ ⟩ + italic_C start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT | italic_ψ start_POSTSUPERSCRIPT italic_l + 1 end_POSTSUPERSCRIPT ⟩ end_CELL end_ROW start_ROW start_CELL italic_a start_POSTSUPERSCRIPT italic_l + 1 end_POSTSUPERSCRIPT | italic_ψ start_POSTSUPERSCRIPT italic_l + 1 end_POSTSUPERSCRIPT ⟩ end_CELL end_ROW end_ARG ) end_CELL end_ROW , (22)

where |ϕketitalic-ϕ\left|{\phi}\right>| italic_ϕ ⟩ is an unknown state. We set

Ml|ϕal+1|ϕ=(Mlal+1)|ϕ=Cl|ψl+1,subscript𝑀𝑙ketitalic-ϕsuperscript𝑎𝑙1ketitalic-ϕsubscript𝑀𝑙superscript𝑎𝑙1ketitalic-ϕsubscript𝐶𝑙ketsuperscript𝜓𝑙1M_{l}\left|{\phi}\right>-a^{l+1}\left|{\phi}\right>=(M_{l}-a^{l+1}\mathcal{I})% \left|{\phi}\right>=C_{l}\left|{\psi^{l+1}}\right>,italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT | italic_ϕ ⟩ - italic_a start_POSTSUPERSCRIPT italic_l + 1 end_POSTSUPERSCRIPT | italic_ϕ ⟩ = ( italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT - italic_a start_POSTSUPERSCRIPT italic_l + 1 end_POSTSUPERSCRIPT caligraphic_I ) | italic_ϕ ⟩ = italic_C start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT | italic_ψ start_POSTSUPERSCRIPT italic_l + 1 end_POSTSUPERSCRIPT ⟩ , (23)

where \mathcal{I}caligraphic_I is an identity matrix. Since al+1superscript𝑎𝑙1a^{l+1}italic_a start_POSTSUPERSCRIPT italic_l + 1 end_POSTSUPERSCRIPT is not an eigenvalue of Mlsubscript𝑀𝑙M_{l}italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT, (Mlal+1)subscript𝑀𝑙superscript𝑎𝑙1(M_{l}-a^{l+1}\mathcal{I})( italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT - italic_a start_POSTSUPERSCRIPT italic_l + 1 end_POSTSUPERSCRIPT caligraphic_I ) is reversible, and Eq. 23 must have nontrivial solutions. In other words, we can always find |ϕketitalic-ϕ\left|{\phi}\right>| italic_ϕ ⟩ such that Eq. 22 is satisfied, yielding the right eigenstate of Ml+1subscript𝑀𝑙1M_{l+1}italic_M start_POSTSUBSCRIPT italic_l + 1 end_POSTSUBSCRIPT.

In summary, the eigenvalue of Mlsubscript𝑀𝑙M_{l}italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT and Al+1,l+1subscript𝐴𝑙1𝑙1A_{l+1,l+1}italic_A start_POSTSUBSCRIPT italic_l + 1 , italic_l + 1 end_POSTSUBSCRIPT are the eigenvalues of Ml+1subscript𝑀𝑙1M_{l+1}italic_M start_POSTSUBSCRIPT italic_l + 1 end_POSTSUBSCRIPT. In other words, our statement is also valid for m=l+1𝑚𝑙1m=l+1italic_m = italic_l + 1. We have therefore proved our statement by induction, that the eigenvalues of the block upper-triangular matrix M𝑀Mitalic_M is the union of the spectra of the diagonal blocks Ai,isubscript𝐴𝑖𝑖A_{i,i}italic_A start_POSTSUBSCRIPT italic_i , italic_i end_POSTSUBSCRIPT.

Furthermore, we notice that the right eigenstates of the block upper-triangular matrix are usually localized in the Hilbert space. Here we provide a simple explanation. We set

Mk=(A1,1A1,2A1,3A1,k0A2,2A2,3A2,k00A3,3000Ak,k),Dk=(A1,l+1A1,l+1A2,l+1A1,l+1Al,l+1A1,l+1),Mk=(Ak+1,k+1Ak+1,k+2Ak+1,k+3Ak+1,m0Ak+2,k+2Ak+2,k+3Ak+2,m00Ak+3,k+3000Am,m).formulae-sequencesubscript𝑀𝑘matrixsubscript𝐴11subscript𝐴12subscript𝐴13subscript𝐴1𝑘0subscript𝐴22subscript𝐴23subscript𝐴2𝑘00subscript𝐴33000subscript𝐴𝑘𝑘formulae-sequencesubscript𝐷𝑘matrixsubscript𝐴1𝑙1subscript𝐴1𝑙1subscript𝐴2𝑙1subscript𝐴1𝑙1subscript𝐴𝑙𝑙1subscript𝐴1𝑙1subscriptsuperscript𝑀𝑘matrixsubscript𝐴𝑘1𝑘1subscript𝐴𝑘1𝑘2subscript𝐴𝑘1𝑘3subscript𝐴𝑘1𝑚0subscript𝐴𝑘2𝑘2subscript𝐴𝑘2𝑘3subscript𝐴𝑘2𝑚00subscript𝐴𝑘3𝑘3000subscript𝐴𝑚𝑚\begin{split}&M_{k}=\begin{pmatrix}A_{1,1}&A_{1,2}&A_{1,3}&\cdots&A_{1,k}\\ 0&A_{2,2}&A_{2,3}&\cdots&A_{2,k}\\ 0&0&A_{3,3}&\cdots&\vdots\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&0&\cdots&A_{k,k}\\ \end{pmatrix},\\ &D_{k}=\begin{pmatrix}A_{1,l+1}&\cdots&A_{1,l+1}\\ A_{2,l+1}&\cdots&A_{1,l+1}\\ \vdots\\ A_{l,l+1}&\cdots&A_{1,l+1}\\ \end{pmatrix},\\ &M^{\prime}_{k}=\begin{pmatrix}A_{k+1,k+1}&A_{k+1,k+2}&A_{k+1,k+3}&\cdots&A_{k% +1,m}\\ 0&A_{k+2,k+2}&A_{k+2,k+3}&\cdots&A_{k+2,m}\\ 0&0&A_{k+3,k+3}&\cdots&\vdots\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&0&\cdots&A_{m,m}\\ \end{pmatrix}\end{split}.start_ROW start_CELL end_CELL start_CELL italic_M start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT = ( start_ARG start_ROW start_CELL italic_A start_POSTSUBSCRIPT 1 , 1 end_POSTSUBSCRIPT end_CELL start_CELL italic_A start_POSTSUBSCRIPT 1 , 2 end_POSTSUBSCRIPT end_CELL start_CELL italic_A start_POSTSUBSCRIPT 1 , 3 end_POSTSUBSCRIPT end_CELL start_CELL ⋯ end_CELL start_CELL italic_A start_POSTSUBSCRIPT 1 , italic_k end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL italic_A start_POSTSUBSCRIPT 2 , 2 end_POSTSUBSCRIPT end_CELL start_CELL italic_A start_POSTSUBSCRIPT 2 , 3 end_POSTSUBSCRIPT end_CELL start_CELL ⋯ end_CELL start_CELL italic_A start_POSTSUBSCRIPT 2 , italic_k end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL italic_A start_POSTSUBSCRIPT 3 , 3 end_POSTSUBSCRIPT end_CELL start_CELL ⋯ end_CELL start_CELL ⋮ end_CELL end_ROW start_ROW start_CELL ⋮ end_CELL start_CELL ⋮ end_CELL start_CELL ⋮ end_CELL start_CELL ⋱ end_CELL start_CELL ⋮ end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL ⋯ end_CELL start_CELL italic_A start_POSTSUBSCRIPT italic_k , italic_k end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ) , end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL italic_D start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT = ( start_ARG start_ROW start_CELL italic_A start_POSTSUBSCRIPT 1 , italic_l + 1 end_POSTSUBSCRIPT end_CELL start_CELL ⋯ end_CELL start_CELL italic_A start_POSTSUBSCRIPT 1 , italic_l + 1 end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL italic_A start_POSTSUBSCRIPT 2 , italic_l + 1 end_POSTSUBSCRIPT end_CELL start_CELL ⋯ end_CELL start_CELL italic_A start_POSTSUBSCRIPT 1 , italic_l + 1 end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL ⋮ end_CELL end_ROW start_ROW start_CELL italic_A start_POSTSUBSCRIPT italic_l , italic_l + 1 end_POSTSUBSCRIPT end_CELL start_CELL ⋯ end_CELL start_CELL italic_A start_POSTSUBSCRIPT 1 , italic_l + 1 end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ) , end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL italic_M start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT = ( start_ARG start_ROW start_CELL italic_A start_POSTSUBSCRIPT italic_k + 1 , italic_k + 1 end_POSTSUBSCRIPT end_CELL start_CELL italic_A start_POSTSUBSCRIPT italic_k + 1 , italic_k + 2 end_POSTSUBSCRIPT end_CELL start_CELL italic_A start_POSTSUBSCRIPT italic_k + 1 , italic_k + 3 end_POSTSUBSCRIPT end_CELL start_CELL ⋯ end_CELL start_CELL italic_A start_POSTSUBSCRIPT italic_k + 1 , italic_m end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL italic_A start_POSTSUBSCRIPT italic_k + 2 , italic_k + 2 end_POSTSUBSCRIPT end_CELL start_CELL italic_A start_POSTSUBSCRIPT italic_k + 2 , italic_k + 3 end_POSTSUBSCRIPT end_CELL start_CELL ⋯ end_CELL start_CELL italic_A start_POSTSUBSCRIPT italic_k + 2 , italic_m end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL italic_A start_POSTSUBSCRIPT italic_k + 3 , italic_k + 3 end_POSTSUBSCRIPT end_CELL start_CELL ⋯ end_CELL start_CELL ⋮ end_CELL end_ROW start_ROW start_CELL ⋮ end_CELL start_CELL ⋮ end_CELL start_CELL ⋮ end_CELL start_CELL ⋱ end_CELL start_CELL ⋮ end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL ⋯ end_CELL start_CELL italic_A start_POSTSUBSCRIPT italic_m , italic_m end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ) end_CELL end_ROW . (24)

the entire matrix follows

Ml=(MkDk0Mk).subscript𝑀𝑙matrixsubscript𝑀𝑘subscript𝐷𝑘0subscriptsuperscript𝑀𝑘M_{l}=\begin{pmatrix}M_{k}&D_{k}\\ 0&M^{\prime}_{k}\end{pmatrix}.italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = ( start_ARG start_ROW start_CELL italic_M start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT end_CELL start_CELL italic_D start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL italic_M start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ) . (25)

Following the previous discussion, for any k𝑘kitalic_k, the eigenvalues of Mksubscript𝑀𝑘M_{k}italic_M start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT are also the eigenvalues of Mlsubscript𝑀𝑙M_{l}italic_M start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT and the corresponding right eigenstate is equal to 0 in the space of Mksubscriptsuperscript𝑀𝑘M^{\prime}_{k}italic_M start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT. Therefore, these right eigenstates are localized within the space of Mksubscript𝑀𝑘M_{k}italic_M start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT.

Appendix B Liouvillian gap

Here we provide analytic expressions for the Liouvillian gap in the main text.

Under the OBC, when γ0=0subscript𝛾00\gamma_{0}=0italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0, ωl=0subscript𝜔𝑙0\omega_{l}=0italic_ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = 0 and Ωl=ΩsubscriptΩ𝑙Ω\Omega_{l}=\Omegaroman_Ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = roman_Ω, the Liouvillian gap follows

ΔOBC={14(γ1γ1216Ω2),forΩγ1<14,14γ1,forΩγ114.subscriptΔ𝑂𝐵𝐶cases14subscript𝛾1superscriptsubscript𝛾1216superscriptΩ2forΩsubscript𝛾11414subscript𝛾1forΩsubscript𝛾114\Delta_{OBC}=\begin{cases}\frac{1}{4}(\gamma_{1}-\sqrt{\gamma_{1}^{2}-16\Omega% ^{2}}),&\text{for}\,\,\frac{\Omega}{\gamma_{1}}<\frac{1}{4},\\ \frac{1}{4}\gamma_{1},&\text{for}\,\,\frac{\Omega}{\gamma_{1}}\geq\frac{1}{4}.% \end{cases}roman_Δ start_POSTSUBSCRIPT italic_O italic_B italic_C end_POSTSUBSCRIPT = { start_ROW start_CELL divide start_ARG 1 end_ARG start_ARG 4 end_ARG ( italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - square-root start_ARG italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 16 roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ) , end_CELL start_CELL for divide start_ARG roman_Ω end_ARG start_ARG italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG < divide start_ARG 1 end_ARG start_ARG 4 end_ARG , end_CELL end_ROW start_ROW start_CELL divide start_ARG 1 end_ARG start_ARG 4 end_ARG italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , end_CELL start_CELL for divide start_ARG roman_Ω end_ARG start_ARG italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG ≥ divide start_ARG 1 end_ARG start_ARG 4 end_ARG . end_CELL end_ROW (26)

If we consider ωl=ωsubscript𝜔𝑙𝜔\omega_{l}=\omegaitalic_ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = italic_ω, the Liouvillian gap becomes

Δ=14(γ1Im[γ124ω24iωγ116Ω2]).Δ14subscript𝛾1Imdelimited-[]superscriptsubscript𝛾124superscript𝜔24𝑖𝜔subscript𝛾116superscriptΩ2\Delta=\frac{1}{4}(\gamma_{1}-\text{Im}\left[\sqrt{\gamma_{1}^{2}-4\omega^{2}-% 4i\omega\gamma_{1}-16\Omega^{2}}\right]).roman_Δ = divide start_ARG 1 end_ARG start_ARG 4 end_ARG ( italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - Im [ square-root start_ARG italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 4 italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 4 italic_i italic_ω italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - 16 roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ] ) . (27)

Under the PBC, when ωl=0subscript𝜔𝑙0\omega_{l}=0italic_ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = 0, Ωl=ΩsubscriptΩ𝑙Ω\Omega_{l}=\Omegaroman_Ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = roman_Ω and γ0=0subscript𝛾00\gamma_{0}=0italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0, we regard 0subscript0\mathcal{L}_{0}caligraphic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT as a four-band, one-dimensional lattice. Due to the lattice translational symmetry, its Hamilton can be written in the k𝑘kitalic_k space as

H4=(γ1iΩiΩγ1eikiΩ12γ10iΩiΩ012γ1iΩ0iΩiΩ0).subscript𝐻4matrixsubscript𝛾1𝑖Ω𝑖Ωsubscript𝛾1superscript𝑒𝑖𝑘𝑖Ω12subscript𝛾10𝑖Ω𝑖Ω012subscript𝛾1𝑖Ω0𝑖Ω𝑖Ω0H_{4}=\begin{pmatrix}-\gamma_{1}&i\Omega&-i\Omega&\gamma_{1}e^{ik}\\ i\Omega&-\frac{1}{2}\gamma_{1}&0&-i\Omega\\ -i\Omega&0&-\frac{1}{2}\gamma_{1}&i\Omega\\ 0&-i\Omega&i\Omega&0\end{pmatrix}.italic_H start_POSTSUBSCRIPT 4 end_POSTSUBSCRIPT = ( start_ARG start_ROW start_CELL - italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_CELL start_CELL italic_i roman_Ω end_CELL start_CELL - italic_i roman_Ω end_CELL start_CELL italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_e start_POSTSUPERSCRIPT italic_i italic_k end_POSTSUPERSCRIPT end_CELL end_ROW start_ROW start_CELL italic_i roman_Ω end_CELL start_CELL - divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_CELL start_CELL 0 end_CELL start_CELL - italic_i roman_Ω end_CELL end_ROW start_ROW start_CELL - italic_i roman_Ω end_CELL start_CELL 0 end_CELL start_CELL - divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_CELL start_CELL italic_i roman_Ω end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL - italic_i roman_Ω end_CELL start_CELL italic_i roman_Ω end_CELL start_CELL 0 end_CELL end_ROW end_ARG ) . (28)

Here k=mπ/N(m=1,2,,N/2)𝑘𝑚𝜋𝑁𝑚12𝑁2k=m\pi/N(m=1,2,\cdots,N/2)italic_k = italic_m italic_π / italic_N ( italic_m = 1 , 2 , ⋯ , italic_N / 2 ) is the lattice momentum. The Liouvillian gap can be calculated from the spectrum of Eq. 28.

We then calculate the Liouvillian gap after introducing the decay term Ll,1subscript𝐿𝑙1L_{l,1}italic_L start_POSTSUBSCRIPT italic_l , 1 end_POSTSUBSCRIPT, under the OBC and with ωl=ωsubscript𝜔𝑙𝜔\omega_{l}=\omegaitalic_ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = italic_ω and Ωl=ΩsubscriptΩ𝑙Ω\Omega_{l}=\Omegaroman_Ω start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = roman_Ω.

In order to derive the Liouvillian gap, we need to calculate the spectrum of each diagonal block in 0subscript0\mathcal{L}_{0}caligraphic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, as well as the blocks ljsubscript𝑙𝑗\mathcal{L}_{lj}caligraphic_L start_POSTSUBSCRIPT italic_l italic_j end_POSTSUBSCRIPT. When ω0𝜔0\omega\neq 0italic_ω ≠ 0, the expression of the Liouvillian gap is extremely complicated. However, we notice that the Liouvillian gap consistently decreases with increasing energy interval ω𝜔\omegaitalic_ω. Thus, we calculate the Liouvillian gap for ω=0𝜔0\omega=0italic_ω = 0 for an upper bound, which is given by

Δ={14(γ0+γ1(γ0+γ1)216Ω2),for  4Ω<γ0+γ1,527γ02+575γ12+1166γ0γ1>9216Ω2,14(γ0+γ1),for  4Ωγ0+γ1,64Ω2(γ1γ0)>3(γ1+γ0)3,γ0+γ1272γ0Ω2+1346656γ02Ω4+(3γ023γ126γ0γ2+48Ω2)332 32/3γ02+γ12+2γ0γ116Ω223372γ0Ω2+1346656γ02Ω4+(3γ023γ126γ0γ1+48Ω2)33,otherwise.\Delta=\begin{cases}\frac{1}{4}(\gamma_{0}+\gamma_{1}-\sqrt{(\gamma_{0}+\gamma% _{1})^{2}-16\Omega^{2}}),&\text{for}\,\,4\Omega<\gamma_{0}+\gamma_{1},527% \gamma_{0}^{2}+575\gamma_{1}^{2}+1166\gamma_{0}\gamma_{1}>9216\Omega^{2},\\ \frac{1}{4}(\gamma_{0}+\gamma_{1}),&\text{for}\,\,4\Omega\geq\gamma_{0}+\gamma% _{1},64\Omega^{2}(\gamma_{1}-\gamma_{0})>3(\gamma_{1}+\gamma_{0})^{3},\\ \frac{\gamma_{0}+\gamma_{1}}{2}-\frac{\sqrt[3]{72\gamma_{0}\Omega^{2}+\frac{1}% {3}\sqrt{46656\gamma_{0}^{2}\Omega^{4}+\left(-3\gamma_{0}^{2}-3\gamma_{1}^{2}-% 6\gamma_{0}\gamma_{2}+48\Omega^{2}\right){}^{3}}}}{2\ 3^{2/3}}\\ -\frac{\gamma_{0}^{2}+\gamma_{1}^{2}+2\gamma_{0}\gamma_{1}-16\Omega^{2}}{2% \sqrt[3]{3}\sqrt[3]{72\gamma_{0}\Omega^{2}+\frac{1}{3}\sqrt{46656\gamma_{0}^{2% }\Omega^{4}+\left(-3\gamma_{0}^{2}-3\gamma_{1}^{2}-6\gamma_{0}\gamma_{1}+48% \Omega^{2}\right){}^{3}}}},&\text{otherwise}.\end{cases}roman_Δ = { start_ROW start_CELL divide start_ARG 1 end_ARG start_ARG 4 end_ARG ( italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - square-root start_ARG ( italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 16 roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ) , end_CELL start_CELL for 4 roman_Ω < italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , 527 italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + 575 italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + 1166 italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT > 9216 roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT , end_CELL end_ROW start_ROW start_CELL divide start_ARG 1 end_ARG start_ARG 4 end_ARG ( italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ) , end_CELL start_CELL for 4 roman_Ω ≥ italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , 64 roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) > 3 ( italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT , end_CELL end_ROW start_ROW start_CELL divide start_ARG italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG - divide start_ARG nth-root start_ARG 3 end_ARG start_ARG 72 italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + divide start_ARG 1 end_ARG start_ARG 3 end_ARG square-root start_ARG 46656 italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT + ( - 3 italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 3 italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 6 italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_γ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT + 48 roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_FLOATSUPERSCRIPT 3 end_FLOATSUPERSCRIPT end_ARG end_ARG end_ARG start_ARG 2 3 start_POSTSUPERSCRIPT 2 / 3 end_POSTSUPERSCRIPT end_ARG end_CELL start_CELL end_CELL end_ROW start_ROW start_CELL - divide start_ARG italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + 2 italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - 16 roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 nth-root start_ARG 3 end_ARG start_ARG 3 end_ARG nth-root start_ARG 3 end_ARG start_ARG 72 italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + divide start_ARG 1 end_ARG start_ARG 3 end_ARG square-root start_ARG 46656 italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT + ( - 3 italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 3 italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 6 italic_γ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + 48 roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_FLOATSUPERSCRIPT 3 end_FLOATSUPERSCRIPT end_ARG end_ARG end_ARG , end_CELL start_CELL otherwise . end_CELL end_ROW (29)

According to Eq. 29, for certain ΩΩ\Omegaroman_Ω and γ1subscript𝛾1\gamma_{1}italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , we have

γ1,max={0,for  4Ω<γ1,4Ωγ1,for  4Ωγ172Ω,4381γ12Ω4+64Ω69γ1Ω23+16Ω2381γ12Ω4+64Ω69γ1Ω23γ1,forγ1<72Ω.subscript𝛾1𝑚𝑎𝑥cases0for4Ωsubscript𝛾14Ωsubscript𝛾1for4Ωsubscript𝛾172Ω43381superscriptsubscript𝛾12superscriptΩ464superscriptΩ69subscript𝛾1superscriptΩ216superscriptΩ23381superscriptsubscript𝛾12superscriptΩ464superscriptΩ69subscript𝛾1superscriptΩ2subscript𝛾1forsubscript𝛾172Ω\gamma_{1,max}=\begin{cases}0,&\text{for}\,\,4\Omega<\gamma_{1},\\ 4\Omega-\gamma_{1},&\text{for}\,\,4\Omega\geq\gamma_{1}\geq\frac{7}{2}\Omega,% \\ -\frac{4}{3}\sqrt[3]{\sqrt{81\gamma_{1}^{2}\Omega^{4}+64\Omega^{6}}-9\gamma_{1% }\Omega^{2}}+\frac{16\Omega^{2}}{3\sqrt[3]{\sqrt{81\gamma_{1}^{2}\Omega^{4}+64% \Omega^{6}}-9\gamma_{1}\Omega^{2}}}-\gamma_{1},&\text{for}\,\,\gamma_{1}<\frac% {7}{2}\Omega.\end{cases}italic_γ start_POSTSUBSCRIPT 1 , italic_m italic_a italic_x end_POSTSUBSCRIPT = { start_ROW start_CELL 0 , end_CELL start_CELL for 4 roman_Ω < italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , end_CELL end_ROW start_ROW start_CELL 4 roman_Ω - italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , end_CELL start_CELL for 4 roman_Ω ≥ italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ≥ divide start_ARG 7 end_ARG start_ARG 2 end_ARG roman_Ω , end_CELL end_ROW start_ROW start_CELL - divide start_ARG 4 end_ARG start_ARG 3 end_ARG nth-root start_ARG 3 end_ARG start_ARG square-root start_ARG 81 italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT + 64 roman_Ω start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT end_ARG - 9 italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + divide start_ARG 16 roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 3 nth-root start_ARG 3 end_ARG start_ARG square-root start_ARG 81 italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT + 64 roman_Ω start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT end_ARG - 9 italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG end_ARG - italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , end_CELL start_CELL for italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT < divide start_ARG 7 end_ARG start_ARG 2 end_ARG roman_Ω . end_CELL end_ROW (30)

The maximum Liouvillian gap is therefore

Δmax(γ1,Ω)={14(γ1γ1216Ω2),for  4Ω<γ1,Ω,for  4Ωγ172Ω,1381γ12Ω4+64Ω69γ1Ω23+4Ω2381γ12Ω4+64Ω69γ1Ω23,forγ1<72Ω.subscriptΔ𝑚𝑎𝑥subscript𝛾1Ωcases14subscript𝛾1superscriptsubscript𝛾1216superscriptΩ2for4Ωsubscript𝛾1Ωfor4Ωsubscript𝛾172Ω13381superscriptsubscript𝛾12superscriptΩ464superscriptΩ69subscript𝛾1superscriptΩ24superscriptΩ23381superscriptsubscript𝛾12superscriptΩ464superscriptΩ69subscript𝛾1superscriptΩ2forsubscript𝛾172Ω\Delta_{max}(\gamma_{1},\Omega)=\begin{cases}\frac{1}{4}(\gamma_{1}-\sqrt{% \gamma_{1}^{2}-16\Omega^{2}}),&\text{for}\,\,4\Omega<\gamma_{1},\\ \Omega,&\text{for}\,\,4\Omega\geq\gamma_{1}\geq\frac{7}{2}\Omega,\\ -\frac{1}{3}\sqrt[3]{\sqrt{81\gamma_{1}^{2}\Omega^{4}+64\Omega^{6}}-9\gamma_{1% }\Omega^{2}}+\frac{4\Omega^{2}}{3\sqrt[3]{\sqrt{81\gamma_{1}^{2}\Omega^{4}+64% \Omega^{6}}-9\gamma_{1}\Omega^{2}}},&\text{for}\,\,\gamma_{1}<\frac{7}{2}% \Omega.\end{cases}roman_Δ start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT ( italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , roman_Ω ) = { start_ROW start_CELL divide start_ARG 1 end_ARG start_ARG 4 end_ARG ( italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - square-root start_ARG italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 16 roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ) , end_CELL start_CELL for 4 roman_Ω < italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , end_CELL end_ROW start_ROW start_CELL roman_Ω , end_CELL start_CELL for 4 roman_Ω ≥ italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ≥ divide start_ARG 7 end_ARG start_ARG 2 end_ARG roman_Ω , end_CELL end_ROW start_ROW start_CELL - divide start_ARG 1 end_ARG start_ARG 3 end_ARG nth-root start_ARG 3 end_ARG start_ARG square-root start_ARG 81 italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT + 64 roman_Ω start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT end_ARG - 9 italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + divide start_ARG 4 roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 3 nth-root start_ARG 3 end_ARG start_ARG square-root start_ARG 81 italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT + 64 roman_Ω start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT end_ARG - 9 italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG end_ARG , end_CELL start_CELL for italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT < divide start_ARG 7 end_ARG start_ARG 2 end_ARG roman_Ω . end_CELL end_ROW (31)

When ΩΩ\Omega\to\inftyroman_Ω → ∞, the maximum possible Liouvillian gap is γ1/2subscript𝛾12\gamma_{1}/2italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / 2.

References

  • [1] W. Franzen and A. G. Emslie, Phys. Rev. 108, 1453 (1957).
  • [2] C. Cohen-Tannoudji, and A. Kastler, Pro. Opt. 5, 1 (1966).
  • [3] W. Happer, Rev. Mod. Phys. 44, 169 (1972).
  • [4] W. Happer and B. S. Mathur, Phys. Rev. 163, 12 (1967).
  • [5] W. Happer, E. A. Miron, S. Schaefer, D. Schreiber, W. A. van Wijngaarden, X. Zeng, Phys. Rev. A 29, 3092 (1984).
  • [6] T. G. Walker, and W. Happer, Rev. Mod. Phys. 69, 629 (1997).
  • [7] S. Appelt, A. Ben-Amar Baranga, C. J. Erickson, M. V. Romalis, A. R. Young, and W. Happer, Phys. Rev. A 58, 1412 (1998).
  • [8] J. Han, M.C. Heaven, Opt. Lett. 37, 2157 (2012).
  • [9] R. E. Drullinger, R. N. Zare, J. Chem. Phys. 51, 5532 (1969).
  • [10] M. Broyer, G. Gouedard, J.C. Lehmann, J. Vigue, Adv. Atom. Mole. Phys. 12, 165 (1976).
  • [11] M. Viteau, A. Chotia, M. Allegrini, N. Bouloufa, O. Dulieu, D. Comparat, and P. Pillet, Science 321, 232 (2008).
  • [12] L. C. Balling, R. J. Hanson, and F. M. Pipkin, Phys. Rev. 133, A607 (1964).
  • [13] B. A. Olsen, B. Patton, Y.-Y. Jau, and W. Happer, Phys. Rev. A 84, 063410 (2011).
  • [14] G. A. Pitz, and M. D. Anderson, Appl. Phys. Rev. 4, 041101 (2017).
  • [15] Y.-Y. Jau, E. Miron, A. B. Post, N. N. Kuzma, and W. Happer, Phys. Rev. Lett. 93, 160802 (2004).
  • [16] E. W. Weber, Phys. Rep. 32, 123 (1977).
  • [17] S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. Büchler,and P. Zoller, Nat. Phys. 4, 878 (2008).
  • [18] S. Diehl, E. Rico, M. A. Baranov, and P. Zoller, Nat. Phys. 7, 971 (2011).
  • [19] K. Stannige, P. Rabl, and P. Zoller, New J. Phys. 14, 063014 (2012).
  • [20] A. Tomadin, S. Diehl, and P. Zoller, Phys. Rev. A 83, 013611 (2011).
  • [21] S. Zhang, J.-Q. Zhang, W. Wu, W.-S. Bao and C. Guo, New J. Phys. 23, 023018 (2021).
  • [22] Z. Lin, Y. Lin, and W. Yi, Phys. Rev. A 106, 063112 (2022).
  • [23] Z. Wang, Y. Lu, Y. Peng, R. Qi, Y. Wang, and J. Jie, Phys. Rev. B 108, 054313 (2023).
  • [24] D. J. Wineland, C. Monroe, W. M. Itano, D. Leibfried, B. E. King and D. M. Meekhof, J. Res. Natl Inst. Stand. Technol. 103, 259 (1998).
  • [25] F. Diedrich, J. Bergquist, W. M. Itano, and D. Wineland, Phys. Rev. Lett. 62, 403 (1989).
  • [26] C. Monroe, D. Meekhof, B. King, S. R. Jefferts, W. M. Itano, D. J. Wineland, and P. Gould, Phys. Rev. Lett. 75, 4011 (1995).
  • [27] Ch. Roos, Th. Zeiger, H. Rohde, H. C. Nägerl, J. Eschner, D. Leibfried, F. Schmidt-Kaler, and R. Blatt, Phys. Rev. Lett. 83, 4713 (1999).
  • [28] D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, Rev. Mod. Phys. 75, 281 (2003).
  • [29] S. Yao and Z. Wang, Phys. Rev. Lett. 121, 086803 (2018).
  • [30] S. Yao, F. Song, and Z. Wang, Phys. Rev. Lett. 121, 136802 (2018).
  • [31] F. Song, S. Yao, and Z. Wang, Phys. Rev. Lett. 123, 246801 (2019).
  • [32] F. K. Kunst, E. Edvardsson, J. C. Budich, and E. J. Bergholtz, Phys. Rev. Lett. 121,026808 (2018).
  • [33] K. Yokomizo and S. Murakami, Phys. Rev. Lett. 123, 066404 (2019).
  • [34] C. H. Lee and R. Thomale, Phys. Rev. B 99, 201103(R)(2019).
  • [35] T.-S. Deng and W. Yi, Phys. Rev. B 100, 035102 (2019).
  • [36] S. Longhi, Phys. Rev. Research 1, 023013 (2019).
  • [37] T. Li, J.-Z. Sun, Y.-S. Zhang, and W. Yi, Phys. Rev. Research 3, 023022 (2021).
  • [38] S. Longhi, Phys. Rev. B 102, 201103(R) (2020).
  • [39] K. Zhang, Z. Yang, and C. Fang, Phys. Rev. Lett. 125, 126402 (2020).
  • [40] N. Okuma, K. Kawabata, K. Shiozaki, and M. Sato, M. Phys. Rev. Lett. 124, 086801 (2020).
  • [41] H. Li, H. Wu, W. Zheng, and W. Yi, Phys. Rev. Research 5, 033173 (2023)
  • [42] H.-Y. Wang, F. Song, and Z. Wang, Phys. Rev. X 14, 021011 (2024).
  • [43] S. Longhi, Phys. Rev. B 105, 245143 (2022).
  • [44] S. Longhi, and E. Pinotti, Phys. Rev. B 106, 094205 (2022).
  • [45] S. Guo, C. Dong, F. Zhang, J. Hu, and Z. Yang, Phys. Rev. A 106, L061302 (2022).
  • [46] F. Song, S. Yao, and Z. Wang, Phys. Rev. Lett. 123, 170401 (2019).
  • [47] K. Wang, T. Li, L. Xiao, Y. Han, W. Yi, and P. Xue, Phys. Rev. Lett. 127, 270602 (2021).
  • [48] T. Helbig, T. Hofmann, S. Imhof, M. Abdelghany, T. Kiessling, L. W. Molenkamp, C. H. Lee, A. Szameit, M. Greiter and R. Thomale Nat. Phys. 16, 747(2020).
  • [49] A. Ghatak, M. Brandenbourger, J. Wezel, and C. Coulais, PNAS 117(47), 29561 (2020).
  • [50] L. Palacios, S. Tchoumakov, M. Guix, I. Pagonabarraga, S. Sánchez, and A. Grushin, Nat. Commun. 12, 4691 (2021).
  • [51] X. Zhang, Y. Tian, J. Jiang, M. Lu, and Y. Chen, Nat. Commun. 12, 5377 (2021).
  • [52] D. Zou, T. Chen, W. He, J. Bao, C. Lee, H. Sun, and X. Zhang, Nat. Commun. 12, 7201 (2021).
  • [53] L. Xiao, T. Deng, K. Wang, G. Zhu, Z. Wang, W. Yi, and P. Xue, Nat. Phys. 16, 761(2020).
  • [54] L. Xiao, T. Deng, K. Wang, Z. Wang, W. Yi, and P. Xue, Phys. Rev. Lett. 126, 230402(2021).
  • [55] Z. Gu, H. Gao, H. Xue, J. Li, Z. Su, and J. Zhu, Nat. Commun. 13, 7668 (2022).
  • [56] T. Li, Y. S. Zhang, and W. Yi, Phys. Rev. B 105, 125111 (2022).
  • [57] T. Haga, M. Nakagawa, R. Hamazaki, and M. Ueda, Phys. Rev. Lett. 127, 070402 (2021).
  • [58] F. Yang, Q. Jiang, and E. Bergholtz, Phys. Rev. Research 4, 023160 (2022).
  • [59] S. Hamanaka, K. Yamamoto, and T. Yoshida, Phys. Rev. B 108, 155114 (2023).
  • [60] S. Begg, and R. Hanai, Phys. Rev. Lett. 132, 120401 (2024).
  • [61] X. Feng, and S. Chen, Phys. Rev. B 109, 014313 (2024).
  • [62] G. Lindblad, Commun. Math. Phys. 48, 119 (1976).
  • [63] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, J. Math. Phys. 17, 821 (1976).
  • [64] F. Minganti, A. Biella, N. Bartolo, and C. Ciuti, Phys. Rev. A 98, 042118(2018).
  • [65] Z. Cai and T. Barthel, Phys. Rev. Lett. 111, 150403(2013).