Operator size growth in Lindbladian SYK

The Lindblad master equation, or Lindbladian, for a density matrix in the Schrödinger picture is

$$\partial_{t}\rho=\mathcal{L}_{S}[\rho]=-i[H,\rho]+\nu\sum_{j}\left[L_{j}\rho L_{j}^{\dagger}-\frac{1}{2}\left\{L_{j}^{\dagger}L_{j},\rho\right\}\right],$$

(1)

which is interpreted as a super-operator $\mathcal{L}_{S}$. The Lindbladian with $\nu\geq 0$ can effectively describe the microscopic model for a system coupled to an environment (bath) in certain limits. On one hand, the Lindbladian can be derived in the weak coupling limit, where the system-bath coupling is weak compared to the characteristic energy scales of both the system and the bath. More precisely, the weak coupling limit encompasses the Born approximation, the Markov approximation, and the rotating wave approximation. On the other hand, the Lindbladian can also be derived in the singular coupling limit, where the bath Hamiltonian dominates over the system Hamiltonian and the system-bath interaction breuer2002theory ; lidar2019lecture ; manzano2020short . Therefore, we will not assume any hierarchy between the strengths of the two terms in the Lindbladian (1).

One can also write down the Lindbladian for an operator in the Heisenberg picture

$$\partial_{t}O=\mathcal{L}_{H}[O]=\mathcal{L}_{U}[O]+\mathcal{L}_{D}[O]=i[H,O]+\nu\sum_{j}\left[(-1)^{\eta}L_{j}^{\dagger}OL_{j}-\frac{1}{2}\left\{L_{j}^{\dagger}L_{j},O\right\}\right],$$

(2)

where $\eta=0$ when any one of $O$ and $L_{j}$ is bosonic and $\eta=1$ when both are fermionic Liu:2022god . In Liu:2022god , the Lindbladian for an operator is derived by tracing out the Markovian bath and imposing the white noise approximation on the bath correlation in an integral equation of the unitary Heisenberg equation of the operator.

The Lindbladian itself is a consistent theory for any real value of $\nu$, including $\nu<0$, mathematically. We will provide an interpretation of the case where $\nu<0$ later in our Lindbladian SYK model. In this paper, we will use the Lindbladian in (2) as our starting point and will not consider its microscopic origin in the system-plus-bath composite system.

Both Lindbladians preserve the trace and the maximally mixed state, which is proportional to the identity, namely $\mathcal{L}_{S}[\mathbb{I}]=\mathcal{L}_{H}[\mathbb{I}]=0$. When the Lindblad operators are Hermitian, $L_{j}^{\dagger}=L_{j}$, the two equations (1) and (2) are identical by mapping $H\to-H$. It is exactly this case we consider in this paper, and we will adopt the Lindbladian in the Heisenberg picture (2) throughout out this paper. The first term $\mathcal{L}_{U}[O]=i[H,O]$ describes the unitary evolution and the second term $\mathcal{L}_{D}[O]$ introduces the effect of physical non-unitary. The coefficient $\nu$ could be interpreted as the error rate of quantum gates in quantum circuits Schuster:2022bot .

In this paper, we consider an even number $N$ Majorana fermions $\left\{\psi_{j}\,|\,j=1,2,\cdots N\right\}$, with anti-commutation relation $\{\psi_{j},\psi_{k}\}=2\delta_{jk}$. The Hamiltonian $H$ is taken to be the Sachdev-Ye-Kitaev (SYK) model Kitaev:2015a ; Maldacena:2016remarks

$$H=i^{q/2}\sum_{1\leq j_{1}<\cdots<j_{q}\leq N}J_{j_{1}\cdots j_{q}}\psi_{j_{1}}\cdots\psi_{j_{q}},\quad\left\langle J_{j_{1}\cdots j_{q}}^{2}\right\rangle=\frac{J^{2}}{\binom{N-1}{q-1}}=\frac{\mathcal{J}^{2}}{2q\binom{N-1}{q-1}},$$

(3)

with even number $q\geq 2$ and $J\geq 0$. The Lindblad operators $L_{j}$ are taken as the linear jump operators Kulkarni:2021gtt

$$L_{j}=\psi_{j}/\sqrt{2},\quad j=1,2,\cdots,N.$$

(4)

The dimension of Hilbert space is $D=2^{N/2}$. To derive the Lindbladian with linear jump operators, one could couple the SYK system with a Markovia bath, where the bath correlation follows white noise approximation $\left\langle\chi_{j}(t)\chi_{k}\right\rangle=\delta(t)\delta_{jk}$, the hopping term is $H_{I}=i\sum_{j=1}^{N}\sqrt{\nu}\psi_{j}\chi_{j}$ and $\chi_{j}$ is the Majorana fermion of the bath Liu:2022god . Last but not least, the Lindbladian SYK model we consider in this paper is theoretically well-defined completely positive trace preserving map for all the parameters $\mathcal{J}>0,\,\nu\in\mathbb{R}$, although is still debatable whether one can find the original from a microscopic dynamics.

Following Qi:2018quantum ; Kulkarni:2021gtt , we will use the Choi-Jamiolkowski isomorphism with the double-copy Hilbert space $\mathcal{H}_{LR}=\mathcal{H}_{L}\otimes\mathcal{H}_{R}$ and $2N$ Majorana fermions satisfying $\left\{\psi_{j}^{a},\psi_{k}^{b}\right\}=2\delta_{jk}\delta_{ab}$ with $a,b=L,R$. Following (3), we use $\psi_{j}^{L}$ to construct $H^{L}$ and use $\psi_{j}^{R}$ to construct $H^{R}$. Given a bosonic (fermionic) operator $O$ acting on $\mathcal{H}$, we can construct two bosinic (fermionic) operators $O^{L}=O\otimes\mathbb{I},\,O^{R}=\mathbb{I}\otimes O$ ($O^{L}=O\otimes\mathbb{I},\,O^{R}=S\otimes O$) acting on $\mathcal{H}_{L}\otimes\mathcal{H}_{R}$ where Hermitian operator $S=i^{N(N-1)/2}\psi_{1}\psi_{2}\cdots\psi_{N}$. One can check that $\left[O,S\right]=0$ ($\left\{O,S\right\}=0$) for any bosonic (fermionic) linear operator $O$.

We define a maximally entangled state $\left|0\right\rangle$ in $\mathcal{H}\otimes\mathcal{H}_{R}$ as

$$\left(\psi_{j}^{L}+i\psi_{j}^{R}\right)\left|0\right\rangle=0,\ \forall j,\quad\left\langle 0|0\right\rangle=1.$$

(5)

One can check that $H^{L}\left|0\right\rangle=i^{q}H^{R}\left|0\right\rangle$. The maximally entangled state $\left|0\right\rangle$ induces a map from a linear operator acting on the single Hilbert space $\mathcal{H}$ to a state in the double-copy Hilbert space $\mathcal{H}_{LR}$ via

$$O\mapsto O^{L}\left|0\right\rangle=\left|O\right\rangle.$$

(6)

Then the identity $\mathbb{I}$ acting on $\mathcal{H}$ is mapped to $\left|0\right\rangle$. The operator trace in $\mathcal{H}$ is mapped to the inner product in $\mathcal{H}_{LR}$ via $\mathrm{Tr}[O_{1}^{\dagger}O_{2}]=D\left\langle O_{1}|O_{2}\right\rangle$.

The Lindbladian $\mathcal{L}_{H}$ is mapped to a non-Hermitian Liouvillian

$$\mathcal{L}=i\left(H^{L}-i^{q}H^{R}\right)-\nu n=i\mathcal{P}+\mathcal{X}$$

(7)

acting on $\mathcal{H}_{LR}$ via $\mathcal{L}_{H}[O]\mapsto\mathcal{L}_{H}[O]\otimes\mathbb{I}\left|0\right\rangle=\mathcal{L}\left|O\right\rangle$, where

$$n=\frac{1}{2}\sum_{j=1}^{N}\left(1+i\psi_{j}^{L}\psi_{j}^{R}\right)$$

(8)

is called the size operator, whose meaning will be explained in the next subsection. The factor $(-1)^{\eta}$ in (2) is canceled out in (7) when $O^{L}$ and $\psi_{j}^{R}$ are exchanged. At the last step of (7), we decompose the Liouvillian into anti-Hermitian part $i\mathcal{P}=i\left(H^{L}-i^{q}H^{R}\right)=iP^{\dagger}$ and Hermitian part $\mathcal{X}=-\nu n=\mathcal{X}^{\dagger}$. Thus, $\mathcal{L}^{\dagger}=-i\left(H^{L}-i^{q}H^{R}\right)-\nu n=-i\mathcal{P}+\mathcal{X}.$ One can check that $\mathcal{P}\left|0\right\rangle=\mathcal{X}\left|0\right\rangle=\mathcal{L}\left|0\right\rangle=\mathcal{L}^{\dagger}\left|0\right\rangle=0$. The Liouvillian $\mathcal{L}$ is reminiscent of the Hamiltonian for eternal traversable wormhole Gao:2016bin ; Maldacena:2018lmt ; Milekhin:2022bzx with the same second term, but here the Hamiltonians on each side have the opposite signs.

To study the information scrambling in the Lindbladian SYK model, we introduce the operator size following Qi:2018quantum . First, we define a local, complete, operator basis acting on $\mathcal{H}$, consisting of Majorana strings, namely

$$\left\{\Gamma_{I}=\Gamma_{j_{1}j_{2}\cdots j_{k}}=i^{k(k-1)/2}\psi_{j_{1}}\psi_{j_{2}}\cdots\psi_{j_{k}},\ |\ 1\leq j_{1}<j_{2}<\cdots<j_{k}\leq N\right\}$$

(9)

with notation $\left|I\right|=k$, where the factor $i^{k(k-1)/2}$ is introduced to make $\Gamma_{I}$ Hermitian. The basis is orthogonal and normalized under the inner product $\frac{1}{D}\mathrm{Tr}[\Gamma_{I}\Gamma_{J}]=\left\langle\Gamma_{I}|\Gamma_{J}\right\rangle=\delta_{IJ}$. The basis is complete since the number of elements $2^{N}$ is equal to the square of the dimension of Hilbert space $D^{2}$, even in the $N\rightarrow\infty$ limit. Second, we define the size superoperator $n[\cdots]$ so that it is diagonal on the Majorana strings basis with the eigenvalues measuring the numbers of Majorana fermions

$$n[\Gamma_{I}]=\left|I\right|,$$

(10)

which defines the size of all the linear operators acting on $\mathcal{H}$. Obviously, the identity $\mathbb{I}$ has smallest size $0$ and $\Gamma_{12\cdots N}=S$ has largest size $N$.

Based on the mapping (6) to the double-copy Hilbert space $H_{LR}$, the assignment (10) simply corresponds to the size operator defined in (8), as one can check from the eigensystem $n\left|\Gamma_{I}\right\rangle=\left|I\right|\left|\Gamma_{I}\right\rangle$ and diagonalization $n=\sum_{I}\left|\Gamma_{I}\right\rangle|I|\left\langle\Gamma_{I}\right|$. Obviously, $n\left|0\right\rangle=0,\,n\left|\Gamma_{1\cdots N}\right\rangle=N\left|\Gamma_{1\cdots N}\right\rangle$ and we denote $\left|N\right\rangle\equiv\left|\Gamma_{1\cdots N}\right\rangle$. We define the size subspace and its projection operator

$$\mathcal{H}_{n}=\text{span}\left\{\left|\Gamma_{I}\right\rangle|\left|I\right|=n,\forall I\right\},\quad\pi_{n}=\sum_{\left|I\right|=n}\left|\Gamma_{I}\right\rangle\left\langle\Gamma_{I}\right|,\quad n=0,1,\cdots,N.$$

(11)

The dimension of $\mathcal{H}_{n}$ is $\binom{N}{n}$. The size operator can be expanded as $n=\sum_{n}n\pi_{n}$. The Lindblad term $-\nu n$ in (7) with $\nu>0$ ($\nu<0$) suppresses (enhances) eigenstates with large sizes. So, we refer to the effect induced by the Lindblad term with $\nu>0$ as dissipation, and that with $\nu<0$ as production. One can measure the size of any operator $O$ acting on $\mathcal{H}$ by the size operator

$$n[O]=\frac{\left\langle O\right|n\left|O\right\rangle}{\left\langle O|O\right\rangle}=-\partial_{\mu}\ln\mathcal{G}_{\mu}[O]|_{\mu=0}$$

(12)

with $n$ the size operator in (8) and

$$\mathcal{G}_{\mu}[O]=\left\langle O\right|e^{-\mu n}\left|O\right\rangle$$

(13)

the generating function of the operator size Qi:2018quantum . The generating function also gives the size distribution via

$$\mathcal{G}_{\mu}[O]=\left\langle O|O\right\rangle\sum_{n=0}^{N}e^{-\mu n}P_{n}[O],\quad P_{n}[O]=\frac{\left\langle O\right|\pi_{n}\left|O\right\rangle}{\left\langle O|O\right\rangle},$$

(14)

where the distribution $P_{n}[O]$ is normalized as $\sum_{n}P_{n}[O]=1$.

For the sake of analytic solvability, in the main text of this paper, we focus on the case of $O=\psi_{1}(t)$ in the Heisenberg picture and calculate the size $n[\psi_{1}(t)]$ and the generating function

$$\mathcal{G}_{\mu}(t)\equiv\mathcal{G}_{\mu}[\psi_{1}(t)]=\left\langle\psi_{1}(t)\right|e^{-\mu n}\left|\psi_{1}(t)\right\rangle.$$

(15)

Obviously, $P_{n}[\psi_{1}(t)]=0$ for all even $n$. So $1\leq n[\psi_{1}(t)]\leq N-1$, in contrast to the case of a non-Markovian reservoir Zhang:2023BrownianSize . In the App. B, we consider two initial operators $e^{-\beta H/2}$ and $\psi_{1}e^{-\beta H/2}$ and numerically study their sizes under the Lindbladian evolution. For later convenience, we further introduce the double-time generating function

$$\mathcal{G}_{\mu}(t_{1},t_{2})\equiv\left\langle\psi_{1}(t_{1})\right|e^{-\mu n}\left|\psi_{1}(t_{2})\right\rangle.$$

(16)

It is the transition amplitude between the states at two different times $t_{1}$ and $t_{2}$ under the influence of the insertion $e^{-\mu n}$. Obviously, $\mathcal{G}_{\mu}(t)=\mathcal{G}_{\mu}(t,t)$.

The operator size $n[\psi_{1}(t)]$ for small enough $\nu>0$ was studied in Schuster:2022bot semi-classically, and also calculated in the same limit in Bhattacharjee:2023uwx from the Krylov complexity perspective. Here, we are able to solve this problem for any value of $\nu$ numerically at finite $q$ and analytically at large $q$, so that we can investigate the operator size and distribution dynamics in the Lindbladian SYK model at finite dissipation or production and observe fast saturation of the size and the distribution.

We will utilize the path integral to calculate the generating function of size (15). Let us consider the partition function

$$Z_{\mu}(t_{1},t_{2})=\left\langle 0\right|e^{t_{1}\mathcal{L}^{\dagger}}e^{-\mu n}e^{t_{2}\mathcal{L}}\left|0\right\rangle=\int\prod_{j}\prod_{a=L,R}{\mathcal{D}}\psi_{j}^{a}e^{-S}.$$

(17)

Since $\mathcal{L}\left|0\right\rangle=n\left|0\right\rangle=0$, we know that $Z_{\mu}(t_{1},t_{2})=1$ for any $t_{1},t_{2},\mu$. So we are free to choose any $t_{1}$ and $t_{2}$ in (17). Furthermore, the two-point function should be independent of the choice of $t_{1},t_{2}$ since $\left\langle 0\right|e^{t_{1}\mathcal{L}^{\dagger}}\psi_{j}^{a}=\left\langle 0\right|\psi_{j}^{a}$ and $\psi_{j}^{a}e^{t_{2}\mathcal{L}}\left|0\right\rangle=\psi_{j}^{a}\left|0\right\rangle$. On the r.h.s. of (17), we write the partition function as the path integral of two real Grassmann fields $\psi_{j}^{L}(\tau),\psi_{j}^{R}(\tau)$ along the double-copy contour in Fig. 1 with $t_{1},t_{2}\geq 0$ Garcia-Garcia:2022adg ; Kawabata:2022osw . The action is

	$\displaystyle-S=\int_{-t_{2}}^{t_{1}}{\rm d}\tau\Big{[}$	$\displaystyle-\frac{1}{4}\sum_{j,a}\psi_{j}^{a}(\tau)\partial_{\tau}\psi_{j}^{a}(\tau)-i{\rm sgn}(\tau)\left(H^{L}(\tau)-i^{q}H^{R}(\tau)\right)$
		$\displaystyle-\left(\nu+\mu\delta(\tau)\right)n(\tau)\Big{]},$		(18)

where $H^{L}(\tau),\,H^{R}(\tau)$ and $n(\tau)$ are functions of real Grassmann variables $\psi_{j}^{a}(\tau)$ in the forms of (3) and (8). Notice that the unusual sign function before the Hamiltonians originates from the piecewise evolution $e^{t_{1}\mathcal{L}^{\dagger}}e^{t_{2}\mathcal{L}}$ along the contour. The $\delta(\tau)$ term corresponds to the insertion of $e^{-\mu n}$ at $\tau=0$. Due to the initial and final state (5), the path integral is subjected to the boundary conditions

$$\psi_{j}^{L}(-t_{2})+i\psi_{j}^{R}(-t_{2})=0,\quad\psi_{j}^{L}(t_{1})-i\psi_{j}^{R}(t_{1})=0.$$

(19)

As in the pure SYK model Maldacena:2016remarks , we take the disorder average, keep the dominating replica diagonal part and obtain the action

$$\begin{split}-S=&\leavevmode\nobreak\ \int_{-t_{2}}^{t_{1}}{\rm d}\tau\sum_{j}\Big{[}-\frac{1}{4}\sum_{a}\psi_{j}^{a}(\tau)\partial_{\tau}\psi_{j}^{a}(\tau)-\frac{1}{2}\left(\nu+\mu\delta(\tau)\right)(i\psi_{j}^{L}(\tau)\psi_{j}^{R}(\tau)+1)\Big{]}\\ &\leavevmode\nobreak\ -\frac{i^{q}J^{2}}{2qN^{q-1}}\int_{-t_{2}}^{t_{1}}{\rm d}\tau_{1}{\rm d}\tau_{2}\sum_{j_{1}\cdots j_{N},ab}s_{ab}(\tau_{1},\tau_{2})\psi_{j_{1}}^{a}(\tau_{1})\cdots\psi_{j_{q}}^{a}(\tau_{1})\psi_{j_{1}}^{b}(\tau_{2})\cdots\psi_{j_{q}}^{b}(\tau_{2}),\end{split}$$

(20)

where $s_{ab}(\tau_{1},\tau_{2})$ comes from the disorder average between two different locations and is defined as

$$s_{ab}(\tau_{1},\tau_{2})\equiv s_{a}(\tau_{1})s_{b}(\tau_{2}),\quad s^{L}(\tau)=-i^{q}s_{R}(\tau)={\rm sgn}(\tau).$$

(21)

We introduce the bi-local fields

$$G_{ab}(\tau_{1},\tau_{2})=\frac{1}{N}\sum_{j}\psi_{j}^{a}(\tau_{1})\psi_{j}^{b}(\tau_{2}),\quad a,b=L,R,$$

(22)

and $\Sigma_{ab}(\tau_{1},\tau_{2})$ via the Lagrange multiplier method

$$1=\int{\rm D}\Sigma{\rm D}G\exp\left[-\frac{1}{2}\int_{-t_{2}}^{t_{1}}{\rm d}\tau_{1}{\rm d}\tau_{2}\sum_{ab}\Sigma_{ab}(\tau_{1},\tau_{2})\left(G_{ab}(\tau_{1},\tau_{2})-\frac{1}{N}\sum_{j}\psi_{j}^{a}(\tau_{1})\psi_{j}^{b}(\tau_{2})\right)\right].$$

(23)

By integrating out the Grassmann variables $\psi_{j}^{a}(\tau)$, we obtain the effective action for bi-local fields

	$\displaystyle-S/N=$	$\displaystyle\leavevmode\nobreak\ \frac{1}{2}\log\det\left(\frac{1}{2}\delta_{ab}\partial_{\tau}-\Sigma_{ab}\right)$
		$\displaystyle\leavevmode\nobreak\ -\frac{1}{2}\int_{-t_{2}}^{t_{1}}{\rm d}\tau_{1}{\rm d}\tau_{2}\sum_{ab}\left(\Sigma_{ab}(\tau_{1},\tau_{2})G_{ab}(\tau_{1},\tau_{2})+s_{ab}(\tau_{1},\tau_{2})\frac{J^{2}}{q}G_{ab}(\tau_{1},\tau_{2})^{q}\right)$
		$\displaystyle\leavevmode\nobreak\ -\frac{1}{4}\int_{-t_{2}}^{t_{1}}{\rm d}\tau\left(\nu+\mu\delta(\tau)\right)\left(iG_{LR}(\tau,\tau)-iG_{RL}(\tau,\tau)+2\right).$		(24)

The boundary condition (19) becomes

$$\begin{split}&G_{aL}(\tau,-t_{2})+iG_{aR}(\tau,-t_{2})=G_{La}(-t_{2},\tau)+iG_{Ra}(-t_{2},\tau)=0,\\ &G_{aL}(\tau,t_{1})-iG_{aR}(\tau,t_{1})=G_{La}(t_{1},\tau)-iG_{Ra}(t_{1},\tau)=0.\end{split}$$

(25)

In the large $N$ limit, by taking the variation on the bi-local fields, we obtain the Schwinger-Dyson (SD) equation

	$\displaystyle\frac{1}{2}\partial_{\tau_{1}}G_{ab}(\tau_{1},\tau_{2})-\int_{-t_{2}}^{t_{1}}{\rm d}\tau_{3}\sum_{c}\Sigma_{ac}(\tau_{1},\tau_{3})G_{cb}(\tau_{3},\tau_{2})=\delta_{ab}\delta(\tau_{1}-\tau_{2}),$		(26a)
	$\displaystyle\Sigma_{ab}(\tau_{1},\tau_{2})=-J^{2}s_{ab}(\tau_{1},\tau_{2})G_{ab}(\tau_{1},\tau_{2})^{q-1}-\frac{i}{2}\epsilon_{ab}\left(\nu+\mu\delta(\tau_{1})\right)\delta(\tau_{1}-\tau_{2}),$		(26b)

where $\epsilon_{LL}=\epsilon_{RR}=0$ and $\epsilon_{LR}=-\epsilon_{RL}=1$. The on-shell solution should be independent of $t_{1}$ and $t_{2}$, such that the boundary condition (25) actually holds for any $t$. This fact leads to the following simplification:

$$\begin{split}G_{LL}(\tau_{1},\tau_{2})={\rm sgn}(\tau_{21})iG_{LR}(\tau_{1},\tau_{2})={\rm sgn}(\tau_{12})iG_{RL}(\tau_{1},\tau_{2})=G_{RR}(\tau_{1},\tau_{2}),\end{split}$$

(27)

where $\tau_{12}=\tau_{1}-\tau_{2}$. To establish similar relations between the components of the self energy, we need to isolate the contribution from the anti-symmetric $\epsilon_{ab}$ term in (26b) by defining $\tilde{\Sigma}_{ab}(\tau_{1},\tau_{2})=-J^{2}s_{ab}(\tau_{1},\tau_{2})G_{ab}(\tau_{1},\tau_{2})^{q-1}$, which follows the same relations as (27), namely

$$\tilde{\Sigma}_{LL}(\tau_{1},\tau_{2})={\rm sgn}(\tau_{21})i\tilde{\Sigma}_{LR}(\tau_{1},\tau_{2})={\rm sgn}(\tau_{12})i\tilde{\Sigma}_{RL}(\tau_{1},\tau_{2})=\tilde{\Sigma}_{RR}(\tau_{1},\tau_{2}).$$

(28)

So, only one of the four components in $G_{ab}$ and $\tilde{\Sigma}_{ab}$ are independent. We choose $G_{LL}$ and $\tilde{\Sigma}_{LL}$, and simplify the SD equations as

	$\displaystyle\delta(\tau_{1}-\tau_{2})=$	$\displaystyle\leavevmode\nobreak\ \frac{1}{2}\partial_{\tau_{1}}G_{LL}(\tau_{1},\tau_{2})-2{\rm sgn}(\tau_{12})\int_{\tau_{2}}^{\tau_{1}}{\rm d}\tau_{3}\tilde{\Sigma}_{LL}(\tau_{1},\tau_{3})G_{LL}(\tau_{3},\tau_{2})$
		$\displaystyle+\frac{1}{2}{\rm sgn}(\tau_{12})\left(\nu+\mu\delta(\tau_{1})\right)G_{LL}(\tau_{12}),$		(29a)
	$\displaystyle\tilde{\Sigma}_{LL}(\tau_{1},\tau_{2})=$	$\displaystyle\leavevmode\nobreak\ -J^{2}{\rm sgn}(\tau_{1}){\rm sgn}(\tau_{2})G_{LL}(\tau_{1},\tau_{2})^{q-1},$		(29b)

such that $t_{1}$ and $t_{2}$ disappear. Now the equation is real, so we expect real solutions of $G_{LL}$ and $\tilde{\Sigma}_{LL}$. Other components could be reconstructed via (27) and (28).

The two-point function $G_{ab}(\tau_{1},\tau_{2})$ derived from the path integral, or from solving the SD equation (26) in the large $N$ limit, can be expressed as the following expectation value

$$G_{ab}(\tau_{1},\tau_{2})=\frac{1}{NZ_{\mu}(t_{1},t_{2})}\sum_{j}\left\langle 0\right|\mathcal{T}\left[e^{-\mu n(0)}\psi_{j}^{a}(\tau_{1})\psi_{j}^{b}(\tau_{2})\right]\left|0\right\rangle,$$

(30)

where the operator evolution is defined as

$$O(\tau)=\begin{cases}e^{-\tau\mathcal{L}}Oe^{\tau\mathcal{L}},&\tau\leq 0\\ e^{-\tau\mathcal{L}^{\dagger}}Oe^{\tau\mathcal{L}^{\dagger}},&\tau>0\\ \end{cases}$$

(31)

and the time ordering $\mathcal{T}$ compatible with the path integral (17) and (2.3) is

$$\mathcal{T}\left[O_{1}(\tau_{1})O_{2}(\tau_{2})\right]=\begin{cases}O_{1}(\tau_{1})O_{2}(\tau_{2}),&\tau_{1}>\tau_{2}\\ (-1)^{\eta}O_{2}(\tau_{2})O_{1}(\tau_{1}),&\tau_{1}<\tau_{2}\\ \end{cases},$$

(32)

with $\eta=1$ if both operators are fermionic, and $\eta=0$ in other cases. Under the disorder average, the generating functions (15)(16) will not depend on the choice of the Majorana index and equals the specific two-point function

$$\mathcal{G}_{\mu}(t)=Z_{\mu}(t,t)G_{LL}(t,-t),\quad\mathcal{G}_{\mu}(t_{1},t_{2})=Z_{\mu}(t_{1},t_{2})G_{LL}(t_{1},-t_{2}).$$

(33)

where $Z_{\mu}(t_{1},t_{2})=1$.

The ${\rm sgn}(\tau)$ factor in the path integral (2.3) divides the time domain $[-t_{2},t_{1}]$ into two parts $[-t_{2},0)$ and $(0,t_{1}]$. So, the double-time argument in $G_{ab}(\tau_{1},\tau_{2})$ has $4$ distinct cases. Besides the relationship (27) between the components of $G_{ab}$, we should further identify the independent time domain in $(\tau_{1},\tau_{2})$ based on other symmetries.

We will review some transformations and symmetries in the $PT$-symmetric Linbladian Prosen:2012sn ; Garcia-Garcia:2023yet . They will provide us with the symmetries of the spectrum and the symmetries of the two-point function, which connect the two-point function in different time domains.

• •

  Anti-commutation relation of fermions. Exchanging the $\psi_{j}^{a}(\tau_{1})$ and $\psi_{j}^{b}(\tau_{2})$ in (30) leads to the relation

  $$G_{ab}(\tau_{1},\tau_{2})=-G_{ba}(\tau_{2},\tau_{1}).$$

  (34)

• •

  Hermitian conjugation $\dagger$. Since $\left(e^{\tau_{1}\mathcal{L}^{\dagger}}e^{-\mu n}e^{\tau_{2}\mathcal{L}}\right)^{\dagger}=e^{\tau_{2}\mathcal{L}^{\dagger}}e^{-\mu n}e^{\tau_{1}\mathcal{L}}$, the conjugate of the two-point function will change the time argument,

  $$G_{ab}(\tau_{1},\tau_{2})=G_{ba}^{*}(-\tau_{2},-\tau_{1}).$$

  (35)

• •

  Parity $P$ and time reversal $T$. They are defined as

  $$P=e^{-\frac{\pi}{4}\sum_{j}\psi_{j}^{L}\psi_{j}^{R}}=\left(e^{-i\frac{\pi}{4}\sigma_{z}}\right)^{\otimes N},\quad TOT=O^{*},$$

  (36)

  where we have expressed the transformation acting on the representation of Majorana fermions of the Jordan–Wigner (JW) transformation,

  $$\psi^{L}_{j}=(-1)^{j-1}\sigma_{z}^{\otimes(j-1)}\otimes\sigma_{x}\otimes 1^{\otimes(N-j)},\quad\psi^{R}_{j}=(-1)^{j-1}\sigma_{z}^{\otimes(j-1)}\otimes\sigma_{y}\otimes 1^{\otimes(N-i)},$$

  (37)

  in this representation $T$ performs the complex conjugation of a matrix. The action of $P,\,T$ transformations are listed in Tab. 1. One can check that $TPT=P^{-1}$ and the Liouvillian has $PT$ symmetry

  $$PT\mathcal{L}PT=\mathcal{L}.$$

  (38)

  The $PT$ symmetry ensures that the spectrum of $\mathcal{L}$ is invariant under conjugation Bender:1998PT ; Mostafazadeh:2001jk ; Mostafazadeh:2001nr ; Mostafazadeh:2002id ; Zhang:2019gyc ; Prosen:2012sn . Combining the $PT$ transformation, we have the relation

  $$G_{ab}(\tau_{1},\tau_{2})=G_{\bar{a}\bar{b}}^{*}(\tau_{1},\tau_{2}),$$

  (39)

  where $\bar{L}=R,\,\bar{R}=L$.

  $O$	$POP^{-1}$	$TOT$	$S^{L}OS^{L}$
  $i$	$i$	$-i$	$i$
  $\psi_{j}^{L}$	$\psi_{j}^{R}$	$\psi_{j}^{L}$	$-\psi_{j}^{L}$
  $\psi_{j}^{R}$	$-\psi_{j}^{L}$	$-\psi_{j}^{R}$	$\psi_{j}^{R}$
  $n$	$n$	$n$	$N-n$
  $H^{L}$	$H^{R}$	$i^{q}H^{L}$	$H^{L}$
  $H^{R}$	$H^{L}$	$i^{q}H^{R}$	$H^{R}$

• •

  $S^{L}$ transformation. It is defined as Garcia-Garcia:2023yet

  $$S^{L}=i^{N(N-1)/2}\psi_{1}^{L}\psi_{2}^{L}\cdots\psi_{N}^{L}=\Gamma_{1\cdots N}^{L}$$

  (40)

  with $S^{L}S^{L}=1$ and $(S^{L})^{\dagger}=S^{L}$. Some results of $S^{L}$ transformation are listed in Tab. 1. The particle-hole conjugation is $TS^{L}$ Cotler:2016fpe . We can check that

  $$S^{L}(\mathcal{L}+\frac{1}{2}\nu N)S^{L}=(\mathcal{L}+\frac{1}{2}\nu N)|_{\nu\to-\nu}=-(\mathcal{L}+\frac{1}{2}\nu N)^{\dagger},\quad S^{L}\left|0\right\rangle=\left|N\right\rangle.$$

  (41)

  So the spectrum of $(\mathcal{L}+\frac{1}{2}\nu N)$ is symmetric respective to the origin. Since $S^{L}$ does not preserve the state $\left|0\right\rangle$, it can not be a new symmetry of the two-point function. However, it can relate the two-point function defined on state $\left|0\right\rangle$ (30) and the following two-point function defined on state $\left|N\right\rangle$

  $$G_{ab}(\tau_{1},\tau_{2})=\frac{1}{N\left\langle N|N\right\rangle}\sum_{j}\left\langle N\right|\mathcal{T}\left[e^{-\mu(N-n(0))}\psi_{j}^{a}(\tau_{1})\psi_{j}^{b}(\tau_{2})\right]\left|N\right\rangle\Big{|}_{\nu\to-\nu},$$

  (42)

  where $\left\langle N\right|e^{-\mu n}\left|N\right\rangle=e^{-\mu N}$ is used. The relation between two-point functions results in a simple relation between two sizes

  $$n[\psi_{1}(t)]=N-n[\Gamma_{2\cdots N}(t)]|_{\nu\to-\nu}.$$

  (43)

  Similarly, for the other $\Gamma_{I}$ operators, we will have $n[\Gamma_{I}(t)]=N-n[\Gamma_{I^{\complement}}(t)]|_{\nu\to-\nu}$, where $I^{\complement}$ is the complement of $I$. This relation aligns with our intuition: a jump term in the Linbladian with strength $|\nu|$, which suppresses the size of $\Gamma_{I}(t)$, has the same effect as the jump term with $-|\nu|$ that enhances the size of $\Gamma_{I^{\complement}}(t)$. In this sense, the Lindbladian with $\nu<0$ is also meaningful for the system coupled to a Markovian reservoir.

In summary, the spectrum of $(\mathcal{L}+\frac{1}{2}\nu N)$ is symmetric with respect to both the real axis and the imaginary axis. The two-point function has the symmetries:

$$G_{ab}(\tau_{1},\tau_{2})=-G_{ba}(\tau_{2},\tau_{1})=G_{ba}^{*}(-\tau_{2},-\tau_{1})=G_{\bar{a}\bar{b}}^{*}(\tau_{1},\tau_{2}).$$

(44)

Combining them with the simplification (27), we can reduce the two-point function to a single component in an independent time domain. We choose

$$G_{LL}(\tau_{1},\tau_{2})\,\text{ with }-\tau_{2}\geq\tau_{1}\geq\tau_{2}.$$

(45)

The time domain forms a triangle, as shown in Fig. 2. With this analysis on our hands, we will now numerically and analytically solve the problem in the following two sections.

We will numerically diagonalize the Liouvillian (7) at $q=4$ and $N=10$ and then study the evolution of the operator size and size distribution. We observe that, when $\nu>0$ ($\nu<0$), the operator size growth is suppressed (enhanced) and the size reaches a stable value smaller (bigger) than $N/2$ at the late time. We skip the case of $q=2$, since the size is trivially $n[\psi_{1}(t)]=1$.

We use the Jordan-Wigner transformation (37) to construct the Liouvillian (7) and numerically diagonalize it to find its spectrum $\left\{E_{k}|k=1,2,\cdots,D\right\}$, see also Kulkarni:2021gtt ; Zhou:2021yyw . We plot the shifted spectrum $\left\{\tilde{E}_{k}=E_{k}+\frac{1}{2}\nu N\right\}$, because it is symmetric with respect to both the real axis and the imaginary axis, according to the symmetry analysis in Subsec. 2.4. The spectrum of one realization at positive $\nu$ is shown in the left panel of Fig. 3, which is identical to the spectrum with $-\nu$. The spectrum of Liouvillian is quite different from the pure SYK spectrum Cotler:2016fpe ; Garcia-Garcia:2016mno ; you2017sachdev . When $\nu/J\gg 1$, the spectrum is real and exhibits energy bands and gaps, where the bands are approximately labeled by their sizes, as shown in the right panel of Fig. 3. When $\nu/J\sim 1$, some real energy pairs collide with each other $E_{k}\to E_{k^{\prime}}$ and then move into the complex plane. The points where such collision happens are called exceptional points Bender:1998PT . When $\nu/J\ll 1$, most but not all the energies become complex.