Scrambling Transition in Free Fermion Systems Induced by a Single Impurity

Generic interacting many-body systems can serve as their own bath, a pivotal element for the manifestation of quantum thermalization in isolated systems [1, 2]. This involves the obscuring of all local initial conditions within the entire system after prolonged evolution, measured by the growth in operator size [3, 4, 5, 6, 7, 8, 9]. In this process, interactions play a crucial role. In the absence of interactions, excitations with infinite lifetimes can carry quantum information, remaining free from dissipation [10]. This occurs because quadratic Hamiltonians conserve the number of field operators, thereby preventing the growth of operator size. However, in the presence of interactions, a single excitation can undergo scattering, giving rise to multiple excitations. This iterative process leads to a rapid increase in complexity for simple initial operators subject to the Heisenberg evolution [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32].

Recently, new insights into this problem have been gained from the study of fermionic systems that only interact through a single impurity [33]. The key question is whether a single impurity can effectively scramble the entire system. While one might typically anticipate that introducing a local impurity into a large many-body system does not lead to significant changes in bulk dynamics, it becomes apparent that even with just a single interaction, the Hamiltonian is no longer quadratic, thus permitting the growth of operator size. Indeed, previous study unveils the emergence of information scrambling in 1D systems with spatial and temporal random short-range hoppings for arbitrary weak interaction strength at the impurity [33]. (Many other papers have investigated similar setups—non-equilibrium dynamics in the presence of boundary perturbations. See Ref. [34] for a review on classical stochastic processes and see Ref. [35] for a review on quantum systems. Here, we study this from the perspective of quantum information. We also note a paper that studies information scrambling in an integrable Kondo model [36].) In this setup, both the operator and entanglement exhibit diffusive scaling, stemming from the random walk characteristics of fermion operators on the 1D lattice. However, it is acknowledged that the properties of random walks are greatly influenced by the connectivity among different sites across various geometries.

Building upon this insight, we delve deeper into investigating the growth of operator size within such systems across arbitrary dimensions. Utilizing a phenomenological percolation model on trees, we identify a notable distinction between systems in 3D and lower dimensions, summarized in Fig. 1. In 3D, interactions do not lead to a persistent growth of operator size unless their strength exceeds a critical value. Within this non-scrambling regime, the operator quickly escapes from the impurity, following the Pólya’s theorem, and we term this as the escape phase. For stronger interactions, the operators near interacting sites can scramble into non-local operators, analogous to the scenario in lower dimensions with short-range hoppings. Additionally, when the hopping range is long enough, we also identify similar dynamical transitions in lower dimensions through the utilization of Lévy flight properties. Our predictions are demonstrated by both numerical simulations based on a small $N$ Brownian circuits [28, 30, 29, 31, 32] and analytical calculations in a solvable large $N$ Brownian Sachdev-Ye-Kitaev (SYK) chain [37, 38, 39, 40, 41, 42, 43, 44]. This suggests universality regardless of the dimensionality of the local Hilbert space.

To be concrete, we first examine Brownian circuits with nearest-neighbor hoppings in a generic dimension $D$, where operator dynamics can be mapped to a classical stochastic process. As we will demonstrate, the physical picture obtained in this model is also applicable to Brownian SYK models with large Hilbert space dimensions. Our focus is on Brownian circuits of Majorana fermions. Each site hosts a single Majorana mode $\chi_{\bm{x}}$ with canonical anti-commutation relations $\{\chi_{\bm{x}},\chi_{\bm{y}}\}=2\delta_{\bm{x}\bm{y}}$. The Hamiltonian reads

where the first term denotes the hopping of the free fermion between neighboring sites and the second term is the interaction term. Here $\square$ labels four sites in a single plaquette near the origin, as illustrated in Fig. 1 with $D\geq 2$. For $D=1$, we can simply pick four contiguous sites near the origin. Independent Brownian variables $dV_{\bm{x},\bm{y}}$ and $dJ$ satisfy the Wiener process, with

In a short time interval $dt$, the unitary evolution is given by $dU=e^{-idH}$. We are interested in the operator dynamics, governed by the Heisenberg evolution $O(t+dt)=dU^{\dagger}O(t)dU$. To study the growth of operator size, we introduce a complete orthonormal basis of Hermitian operators $\{B_{\mu}\}=\{i^{q(q-1)/2}\chi_{\bm{x}_{1}}\chi_{\bm{x}_{2}}...\chi_{\bm{x}_{q}}\}$. Each Majorana string $B_{\mu}$ can be labeled by its height $\mathbf{h}_{\mu}$, defined as $h_{\mu,\bm{x}}=1$ for $\bm{x}\in\{\bm{x}_{1},\bm{x}_{2},...,\bm{x}_{q}\}$ and otherwise $h_{\mu,\bm{x}}=0$. The size of $B_{\mu}$ is further defined as $n_{\mu}=\sum_{{\bm{x}}}h_{\mu,{\bm{x}}}$

We expand $O(t)$ in this set of basis operators as $O(t)=\sum_{\mu}\alpha_{\mu}(t)B_{\mu}$, where $\alpha_{\mu}(t)$ represents the wave function for the operator evolution. In Brownian circuits, the phase of $\alpha_{\mu}(t)$ is averaged out due to the temporal randomness, and the evolution can be formulated in terms of a classical stochastic process described by the probability distribution $f_{\mu}(t)=|\alpha_{\mu}(t)|^{2}$, which is normalized $\sum_{\mu}f_{\mu}(t)=1$ due to unitarity [7, 9]. The size of $O(t)$ is defined as $N(t)\equiv\sum_{\mu}n_{\mu}f_{\mu}(t)$. By generalizing the analysis in Ref. [33], this dynamics describes a symmetric exclusion process (SEP) with a single source term at the origin [45, 46]. The governing master equation for this dynamics, along with its derivation details, are provided in the Appendix. Based on this master equation, we conduct numerical simulation using the following update rules:

1. 1.

   We implement an unbiased random walk for each particle independently.

2. 2.

   Each time when one particle returns to the origin, we branch it by adding $n_{i}$ particles with a probability $p_{i}$.

In the original model (1), we expect $n_{i}=2$ and $p_{i}\propto J/V$. The system is initialized by putting a single particle at the origin, and the operator size growth is studied by counting the number of particles.

We note that in 1D and 2D, an operator originating from the origin exhibits diffusive spreading, independent of the parameter $p_{i}$, eventually encompassing the entire space. Complexity arises in 3D, where for large $p_{i}$, the operator continues to grow over time. However, decreasing $p$ can mitigate its growth. Notably, for $p<p_{c}$, during time evolution, the operator’s size saturates to a finite constant, indicating the presence of a scrambling transition as $p_{i}$ varies.

To understand the presence or absence of the transitions across various dimensions, we recognize that the dynamics can be modeled as percolation on a tree. As shown in Fig. 2 (a), in this model, each vertex is branching into $n$ vertex, with each edge subject to removal with probability $1-p$. It is known that this model exhibits a percolation transition at $pn=1$ [47]. When $pn>1$, the root is connected to an infinite number of the vertices in the tree. Conversely, if $pn<1$, the root connects to only a finite number of vertices in the tree. Our model maps to the percolation model as follows: each edge represents a particle from the origin, with each vertex’s branching number estimated as $n=n_{i}p_{i}+1$. Edges from a vertex are removed together with probability $1-p_{r}$ if the particle doesn’t return to the origin. Here $p_{r}$ represents the return probability of a random walker to return to its starting point. The transition of this model is expected to occur at $p_{r}n=1$.

Given $n=n_{i}p_{i}+1\geq 1$, observing a scrambling transition requires $p_{r}$ to be less than 1. In 1D and 2D, where the random walker returns to its starting point with probability one, there is no scrambling transition, and the model remains in the scrambling phase as long as $p_{i}>0$. On the other hand, according to Pólya’s theorem [48, 49], the 3D random walker will return to its starting point with probability $p_{r}\approx 0.341$. A brief review of the derivation is presented in the supplementary material for clarity. The probability of returning being less than one allows us to induce a scrambling transition at some finite $p_{i}^{c}$, such that $(n_{i}p_{i}^{c}+1)p_{r}=1$. For $p_{i}<p_{i}^{c}$, particles escape from the origin, leading to saturation of the operator size. Specifically, we find $p_{i}^{c}\approx 0.646$ for $n_{i}=3$, which aligns with the results from 3D classical particle model simulations illustrated in Fig. 2 (b).

We can generalize the above discussion to fermion operators that are initially not located at the origin. As demonstrated in the supplementary material, the probability of reaching the origin is inversely proportional to the distance $r$ from the origin. This suggests that these operators have a reduced probability of reaching the origin. Nevertheless, upon returning to the origin, they may undergo scrambling if $p_{i}>p_{i}^{c}$.

We now develop an analytically solvable large-N model capable of demonstrating a scrambling transition. Specifically, we consider the Brownian SYK model [42, 43] with a single interacting impurity. In this model, we have $N$ Majorana fermions $\chi_{\bm{x},i}$ with $i\in\{1,2,...,N\}$ on each site, which satisfies $\{\chi_{\bm{x},i},\chi_{\bm{y},j}\}=2\delta_{\bm{x}\bm{y}}\delta_{ij}$. The Hamiltonian reads

where the second term represents a $q$-body interaction on the impurity with $q\geq 4$. Brownian variables with different indices are independent and satisfy

The model can be analyzed using the large-$N$ expansion. Focusing on its real-time dynamics, we first introduce the retarded Green’s functions $G^{R}_{\bm{x}}(t)\equiv-i\theta(t)\langle\{\chi_{\bm{x},i}(t),\chi_{\bm{x},i}(0)\}\rangle$, where $\theta(t)$ is the Heaviside step function. In SYK-like models, the self-energy is dominated by melon diagrams, which gives $\Sigma^{R}_{\bm{x}}(\omega)=i(zV+J\delta_{\bm{x},\bm{0}})/4\equiv i\Gamma_{\bm{x}}/4$. Here, $z=2D$ is the coordination number for the square lattice and $\Gamma_{\bm{x}}$ represents the quasi-particle lifetime on site $\bm{x}$. Transforming into the time domain, we find $G^{R}_{\bm{x}}(t)=-2ie^{-\Gamma t/2}$.

In large-$N$ systems, due to the large on-site Hilbert space, the operator size usually experiences exponential growth in the early-time regime, characterized by a rate known as the quantum Lyapunov exponent $\varkappa$. Consequently, we are examining whether the system undergoes scrambling by assessing the existence of a non-vanishing $\varkappa$. It is known that the average size is related to out-of-time-order (OTO) commutator as [9, 15]

where we introduce a subscript $\bm{x}$ to denote the location of the operator. In SYK-like models, the self-consistent equation for the OTO commutator arises from the ladder diagram [39, 50]

Here we neglect the inhomogeneous terms since it does not contribute to the asymptotic behavior. In Brownian models, this is equivalently expressed as a differential equation:

The first term describes the diffusive spreading, while the second term is a source term at $x=0$, which can potentially result in an exponential growth of $N_{x}$. This equation takes the form of the imaginary-time Schrödinger equation on the lattice with an attractive delta potential with depth $\sim J$ at the origin. Assuming $N_{\bm{x}}(t)=\exp(\varkappa t)N_{\bm{x}}$ for sufficiently long time $t$, we recognize that a positive Lyapunov exponent $\varkappa$ corresponds to a bound state with energy $-\varkappa$. The solution is given by the Lippmann-Schwinger equation

Here, $\alpha$ labels different spatial directions $\alpha\in\{x,y,...\}$.

We focus on the regime with a shallow bound state $\varkappa/V\ll 1$, where we can transition to the continuum limit by expanding $1-\cos k\approx k^{2}/2$. In this limit, Eq. (8) aligns with the standard Schrödinger equation with a quadratic dispersion. It is well-established that in 3D, a finite depth of the potential is required to sustain a bound state. In contrast, in 1D or 2D, a bound state emerges under infinitely weak attractions, leading to exponential growth in operator size. We can further make a direct connection between the bound state problem and the random walk picture in the last section. To determine the critical point $J_{*}$, we set $\varkappa=0$ on the R.H.S. of (8). Then, as reviewed in the supplementary material, it can be related to the returning probability $p_{r}$ by $\text{R.H.S.}=[2(1-p_{r})D]^{-1}$. In particular, $p_{r}=1$ for $D\leq 2$ originates from the divergence of the integral $\sim\int k^{D-1}dk/k^{2}$. The critical point $J_{*}$ is then given by $J_{*}/V=2D(1-p_{r})/(q-2)$, a close analog of results in Brownian circuits. For $D=3$, this predicts $J(q-2)/2DV=(1-p_{r})=0.659$. We can further obtain the quantum Lyapunov exponent by solving Eq. (8) exactly, using the analytical results for lattice Green’s functions [51]. The result is plotted in Fig. 3 for $D\in\{1,2,3\}$. In 3D, we can clearly observe a transition from escape phase with $\chi=0$ to scrambling phase with $\chi>0$.

We now generalize our results to systems with long-range random hoppings [30, 31, 52, 53, 52, 53, 54, 55, 56]. Since we have demonstrated that both $N=1$ Brownian circuits and the large $N$ Brownian SYK model share the same phase diagram, here we take the Brownian SYK model as an example. The hopping term in the Hamiltonian now becomes

which $dV_{\bm{x},\bm{y}}^{ij}$ still has variance (4). Carrying out the similar calculation as in the last section, we find firstly the decay rate becomes $\Gamma_{\bm{x}}=\sum_{\bm{y}\neq\bm{x}}\frac{V}{|\bm{x}-\bm{y}|^{\alpha+D}}+J\delta_{\bm{x},\bm{0}}$. Therefore, to ensure the convergence of the decay rate, we focus on $\alpha>0$. Otherwise, the model is effectively all-to-all connected, whose scrambling dynamics in the presence of a singular interaction has been analyzed in Ref. [33]. Secondly, the random walk is replaced by a Lévy flight. This is reflected in the Lippmann-Schwinger equation, which now leads to

To determine whether the R.H.S. diverges, we first perform the small $k$ expansion of the denominator, which gives

As a result, the returning probability of the Lévy flight is $1$ for $D\leq\alpha$ when $\alpha\in(0,2)$ and $D\leq 2$ when $\alpha\in[2,\infty)$. This defines the parameter regime referred to as case 1 in Fig. 1. Conversely, when the integral converges for higher dimensions, an escape-to-scrambling transition typically occurs, denoted as case 2 in Fig. 1. In particular, for $\alpha\geq 2$, the phase diagram is the same as the short-range hopping case.

In this study, we investigate the information dynamics in free fermion systems with spatial and temporal randomness that interacts through a single impurity. Our findings reveal that a solitary interaction can trigger a scrambling phase transition depending on the system’s dimensionality and hopping range. We establish a universal phase diagram by employing both $N=1$ Brownian circuits and large $N$ Brownian SYK models. In high dimensions or with a long hopping range, the model undergoes an escape-to-scrambling transition upon tuning the interaction strength $J$. Within the scrambling phase, operators proximate to the interacting site can scramble into an extensive operator under time evolution. On the other hand, in low-dimensional systems with a short hopping range, even arbitrarily weak interactions lead to the scrambling of the entire system. We expect that these problems can also be characterized by other information quantities, such as the entanglement entropy and mutual information  [33].

The dynamical transitions of information scrambling have been observed in various contexts. Specifically, Refs. [57, 58] identify an environment-induced scrambling transition in systems embedded in environments. One can also treat Majorana fermions in bulk as an environment, achieving the transition by increasing the system-environment coupling $V$. In Refs. [57, 58], the information never returns to the system once it enters the environment, corresponding to $p_{r}=0$. Our analysis in this work proposes a refinement of the picture for environments with memory. (A more recent preprint [59] studied the scrambing transition where the information backflow can be tuned. Consistent with our results, they didn’t find scrambling transition in low dimensions when the information backflow was nonzero, corresponding to $p_{r}=1$.) The discussions in our paper are also pertinent to the intriguing question of whether a single thermal island can thermalize the entire system, a query crucial for understanding the existence of many-body localization phases. To address this question, it is imperative to extend the current discussions to models with static hopping strength, where localization is feasible. Nevertheless, the theoretical analysis in such cases becomes considerably more challenging, and we defer this task to future works.