Entanglement growth after inhomogenous quenches

TL;DR

This paper investigates how initial inhomogeneities in conserved quantities influence entanglement growth in quantum systems, revealing a diffusive $ oot t$ scaling of entanglement features through an effective model and numerical simulations.

Contribution

It introduces a generalized minimal cut model incorporating local entropy density and connects entanglement dynamics with KPZ and diffusion equations, supported by numerical evidence.

Findings

Entanglement features grow as $ oot t$ due to initial inhomogeneities.

Diffusive charge transport constrains entanglement growth to $ oot t$.

Features persist even after charge density equilibrates.

Abstract

We study the growth of entanglement in quantum systems with a conserved quantity exhibiting diffusive transport, focusing on how initial inhomogeneities are imprinted on the entropy. We propose a simple effective model, which generalizes the minimal cut picture of \textit{Jonay et al.} in such a way that the `line tension' of the cut depends on the local entropy density. In the case of noisy dynamics, this is described by a Kardar-Parisi-Zhang (KPZ) equation coupled to a diffusing field. We investigate the resulting dynamics and find that initial inhomogeneities of the conserved charge give rise to features in the entanglement profile, whose width and height both grow in time as $\propto\sqrt{t}$. In particular, for a domain wall quench, diffusion restricts entanglement growth to be $S_\text{vN} \lesssim \sqrt{t}$. We find that for charge density wave initial states, these features in the entanglement profile are present even after the charge density has equilibrated. Our conclusions are supported by numerical results on random circuits and deterministic spin chains.

Equations

$$S_1(x,t+1)\le S_1(x-1,t)+s_{x,t+1},$$

$$S_1(x,t+1)$$

$$\displaystyle\leq\min{\big{(}}S_{1}(x-1,t)+s_{x,t+1},S_{1}(x+1,t)+s_{x+1,t+1}{\big{)}}.$$

$$S_1=\sum _a{p}_a{S}_{1,a}-\sum _a{p}_a\ln p_a,$$

$$S$$

$$\le \min \left(S_{1,a}(x-1,t)+\ln {\eta }_a,S_{1,a}(x+1,t)+\ln {\eta }_{Q-a}\right),$$

$$\chi (x,t+1)=\sum _a\min \left(\chi (x-1,t){\eta }_a,\chi (x+1,t){\eta }_{Q-a}\right).$$

$$S_0(x,t+1)=\min \left(S_0(x-1,t),S_0(x+1,t)\right)+f_N(Q,\Delta ),$$

$$f(\mathfrak{q},\delta )\equiv f_{N\to \infty }(2N\mathfrak{q},N\delta )=\{\begin{array}{ll}\frac{h_2({\xi }_{\ast })}{\ln 2}& \text{ if }{\xi }_{\ast }\le 1/2\\ 1& \text{ if }{\xi }_{\ast }>1/2.\end{array}$$

$$S(x,t+1)=\min \left(S(x-1,t),S(x+1,t)\right)+h_2(\mathfrak{q}).$$

$${\partial }_{t}S(x,t)$$

$${\partial }_t\mathfrak{q}(x,t)$$

$$S_1(x,t)\le \sum _{x^{\prime }\le x}s(x^{\prime },t)=\sqrt{Dt}f\left(\frac{x}{\sqrt{Dt}}\right),$$

$$U_{\text{F}}$$

$$H_1$$

$$H_3$$

$$H_2$$

$${\rho }_A=\sum _{Q_A}{p}_{Q_A}{\rho }_A^{(Q_A)},$$

$$S_1[{\rho }_A]$$

$$=S_Q+S_{\text{conf}}.$$

$$S_2^{(a)}(t)\equiv -\log \overline{\mathrm{Tr}({\rho }_A(t{)}^2)}=-\log \overline{\mathcal{P}[{\rho }_A(t)]}.$$

$$\begin{gathered} \overline{U^{\ast }\otimes U\otimes U^{\ast }\otimes U}=\sum _{Q_1,Q_2}\frac{1}{d_{Q_1}{d}_{Q_2}}|{\mathcal{I}}_{Q_1{Q}_2}⟩⟨{\mathcal{I}}_{Q_1{Q}_2}|+ \\ \sum _{Q_1,Q_2}\frac{1}{d_{Q_1}{d}_{Q_2}-{\delta }_{Q_1{Q}_2}}|{\mathcal{J}}_{Q_1{Q}_2}⟩⟨{\mathcal{J}}_{Q_1{Q}_2}|, \end{gathered}$$

$$|{\mathcal{I}}_{Q_1{Q}_2}⟩$$

$$|{\mathcal{J}}_{Q_1{Q}_2}⟩$$

$$\mathcal{P}[{\rho }_A]={\text{tr}}_A({\rho }_A^2)\equiv \text{tr}({\mathcal{S}}_A[\rho \otimes \rho ]),$$

$${\mathcal{S}}_A=\underset{x\in A}{⨂}|\mathcal{I}{⟩}_x\underset{x^{\prime }\notin A}{⨂}|\mathcal{K}{⟩}_{x^{\prime }},$$

$$|\mathcal{I}⟩\equiv \sum _{\alpha ,\beta \in {\mathcal{H}}^1}|\alpha \alpha \beta \beta ⟩$$

$$\sum _{Q_1,Q_2}\frac{⟨{\mathcal{I}}_{Q_1{Q}_2}|\mathcal{I}\mathcal{K}⟩}{d_{Q_1}{d}_{Q_2}}|{\mathcal{I}}_{Q_1{Q}_2}⟩$$

$$\sum _{Q_1,Q_2}\frac{⟨{\mathcal{J}}_{Q_1{Q}_2}|\mathcal{S}\mathcal{I}⟩}{d_{Q_1}{d}_{Q_2}}|{\mathcal{I}}_{Q_1{Q}_2}⟩$$

$$f_{Q_1{Q}_2}\equiv \sum _{a,b_1,b_2=0}^N{\delta }_{a+b_1=Q_1}{\delta }_{a+b_2=Q_2}\frac{{\eta }_a{\eta }_{b_1}{\eta }_{b_2}}{\eta },$$

$$P_Q\rho ={\delta }_{Q\bar{Q}}\rho ,$$

$$\frac{\eta f_{\bar{Q}\bar{Q}}}{d_{\bar{Q}}(d_{\bar{Q}}+1)}\left(|\mathcal{II}⟩+|\mathcal{KK}⟩\right).$$

$$\mathcal{P}(x,t+1)=e^{s_2(\bar{Q})}\left(\mathcal{P}(x-1,t)+\mathcal{P}(x+1,t)\right).$$

$$S_2^{(a)}(x,t+1)=\text{min}(S_2^{(a)}(x-1,t),S_2^{(a)}(x+1,t))+s_2(\bar{Q}).$$

Full text

I Introduction

Understanding how, and under what conditions, closed quantum systems approach thermal equilibrium due to their own unitary dynamics has been at the center of much recent attention, both theoretically and experimentally Rigol et al. (2008); Calabrese and Cardy (2006); D’Alessio et al. (2016); Gogolin and Eisert (2016); Kaufman et al. (2016). One important insight is to focus on the reduced density matrices of sufficiently small subsystems, which relax to a thermal Gibbs state with temperature, chemical potential, etc. set by the initial conditions, in systems where the Eigenstate Thermalization Hypothesis is satisfied Deutsch (1991); Srednicki (1994); D’Alessio et al. (2016). These long-time states have a von Neumann entropy that is extensive in subsystem size, in accordance with the prediction of thermodynamics. For a closed system, this entropy comes entirely from entanglement between the subsystem and its environment - in this sense we can say that a thermalizing quantum system acts ‘as its own bath’. Therefore an important part of understanding the mechanisms of thermalization is to describe how entanglement between subsystems builds up as a consequence of unitary dynamics, a question that has recently become partially amenable to experimental probes in systems of cold atoms, through the measurement of so-called Rényi entropies Abanin and Demler (2012); Daley et al. (2012); Islam et al. (2015); Kaufman et al. (2016); Elben et al. (2018).

A simple picture of entanglement growth was developed for non-interacting and critical one-dimensional systems Calabrese and Cardy (2005, 2007) and then generalized to generic quantum integrable systems Alba and Calabrese (2017, 2018). In this picture, entanglement is carried by pairs of quasi-particles with opposite momenta that are present in the initial state, leading to linear growth of entanglement entropies for a contiguous subsystem with eventual saturation to a volume law. This description crucially relies on the notion of infinitely long-lived quasi-particle excitations. However, linear growth of entanglement is in fact more generic and has been shown to apply quite generally to so-called ‘global quenches’ Calabrese and Cardy (2006), even in systems which do not admit a quasi-particle description, and exhibit diffusive transport Kim and Huse (2013). More recently, important advances have been made in explaining this linear growth from a a coarse-grained hydrodynamic or ‘minimal cut’ picture Nahum et al. (2017); Jonay et al. (2018). This picture is supported by calculations in random circuits Nahum et al. (2018a); Zhou and Nahum (2018), holography Mezei (2018), and results in spin chains, both numerical Jonay et al. (2018) and analytical Bertini et al. (2018a). As an aside, however, we note that not all measures of entropy grow ballistically or are directly amenable to the minimal cut picture, even in ergodic systems; in particular a recent work by the present authors Rakovszky et al. (2019) argued that in systems with diffusive transport, such as the ones studied below, higher Rényi entropies grow diffusively after a global quench.

Most of these results concern states that are invariant under spatial translations (hence the qualifier ‘global’). It is interesting to ask how the above picture of linearly growing entanglement entropy generalizes to the case of initial states with long-wavelength inhomogeneities. These are especially important if they involve changes in the densities of conserved quantities - these inhomogeneities than have to smooth out during time evolution through transport processes, implying that equilibrium can only be approached on time scales that are long compared to the associated wavelengths that characterize the initial state. An extreme case is that of a domain wall: an initial state where the density of one or more conserved quantity jumps between two extremal values in the middle of the system. Entanglement growth starting from such a domain wall was investigated for integrable systems, where both logarithmic Eisler et al. (2009); Alba and Heidrich-Meisner (2014); Vidmar et al. (2017) and power law Ljubotina et al. (2017); B. Bulchandani and Karrasch (2018) growth have been found for the entanglement. More recently, the quasi-particle picture has been extended to a variety of inhomogeneous initial states, using notions of generalized hydrodynamics, applicable to integrable systems Bertini et al. (2018b); Alba et al. (2019); Mestyán and Alba (2019). However, the case of generic interacting systems remains largely unexplored.

The goal of the present paper is to investigate the question of entanglement production, starting from inhomogenous initial states in generic non-integrable systems, and elucidate how it is affected by the diffusive transport that is expected to be generic for the high temperature regime of non-integrable lattice systems Bloembergen (1949); Gennes (1958); Kadanoff and Martin (1963); Bohrdt et al. (2017). We approach this problem by combining the study of simple toy models with general considerations arising from hydrodynamics. In particular we consider local random unitary circuits with a U(1) symmetry, which have been used recently as a minimal model for dynamics in ergodic quantum many-body systems both with Khemani et al. (2018); Rakovszky et al. (2018, 2019) and without Nahum et al. (2018a); von Keyserlingk et al. (2018) conserved quantities. Our approach here can be considered as a generalization of the hydrodynamic picture of Ref. Jonay et al., 2018, by considering how the entanglement couples to the simpler, diffusive hydrodynamic equations obeyed by the conserved densities such as particle number, energy etc.

Using simple considerations, we conjecture an effective model that describes the coupling between entanglement growth and transport in our noisy random circuit setup, which can be interpreted as a Kardar-Parisi-Zhang (KPZ) equation Kardar et al. (1986), with a growth rate that is coupled to the diffusive variable. This implies that inhomogeneities in the original charge get imprinted into the growing ‘height’ profile of the entanglement. In particular, we show that whenever there is a domain wall between two regions with different charge densities, the entanglement entropy near the interface has a component that grows in time as $\sqrt{t}$. The size of the region where this scaling applies also grows in time as $\sqrt{t}$ as the charge density smooths out due to diffusion. This is sketched in Fig. 1. For a charge density wave initial state this process leads to an inhomogeneous entanglement profile across the system with a periodicity that is half the wavelength that characterizes the initial state. At longer times, as the charge distribution flattens, the entanglement profile also smooths out exponentially, but at a timescale that is parametrically larger than the time needed for the charge to equilibrate. We illustrate this general picture on a variety of different initial states, both for a random circuit model and for a deterministic system.

We also briefly discuss how the effects of diffusive charge transport show up in the von Neumann entropy even for states without large scale inhomogeneities. In particular, we argue that the number entropy, which is the component of the von Neumann entropy associated with the probability distribution of the conserved charge in a subsystem, has a characteristic growth, $\propto\log{t^{1/4}}$, which we associate to the diffusive spreading of correlations. This quantity is itself measurable in cold-atomic experiments Lukin et al. (2018).

The remainder of the paper is organized as follows. In Sec. II we introduce a charge-conserving local random circuit model that we use in the subsequent discussion. In Sec. III we analyze the effect of a single unitary gate to motivate a simple local update rule that couples entanglement growth with charge transport. This update rule is used to build an effective random surface growth model for entanglement dynamics, whose behavior we explore for a variety of inhomogeneous initial states. The results of this surface growth model are compared to numerics on both the random circuit model and on a deterministic spin chain in Sec. IV. In Sec. V we discuss the dynamics of the number entropy, which applies even in the case of initial states without long-range inhomogeneities. We conclude in Sec. VI. In App. A we review the mapping of the calculation of the annealed average second Rényi entropy in the random circuit to a classical parition function and use it to arrive at a simple update rule in the absence of charge fluctuations. In App. B we show further numerical data, complementing the results of Sec. IV.

II Local random circuit with charge conservation

To investigate the relationship between entanglement growth and transport of conserved quantities, we turn to a simple minimal model that possesses all the main ingredients: a time evolution which is unitary, local in space, and has a conserved charge that obeys diffusive dynamics. This model was originally introduced in Ref. Rakovszky et al., 2018. Its basic building block is a collection of $N$ qubits (spin-$\frac{1}{2}$ degrees of freedom), which we refer to as a ‘cell’ in analogy with the notion of a ‘fluid cell’ in hydrodynamics or as a ‘site’. The physical system consists of $L$ such cells, arranged in a one-dimensional chain. We choose a specific basis for each qubit, and refer to these two basis states, $|{{0}}\rangle$ and $|{{1}}\rangle$, as ‘empty’ and ‘filled’ - the total number of filled sites is going to be our conserved quantity and we will refer to it as ‘charge’. On each cell, there are $\eta_{a}\equiv\binom{N}{a}$ states with charge $a=0,\ldots,N$. We evolve the system in discrete steps, by applying local unitary gates acting on pairs of neighboring cells in the chain, as indicated in Fig 2a. Each of these two-site gates is chosen such that it conserves the total charge on the two sites, i.e., it is block diagonal in the charge basis with a block of size $d_{Q}=\binom{2N}{Q}$ in the charge $Q$ sector. All blocks are independent and chosen Haar randomly.

We consider two different circuit geometries. For the discussion in Sec. III, where we use the circuit to motivate a random surface growth picture of entanglement growth, we use a model where the two-site gates are applied on a randomly picked bond in each step, resulting in the irregular circuit shown in Fig. 2(b). Later on, in Sec. IV.1, we will use a regular ‘brick wall’ geometry, wherein odd numbered layers act on all the odd bonds of the chain while even numbered layers act on even bonds, as illustrated in Fig. 2(c). All the different gates are independently chosen from the charge conserving (i.e., block diagonal) random ensemble defined above. We denote averages over the different circuit realizations by $\overline{(\ldots)}$.

Such random circuit models have several advantages compared to usual Hamiltonian or Floquet systems, and have been used extensively to investigate a number of questions related to far-from-equilibium many-body systems Nahum et al. (2017, 2018a); von Keyserlingk et al. (2018); Nahum et al. (2018b); Zhou and Nahum (2018); Khemani et al. (2018); Rakovszky et al. (2018); Chan et al. (2018a, b, 2019). On the one hand they often admit exact analytical solutions, at least in certain limits (we will use the large-$N$ limit to simplify certain calculations). On the other hand, even when such exact treatment is not possible, by averaging over the randomness one can map the calculation of many relevant quantities to an effective ‘classical’ model, which is often more tractable numerically than the original quantum problem. We will take advantage of such a mapping, which we outline in App. A, in Sec. IV.1. Last but not least, even when there is no direct computational advantage (e.g., in calculations of the number entropy in Sec. V), by averaging over random time evolutions, one can get access to the relevant universal behavior more easily then in any specific microscopic model - for example the random circuit shows diffusion of conserved quantities at all timescales.

III Surface growth model of entanglement growth

Our goal is to study the effect of an inhomogeneous charge distribution on the evolution of bipartite entanglement. We first discuss how the entanglement entropy $S(x)$ across the bond $x,x+1$ evolves under the application of a single charge-conserving unitary gate of the type introduced in the previous section. Based on these considerations, we conjecture a simple random surface growth model to capture certain universal features of entanglement growth. In Sec. IV we will compare the results of this surface growth model with numerics on various spin-$\frac{1}{2}$ chains.

III.1 Effect of a single gate

We begin by considering the effect of a single charge-conserving 2-site gate, and use it to conjecture an effective update rule for the half-chain entanglement across the bond at these two sites. The main feature of this update rule, established using the subadditivity of entropy, is that the growth of entanglement is limited by the local entropy density, corresponding to the density of conserved charge. We also explore the possible interplay between charge density and the gradient of the half-chain entanglement.

We mainly focus on the half-chain von Neumann entropy, $S_{1}(x,t)\equiv-\text{tr}\left(\rho\ln{\rho}\right)$, where $\rho$ is the reduced density matrix of the subsystem consisting of all sites $\leq x$, at time $t$. The von Neumann entropy is the $\alpha\to 1$ limit of the Rényi entropy, $S_{\alpha}\equiv\frac{1}{1-\alpha}\ln{\text{tr}\left(\rho^{\alpha}\right)}$. We mostly focus on $S_{1}$ in this paper, which captures (the logarithm of) the typical eigenvalue of $\rho$, while later in this section we will briefly consider the case of $\alpha=0$, also known as the Hartley entropy, which measures (the logarithm of) the rank of $\rho$. For this quantity, we are able to establish an exact update rule in the case when a charge-conserving gate is applied to a pair of neighboring sites with a fixed amount of total charge. Later on, in Sec. IV.1 we will also consider $S_{2}$, which is more amenable to numerical calculations in the random circuit setting.

III.1.1 Von Neumann entropy $S_{1}$

We begin by noting that the growth of $S_{1}$ is constrained by the local entropy density, $s_{x,t}$ (the von Neumann entropy of the subsystem consisting of the site $x$ alone), through the subadditivity condition Nahum et al. (2018b),

[TABLE]

where we used the fact that the entanglement across the bond $x-1,x$ is unchanged by the 2-site gate acting on sites $x,x+1$ to replace $S_{1}(x-1,t+1)$ with $S_{1}(x-1,t)$. The same argument holds when replacing $S_{1}(x-1,t)$ with $S_{1}(x+1,t)$ and $s_{x,t+1}$ with $s_{x+1,t+1}$, such that

[TABLE]

What we will conjecture in the following is that replacing the inequality with an equality in Eq. (III.1.1) provides a good qualitative description of the time evolution of the von Neumann entropy, particularly in situations where the charge has locally equilibrated. After local equilibration, we also expect that $s_{x,t+1}\approx s_{x+1,t+1}\approx(x,t)$, so that it can be pulled out of the $\min$ function. The resulting equation is already sufficient to capture the qualitative features we discuss below.

The bound in Eq. (III.1.1) can be improved in certain circumstances, through the following line of reasoning, although there does not appear to be a sharp qualitative change in the physics in considering these subtler effects. In particular, we can consider a case where the two sites on which our unitary gate acts have a fixed total charge. This is motivated by the fact that if the charge distribution on the two sites in question is tightly peaked (as we would expect it to be at times in excess of the local equilibration time, in the large-$N$ limit), then the von Neumann entropy, which measures the size of typical Schmidt values, should be well approximated by focusing on Schmidt values corresponding to the average charge 111Note that this is in contrast with Rényi entropies $S_{\alpha>1}$, which are strongly influenced by fluctuations of the conserved charge, even for homogenous quenches, and can grow sub-ballistically as a consequence Rakovszky et al. (2019)

Thus, for the time being let us assume that the charge on sites $x,x+1$ takes a definite value $Q_{x,x+1}=Q$. In this case, the reduced density matrix on sites $\leq x$ is block diagonal in $Q_{x}$, the charge on site $x$. Let us denote the block with $Q_{x}=a$ as $p_{a}\rho_{a}$, where $p_{a}$ is the probability of having charge $a$ on site $x$, and $\text{tr}(\rho_{a})=1$. Then the von Neumann entropy can be written Lukin et al. (2018)

[TABLE]

where $S_{1,a}=-\text{tr}(\rho_{a}\ln{\rho_{a}})$ is the von Neumann entropy of the density matrix associated to charge $a$. Note that the last term is upper bounded by $\ln(Q+1)$, since $a$ takes values between [math] and $Q$, and therefore should be strongly subleading, since the total von Neumann entropy will eventually increase to an extensive value. This means that to a good approximation, the entropy can be written as a weighted sum over entropies of different charge sectors. One can then apply subadditivity to each block separately

[TABLE]

where $\eta_{a}\equiv\binom{N}{a}$ is the dimension of the on-site Hilbert space with charge $a$. Here the entropies of the neighboring bonds also depend on $a$, since the projection to a given charge sector can in principle change the entanglement even away from the bond $x,x+1$. However, based on a simple numerical experiment involving a 4-site random MPS, we expect that if there is a gradient $S_{1}(x+1,t)-S_{1}(x-1,t)$ in the original state, then there should also be a a gradient $S_{1,a}(x+1,t)-S_{1,a}(x-1,t)$ in the projected state as well – albeit one whose magnitude decreases with $a$. Therefore, one ends up with an upper bound that depends on a particular combination of the charge density and the entropy gradient. The same numerical experiment indicates that this bound is tighter then the bound in Eq. (III.1.1), at least away from half filling. This suggests that there is a possibility of a more complicated coupling between charge and entanglement, where the growth rate of the latter depends not only on the charge density, but also on the local gradient of $S_{1}(x,t)$. We will argue more rigorously for such a coupling in the case of $S_{0}$ below. However, our numerical experiments suggest that the main qualitative features are already captured by the simpler bound in Eq. (III.1.1), ignoring these subtler effects.

III.1.2 Parameter counting argument for $S_{0}$

The above argument, regarding the coupling between charge and entanglement gradient, can be made more precise if we consider the Hartley entropy, $S_{\alpha=0}$, which equals the logarithm of the number of non-zero Schmidt values. Below we derive an update rule for this quantity under the effect of a generic unitary acting on a generic state with well defined charge on $x,x+1$. While the assumption of fixed local charge is not well motivated in the case of $S_{0}$, it is nevertheless useful to consider, as it sheds some light on the nature of the coupling between charge and entanglement gradient mentioned above.

Consider a matrix product state (MPS) with fixed 2-site charge $Q$ on sites $x,x+1$, and a bond dimension $\chi(x)\equiv e^{S_{0}(x)}$ across this bond. Note that $\chi(x)$ is nothing else but the number of non-zero Schmidt values for a decomposition between the two halves of the chain separated by this bond. One can then generalize the parameter counting argument of Ref. Nahum et al., 2017 to estimate the bond dimension $\tilde{\chi}(x)$ across the same bond after a charge-conserving unitary has been applied to these two sites 222While this parameter counting argument is not rigorous, in the sense that fine-tuned states or unitaries could violate it, it is expected to hold for generic states, which we indeed find numerically for random MPS.. Since $Q$ is conserved, we can label the Schmidt values at $x,x+1$ (after applying the unitary) according to the amount of charge on site $x$. Let us denote this charge by $Q_{x}=a$. We will estimate the number of non-zero Schmidt values $\tilde{\chi}_{a}(x)$ in each $a$ sector separately, by equating the number of parameters in the new MPS tensors on sites $x,x+1$ with the number of equations defined by equating the new MPS to the time evolved one.

We want to find new MPS tensors on both $x$ and $x+1$ that describe the state after it has been evolved with the unitary gate on these sites. For a given $a$, the two MPS tensors have $\chi(x-1)\tilde{\chi}_{a}(x)\eta_{a}+\tilde{\chi}_{a}(x)\chi(x+1)\eta_{Q-a}-\tilde{\chi}_{a}(x)^{2}$ parameters in total, where $\chi(x\pm 1)$ are the bond dimensions on neighboring bonds and the last term takes into account the gauge freedom in choosing the MPS tensors. The number of parameters has to be sufficiently large to satisfy all the $\chi(x-1)\chi(x+1)\eta_{a}\eta_{Q-a}$ equations, coming from equating the new tensors with the original MPS evolved by a single 2-site gate. Equating the number of equations with the number of parameters gives a quadratic equation for $\tilde{\chi}_{a}(x)$. Taking the smaller of the two solutions and adding up the Schmidt values for different $a$ sectors gives the total bond dimension across a bond $x,x+1$ with charge $Q$, after applying the gate:

[TABLE]

i.e., we end up with the maximal bond dimension allowed by subadditivity, in each charge sector separately.

Without loss of generality, let $\chi(x+1)=\chi(x-1)e^{2\Delta}$, for some entropy gradient $\Delta\geq 0$ going across the bond (note that $e^{2\Delta}\leq d_{Q}$, due to subadditivity of entanglement). Then, we get the update rule for the Hartley entropy, $S_{0}(x,t)\equiv\log{\chi(x,t)}$, as

[TABLE]

where $f_{N}(Q,\Delta)\equiv\log\left(\sum_{a=0}^{Q}\min(\eta_{a},\eta_{Q-a}e^{2\Delta})\right)$, and we measure logarithms base $2^{N}$. Therefore we find a result that has more structure then the simple subadditivity bound in Eq. (III.1.1): the entropy increase depends not only on the local entropy, but also on the spatial derivative of $S_{0}(x,t)$ across the bond. This arises due to the fact that after fixing $Q_{x,x+1}$, we can count the entropy for each charge sector separately: having a gradient $\Delta$ across the bond allows for putting more charge on one site in order to increase the total entropy.

To get a sense of the behavior of the term $f_{N}(Q,\Delta)$ in Eq. (4), we take its large-$N$ limit. In this case $\eta_{a}={N\choose a}\to e^{Nh_{2}(a/N)}$, using the binary entropy function, $h_{2}(\xi)\equiv-\xi\ln{\xi}-(1-\xi)\ln(1-\xi)$. Defining rescaled variables, $\xi\equiv a/N$, $\mathfrak{q}\equiv Q/2N$ and $\delta=\Delta/N$, the condition $\eta_{a}=\eta_{Q-a}e^{2\Delta}$ translates to $h_{2}(\xi)=h_{2}(2\mathfrak{q}-\xi)+2\delta$. Let $\xi_{*}(\mathfrak{q},\delta)$ denote the solution to this equation when it exists, and $\xi_{*}=2\mathfrak{q}$ otherwise. The sum over $a$ in the definition of $f_{N}$ can then be replaced by a pair of integrals over $\xi$, on intervals $[0,\xi_{*}]$ and $[\xi_{*},2\mathfrak{q}]$. In the limit $N\to\infty$ these integrals can be approximated by the maximal value iside the interval, giving the result

[TABLE]

Fig. 3 shows the resulting function $f(\mathfrak{q},\delta)$ as a function of the charge density $\mathfrak{q}$ for a variety of $\delta$. Its most notable feature is that for $\delta>0$ it becomes constant for sufficiently large values of $\mathfrak{q}$.

The arguments above, in particular the result (4), indicate a complicated coupling between charge and entanglement entropy, where the local growth rate of the bipartite entanglement depends on a particular combination of local charge density and the entropy gradient. However, in the physical situations we consider in the next section, we find that this dependence on the entropy gradient does not change the qualitative features of entanglement growth, and in practice one can replace the constant term in the RHS with just the on-site entropy $h_{2}(\mathfrak{q})$, which is the original subadditivity bound in Eq. (III.1.1). As we argue below, in Sec. III.2, this can be understood from the fact that in a coarse-grained continuum description the difference only shows up in the form of higher order subleading corrections. Nevertheless, it would be interesting to understand whether there are physical situations where such a coupling can still play an important role.

III.2 Surface growth model

The equation of motion in Eq. (4) is a direct generalization of the local update rule for a non-symmetric random gate, derived in Ref. Nahum et al., 2017; the only difference being that the constant growth term on the RHS now depends on space and time, through the local charge, and also on the entropy gradient. We can then construct an effective, ‘surface growth’ model of entanglement growth, along the same lines as in Ref. Nahum et al., 2017: in each time step we apply the update rule on a randomly chosen bond, using the function $f(\mathfrak{q},\delta)$ derived in Eq. (5), to update the half-chain entanglement $S(x,t)$ (the ‘height’). At the same time we update the local charge densities $\mathfrak{q}(x)\equiv Q_{x}/N$ as $\mathfrak{q}(x),\mathfrak{q}(x+1)\to(\mathfrak{q}(x)+\mathfrak{q}(x+1))/2$. This leads to a coupled stochastic evolution between $\mathfrak{q}(x)$ and the entropy $S(x,t)$. Below, we investigate the behavior of this stochastic model, which we compare to simulation of quantum systems in Sec. IV.

III.2.1 General considerations

As stated before, we find numerically that the coupling between the charge and the entropy gradient does not affect the main qualitative features of the evolution, therefore in practice one can replace $f(\mathfrak{q},\delta)$ with $f(\mathfrak{q},0)=h_{2}(\mathfrak{q})$, the local binary entropy associated to the charge density. This leads to a simplified update rule

[TABLE]

In the continuum limit this corresponds to a KPZ equation for the entanglement (as in Ref. Nahum et al., 2017), but one that is coupled to the diffusion equation for the charge as

[TABLE]

where $s(x,t)\equiv h_{2}(\mathfrak{q}(x,t))$ is the entropy density, $\zeta(x,t)$ is uncorrelated white noise, $\zeta_{\mathfrak{q}}(x,t)$ is a noise field consistent with charge conservation (model B dynamics Täuber (2007)) and the diffusion constant $D$ is $1/2$ for the random circuit model Khemani et al. (2018); Rakovszky et al. (2018). Note that is we kept the function $f(\mathfrak{q},\delta)$ in Eq. (6) (instead of replacing it with $f(\mathfrak{q},0)=h_{2}(\mathfrak{q})$), the highest order terms we would need to add to the continuum description would be of the form $(\partial_{x}\mathfrak{q})^{2}(\partial_{x}S)^{2}$, which are subleading in the long wavelength limit, therefore we are justified in dropping them.

The above random surface growth model also admits an interpretation as a directed polymer problem, in the spirit of Ref. Jonay et al., 2018. As shown there, on the longest length and time scales, the entanglement resulting from Eq. (6) can be rewritten as the energy of a minimal energy polymer, characterized by an ‘entanglement line tension’ $\varepsilon(v)$, with $v$ being the slope of the polymer, representing a space-tim cut through the unitary circuit. Since our update rule differs from the one in Ref. Nahum et al., 2017 by having a constant term $s(x,t)$, rather then $1$, on the RHS, one has to rescale $S\to S/s$ to get the same physics locally. This implies that the line tension gets rescaled as $\varepsilon(v)\to s(x,t)\varepsilon(v)$. Therefore the calculation of entanglement growth becomes that of finding a polymer with minimal energy in a space-time dependent background, where the background itself contains a deterministic evolving part, governed by the diffusion equation, as well as random noise. In particular, parts of the systems with very low/high fillings act as bottlenecks for the entanglement growth of nearby regions 333This is similar to the mechanism proposed for entanglement growth in disordered Griffith phases where the bottleneck is provided by localized regions which act as ‘weak links’, see Ref. Nahum et al., 2018b. On the other hand, the charge density itself undergoes diffusion, tending toward a more homogenous distribution. The entropy density follows the charge distribution, leading to a speed-up of entanglement growth in regions where it gets closer to half filling and vice versa.

We expect this generalized minimal cut picture to apply also to systems without noise, much like its original version which did not take conservation laws into account Jonay et al. (2018). In our case this would correspond to an equation of motion analogous to Eq. (III.2.1), but with the noise terms omitted. The original minimal cut picture can also be generalized to higher spatial dimensions, replacing the polymer with a ‘membrane’ Nahum et al. (2017); Jonay et al. (2018); Mezei (2018). We expect that the model we consider here similarly generalizes, with a minimal membrane whose local surface tension depends on the entropy density at a given position and time.

III.2.2 Application to various initial states

To see the effect that the diffusion of charge has on entanglement growth, let us first consider the paradigmatic example of a maximally polarized domain wall initial state, wherein all sites on the left half of the chain are empty ($\mathfrak{q}(x\leq 0,t=0)=0$) and all sites on the right half are filled ($\mathfrak{q}(x>0,t=0)=1$). We imagine working in the thermodynamic limit, so that $x\in\mathbb{Z}$ and we can ignore boundary effects. At long times and large distances we expect the charge profile $\mathfrak{q}(x,t)$ to be well approximated by the solution of the continuum diffusion equation, which for these initial conditions reads $\mathfrak{q}(x,t)=\left(1+\text{erf}(x/\sqrt{Dt})\right)/2$, where $\text{erf}(x)$ is the error function 444This is most easily seen by taking the spatial derivative of the diffusion equation, which results in the same equation for $\partial_{x}\mathfrak{q}$. The initial condition for $\partial_{x}\mathfrak{q}$ is given by a delta function, $\partial_{x}\mathfrak{q}(x,0)=\delta(x)$, for the domain wall, which becomes a Gaussian of width $\sqrt{Dt}$ at time $t$. Integrating it up gives the result for $\mathfrak{q}(x,t)$ stated in the text.. The important property of this solution is that the width of the domain wall scales diffusively, as $\propto\sqrt{t}$, with time, and consequently the local entropy density $s(x,t)$ is also a function solely of $x/\sqrt{Dt}$. This indicates that significant entanglement growth can only occur within a region of order $\mathcal{O}(\sqrt{Dt})$ around the origin. At long times this diffusive growth is very slow compared to the linear growth of the entanglement surface that one would get for a homogeneous charge distribution. This means that the charge diffusion acts as a bottleneck for the entanglement growth: for the entanglement in the middle to grow further it has to ‘wait’ for the width of domain wall to increase. As a consequence, the height itself grows as $S(0,t)\propto\sqrt{t}$, as sketched in Fig. 1(a).

By performing the stochastic surface growth numerically, we find a scaling collapse of the form $S(x,t)=\sqrt{Dt}\,f(x/\sqrt{Dt})$ for some function $f$, as shown in Fig. 4(a). At long enough times, $Dt\gg 1$ (in units of the lattice spacing), the surface profile is always close to being saturated to the maximum value allowed by the charge distribution, with a distance between the two becoming constant. This also means that there is no allowed ‘space’ for the surface to develop large random fluctuations, which would be present for a flat surface - such KPZ fluctuations are suppressed by the slow dynamics of the charge that is coupled to the surface growth. In the opposite limit of $Dt\ll 1$ (which we cannot test in the random circuit model), the growth of the domain wall appears fast compared to the speed of the surface growth: the entanglement profile near the origin is not yet sensitive to finite width of the domain wall and we expect it to grow the same way as it would for a homogenous state at half filling.

While the domain wall constitutes an extreme case, where the charge densities at the two infinities take the two extreme values, [math] and $1$, a similar behavior occurs more generally if we take an initial state with a jump in the average charge at the origin. The simplest generalization is to take an initial state where $\mathfrak{q}(x\leq 0)=n$ and $\mathfrak{q}(x>0)=1-n$ for some $0<n<1/2$. In this case the dynamics of entanglement far from the origin, at positions $x\gg\sqrt{Dt}$ is the same as for a homogenous state with charge density $n$, i.e., it grows linearly with the rate set by $n$ and exhibits KPZ fluctuations. However, in a diffusively growing region around the origin, where the charge density has started equalizing, the entanglement grows faster and exhibits a ‘bump’ on top of the linearly growing background. Both the spatial size and the height of this bump is once again increasing as $\sqrt{t}$ (see Fig. 4(c)) and the KPZ fluctuations are suppressed in this region. This is shown in Fig. 4(b).

Another, slightly different scenario where the diffusion plays an important role, is an interface between two regions of different entropy densities. An extreme example of this is when a region at half filling, such as a Néel state, expands into empty space. In this case, in the middle of the half-filled region, far from the interface, we once again get linear growth with the usual KPZ fluctuations. Far on the other side of the interface the charge density vanishes and so does the entanglement. Near the interface there is a diffusively growing region that interpolates between these two extremes, where the entanglement obeys the scaling $S(x,t)=\sqrt{Dt}\,\,g(x/\sqrt{Dt})$, for some scaling function $g$. The middle region where this scaling is valid then penetrates linearly into the half-filled region. This is shown in Fig. 4(d). A similar situation would occur if the average charge densities on the two halves of the chain are initially $0<n_{L}<n_{R}<1-n_{L}$. In this case there is linear growth on both sides, but with different growth rates, so that there are two plateaus of different heights with an interpolating region in-between.

The initial states described above all have charge-imbalances between the left and right halves of the chain, such that in the thermodynamic limit it would take an infinite amount of time for them to become homogeneous. We can also consider initial states which have charge-inhomogeneities on some parametrically smaller length scale $1\ll\lambda\ll L$. An example of such a state is a step-like charge density wave (CDW) consisting of domains of $\lambda/2$ filled sites followed by $\lambda/2$ empty ones in a regular pattern. At short times, $t\ll\lambda^{2}/4D$, each domain wall in the initial state evolves independently, in the way described above, while far from the domain walls there is no dynamics. This results in an entanglement profile with peaks separated by a distance $\lambda/2$ from one another, whose width and height grows as $\sqrt{t}$. Once the peaks start to overlap, at times $t_{Q}\sim\mathcal{O}(\lambda^{2}/4D)$, i.e., when the domain walls start to melt together, the bottleneck disappears and the entanglement growth speeds up from $\sqrt{t}$ to linear in $t$ (see Fig. 5(b)). At the same time the entanglement profile starts to smooth out: the difference between the top and the bottom of a peak, $\Delta S\equiv S(0,t)-S(\lambda/4,t)$, decreases to zero, as shown by Fig. 5. We find numerically (see Fig. 6) that the average size of the peaks in the entanglement profile decreases exponentially, as $\Delta S\propto e^{-\frac{4\pi^{2}}{\lambda^{2}}t}$. We associate this with the first (diffusive) term on the right hand side of the first equation in Eq. (III.2.1).

Using the results above we can associate two distinct time scales to the evolution of the entanglement starting from a CDW, made out of regions of different charge densities of typical size $\lambda$. First there is the aforementioned time scale $t_{Q}\sim\lambda^{2}/D$, until which each domain wall evolves independently. By this time the entanglement profile develops peaks of size $\Delta S\propto\sqrt{Dt_{Q}}\propto\lambda$. At times $t\approx t_{Q}$ a second time regime begins, wherein the entanglement features smooth out exponentially as $\Delta S\propto\lambda e^{-8\pi^{2}\frac{t-t_{Q}}{t_{Q}}}$. The time needed for the entanglement profile to become sufficiently smooth, i.e., $\Delta S\ll 1$ is therefore $\frac{t-t_{Q}}{t_{Q}}\gg\log{\lambda}$, or equivalently $t\gg\frac{\lambda^{2}}{D}(1+\log{\lambda})$. Note that this is parametrically larger than the time scale it takes for the charge-distribution to become approximately flat, which merely requires $\frac{t-t_{Q}}{t_{Q}}\gg 1$.

To summarize our findings in this section, let us consider a rather broad set of initial states, that consists of different domains of typical size $\lambda$, each domain roughly homogeneous with some specific average charge density. For such a state, we can distinguish four different time scales. I) There is some short initial time scale $t_{\text{loc}}$ associated to reaching local equilibrium. II) This is followed by the transport of charge from domains of high density to neighboring regions of lower density. In the middle of the domains the entanglement grows linearly at the rate set by the local entropy density of the region. At the interfaces of different domains the transport of charge is accompanied by the appearance of ‘bumps’ in the entanglement profile. These bumps grow as $\sqrt{t}$ with respect to the entanglement in the middle of the domains which they separate (the one with the smaller entropy density if those are different). III) This bumpy entanglement profile develops up to times $t_{Q}\sim\mathcal{O}(\lambda^{2}/D)$, at which point the domains start to melt completely into one another. The entanglement growth at the domain walls speeds up to linear in time and the entanglement profile smooths out. There is another timescale associated to this smoothing which scales as $t_{Q}(1+\log{\lambda})$. IV) Eventually there is a fourth time regime, when the entanglement stops growing and saturates to its equilibrium value (assuming the system is finite). Note that the III) regime is only present if there is a separation of scales $\lambda\ll L$. If the initial inhomogeneities are on the scale of the entire system then equilibration of the charge density and the saturation of the entanglement profile happens simultaneously.

IV Spin-$\frac{1}{2}$ chains

We now compare the predictions of the simple effective model of the previous section with results on spin-$\frac{1}{2}$ chains. We first take the a version of the $N=1$ random circuit model defined in Sec. II, one where the arrangement of random gates follows a regular brick-wall pattern (Fig. 2(c)). In this case we can compute efficiently the annealed average of the $2^{\text{nd}}$ Rényi entropy (defined below). While Rényi entropies are expected to have very different dynamics from the von Neumann entropy in the case of a global quench Rakovszky et al. (2019), we find that the effects associated to charge inhomogeneities are similar to those predicted by our surface growth model for the latter. In particular, for a domain wall initial state we find a scaling collapse of the entanglement profile of the form $S(x,t)=\sqrt{t}f(x/\sqrt{t})$. In Sec. IV.2 we compare to exact results on a non-random, periodically driven spin chain.

IV.1 Spin-$\frac{1}{2}$ circuit model

The simplified model of the previous section relied on ignoring local fluctuations of the charge, and using only the average local density as the only relevant variable. Here we show that the features associated with inhomogeneities are quite similar even if we consider the spin-$\frac{1}{2}$ random circuit where fluctuations in local charge density are expected to be even larger. We take now a circuit where the unitary gates are applied in the regular ‘brick wall’ pattern shown in Fig. 2(c), and each site contains a single conserved spin ($N=1$).

We consider the annealed average of the second Rényi entropy, $S_{2}^{(a)}\equiv-\ln\left(\overline{\text{tr}(\rho_{A}^{2})}\right)$. The calculation of this quantity can be mapped the problem of evaluating a classical partition function Rakovszky et al. (2019). For $N=1$ this classical partition function can be evaluated by representing it as a two-dimensional tensor network, as we review in App. A. This allows us to treat the dynamics of $S_{2}^{(a)}$ for much larger systems and longer times then those available in the original quantum problem. The annealed average Rényi entropy provides a lower bound on the von Neumann entropy, $S_{2}^{(a)}(t)\leq\overline{S_{2}(t)}\leq\overline{S_{\text{vN}}(t)}$. While in Ref. Rakovszky et al., 2019 we argued that this bound is not tight, and $S_{2}^{(a)}$ behaves qualitatively differently from $S_{1}$ for homogenous initial states, it is nevertheless interesting to consider the behavior in the former for inhomogenous states. As we now show, several of the features predicted by our surface growth model in Sec. III appear also in the dynamics of $S_{2}^{(a)}$.

We once again start by considering a domain wall initial state. As argued before in Sec. III.2, since the charge density obeys the diffusion equation, it has a profile $\langle\hat{Q}_{x}(t)\rangle$ that depends only on the combination $x/\sqrt{Dt}$, where for the present case the diffusion constant is $D=1/2$. Assuming local equilibration, i.e., that the on-site reduced density matrix takes the form $\rho_{x}(t)\propto e^{-\mu(x,t)\hat{Q}_{x}}$ for some local chemical potential $\mu(x,t)$ (which, as we show in App. B.1, indeed holds after some short-time relaxation process), one can use subadditivity to upper bound for the von Neumann entropy as

[TABLE]

where $s(x^{\prime},t)$ is the one-site von Neumann entropy of site $x^{\prime}$ at time $t$. Therefore the von Neumann entropy cannot grow faster than $\sqrt{t}$. Clearly the same upper bound applies to all higher index Rényi entropies, as well as to $S_{2}^{(a)}$.

The annealed average, $S_{2}^{(a)}$, on the other hand, provides a lower bound for both $\overline{S_{1}}$ and $\overline{S_{2}}$. Therefore if $S_{2}^{(a)}$ also grows as $\propto\sqrt{t}$ then it follows that all of these quantities have to have the same diffusive growth. This is exactly what we find by evaluating the annealed average numerically: after some initial short-time dynamics its behavior at $x=0$ is well described by $S_{2}^{(a)}(0,t)=a\sqrt{t}+b$ for some constants $a,b$ as shown in Fig. 7. In fact we find that the profile of $S_{2}^{(a)}(x,t)$ is well approximated by the sum of the second Rényi entropies of the local one-site density matrices, $\sum_{x^{\prime}\leq x}s_{2}(x,t)$, where $s_{2}(x,t)=-\log(1-2\mathfrak{q}+2\mathfrak{q}^{2})$ is the Rényi entropy density associated to a state in local equilibrium, parametrized by the charge density $\mathfrak{q}=\mathfrak{q}(x,t)\equiv\overline{\langle\hat{Q}_{x}(t)\rangle}$. That is, despite the fact that $S_{2}$ does not obey subadditivity, its annealed average behaves as if it satisfied a modified version of the inequality in Eq. (8), where we replace the von Neumann entropy density with that of the second Rényi entropy. This interpretation is made plausible by a random circuit calculation in App. A, which shows that when charge fluctuations are ignored, $S_{2}^{(a)}$ obeys an equation similar to Eq. (6), but with a modified constant term.

IV.2 Deterministic Floquet spin chain

While random circuits have many advantages, we also wish to compare our predictions to deterministic models of unitary dynamics. While the KPZ fluctuations mentioned in Sec. III are particular to noisy dynamics, we expect that the imprints of charge transport, such as the $\sqrt{t}$ growth of entanglement for a domain wall, do generalize to the deterministic setting. We expect this to be the case since many of the arguments discussed above largely follows from the diffusive dynamics of conserved quantities, which is expected to hold generically at high temperatures even in deterministic systems Bloembergen (1949); Gennes (1958); Kadanoff and Martin (1963); Lux et al. (2014); Bohrdt et al. (2017).

To this end we also perform calculations in a Floquet (periodically driven) spin chain, introduced by the authors in Ref. Khemani et al., 2018, where it was shown to exhibit clear diffusive tails. Considering such a driven system has the advantage of simplicity; U(1) spin is the only conserved quantity. We leave it to future studies to consider other models, in particular those with energy conservation, and the appropriate generalizations of CDWs and domain walls in that context.

The model consists of a spin-$\frac{1}{2}$ chain evolved by a time-periodic Hamiltonian. A single driving sequence consists of four parts, with the so-called Floquet unitary given by

[TABLE]

where $\hat{X}_{r},\hat{Y}_{r},\hat{Z}_{r}$ are Pauli spin operators on site $r$. Every part of the drive individually conserves the spin z component, $[H_{a},\sum_{r}\hat{Z}_{r}]=0$ for $a=1,2,3,4$. We take the period time to be $T\equiv 4\tau=1$ and the couplings to be all order 1, namely $J_{z}^{(1)}=(\sqrt{3}+5)/6$, $J_{z}^{(2)}=\sqrt{5}/2$ and $J_{xy}=(2\sqrt{3}+3)/7$.

As emphasized in the previous section, the scaling observed for domain wall-like initial states, where both the size and the width of the features in the entanglement profile grow as $\propto\sqrt{t}$ is a direct consequence of diffusive transport - the upper bound derived above only relies on the assumptions of a) local equilibration and b) diffusive charge-transport, both of which are expected to apply to generic many-body systems at high temperatures. While this is only an upper bound, and in principle slower than $\sqrt{t}$ growth is possible, we do not expect this to be the case unless there are some further constraints on the dynamics 555Such is the case of integrable systems, which can exhibit logarithmic growth of the entanglement Alba and Heidrich-Meisner (2014); Vidmar et al. (2017), even though their transport properties are typically faster than diffusive.. Indeed, evolving a domain wall initial state with the Floquet unitary (IV.2) we find that the entanglement entropy across the middle of the chain grows slower than ballistically, approximately as $t^{0.5-0.6}$ 666Note that even in the random circuit results at $N=1$, shown in Fig. 7, the $~{}\sqrt{t}$ growth is only true asymptotically and at short times a slightly faster growth is observed.. This result, shown in Fig. 8, reinforces our expectation that the asymptotic growth should be $\propto\sqrt{t}$.

Considering the full entanglement profile, however, we do not observe the scaling collapse seen previously for the random circuit (i.e. a $\propto\sqrt{t}$ width). Our interpretation for this is the following. While the deterministic model we consider is expected to show diffusive transport at ‘high temperatures’, or in this case at filling fractions close to $1/2$, this is not the case when the filling is very small (or very close to $1$), which is the case far from the middle of the domain wall. In these regions there are very few particles (or very few holes), which therefore propagate without many scattering events. This suggest that in the deterministic model there should be two distinct regimes: near the origin where the domain wall has sufficiently melted and the average density is close to $1/2$ the diffusive scaling should apply, while at the tails of the domain wall, where the filling is close to (but different from) [math] or $1$ there is a ballistic region. We show some further data in support of this interpretation in App B.2.

V Logarithmic growth of number entanglement

Before concluding, let us briefly comment on the dynamics of the von Neumann entropy in charge-conserving systems in the absence of large-scale inhomogeneities. The discussion that follows applies to initial states that have a fixed amount of total charge in the entire system. A simple example of such a state is a Néel state, where exactly every second site is occupied. While this state is not completely homogenous, its charge-density becomes uniformly $1/2$ after some quick local equilibration process. The resulting state, however, is still far from global equilibrium, as indicated, among other things, by the fluctuations in the conserved charge. Consider a subsystem $A$, containing $L_{A}$ sites, and the total charge contained in it, as measured by the operator $\hat{Q}_{A}\equiv\sum_{x\in A}\hat{Q}_{x}$. In an equilibrium state, $\rho_{A}\propto e^{-\mu\hat{Q}_{A}}$ the charge variance should obey a volume law, $\text{Var}(Q_{A})\equiv\langle\hat{Q}_{A}^{2}\rangle-\langle\hat{Q}_{A}\rangle^{2}\propto L_{A}$. As suggested by our arguments made previously in Ref. Rakovszky et al., 2019 for the random circuit model, as well as those established more generally in the literature Lux et al. (2014), to leading order the evolution of the variance should be governed by diffusion, indicating that it takes a time $\propto L_{A}^{2}$ for it to reach its final value. Our goal is to understand how this diffusively slow relaxation affect the dynamics of entanglement.

To see the effects of charge diffusion more clearly it is worthwhile to write the von Neumann entropy as a sum of two contributions, one of which is related directly to the distribution of total charge in the subsystem. In order to achieve this, note that the conservation of a fixed total charge in the system implies that the reduced density matrix has a block-diagonal structure

[TABLE]

where $\rho_{A}^{(Q_{A})}$ corresponds to the block with total charge $Q_{A}$ in subsystem $A$, normalized to have unit trace. With this normalization, the prefactor $p_{Q_{A}}$ is exactly the probability of having a total number of $Q_{A}$ charges in $A$, given by the expectation value $\text{tr}(\rho_{A}\hat{P}_{Q_{A}})$ where $\hat{P}_{Q_{A}}$ is the projector appearing in the spectral decomposition $\hat{Q}_{A}=\sum_{Q_{A}}Q_{A}\hat{P}_{Q_{A}}$.

With this block-decomposition in hand one can rewrite the von Neumann entropy as Lukin et al. (2018)

[TABLE]

The first term is the entropy associated to the probability distribution $p_{Q_{A}}$, which we will refer to as the number entropy. The second term is the average entropy of the blocks which, following Ref. Lukin et al., 2018, we call the configurational entropy, as it is associated to superpositions between different ways of arranging the same number of total charges within the subsystem. Importantly, the number entropy is independently measurable in cold atom experiments Lukin et al. (2018).

For the types of initial states considered here, with roughly homogeneous charge-distributions, the mean value of the distribution $p_{Q_{A}}$ is close to its equilibrium value even after short times. However, this distribution is initially very narrow, with a variance $\text{Var}(Q_{A})$ that should grow as $\propto\sqrt{t}$, due to diffusion, as argued in Ref. Rakovszky et al., 2019. This provides an upper bound on the number entropy, as $S_{Q}\leq\ln\left(\sqrt{2\pi e\text{Var}(Q_{A})}\right)\approx\log{t^{1/4}}+\text{constant}$. We indeed confirm that the growth of $S_{Q}$ has this form, by performing the exact random circuit time evolution for a small system, as shown in Fig. 9. Note that approximating the entropy as the logarithm of the standard deviation is also consistent with the long-time saturation value, which should be $S_{Q}(t\to\infty)=\log{\sqrt{L_{A}}}+\text{const.}$ for the binomial distribution. As emphasized earlier, $S_{Q}$ itself is a measurable quantity and therefore this prediction could in principle be confirmed in experiments. We also expect this to manifest as a logarithmically growing component of the full von Neumann entropy (V), although in principle it is possible that this contribution is cancelled by some part of the configurational entropy.

VI Discussion

We investigated the dynamics of entanglement in situations where a one-dimensional system approaches equilibrium from an initial state with large scale inhomogeneities in some conserved quantity, such as particle number. We argued that all other things being equal, entanglement grows more slowly at very low and very high fillings: in particular, entanglement does not grow in empty regions, leading to a bottleneck which strongly limits the overall global growth of entropy. We made this intuition more precise by noting (using subadditivity) that the entanglement growth at each time step can be bounded by the local entropy density, which in turn is related to the local filling/chemical potential, once the system has locally equilibrated. Using this subadditivity bound as a heuristic equation of motion for the von Neumann entropy under a random circuit evolution Eq. (III.2.1), we constructed a stochastic surface growth model for the von Neumann entropy, resembling a noisy Kardar-Parisi-Zhang equation coupled to a diffusively evolving charge distribution. This can be recast as a modified minimal cut variational principle similar to that of Ref. Jonay et al., 2018, which suggests generalizations both to deterministic systems, and to higher dimensions.

We investigated the behavior of the surface growth equation of motion for various initial states. We showed that domain wall-like jumps in the initial filling lead to diffusively growing features in the entanglement profile. When such domain walls are combined into a Néel-like charge density wave of wavelength $\lambda$, this leads to an inhomogenous entanglement profile, with peaks at the positions of the initial domain walls. These peaks grow (and spread in space) as $\sim\sqrt{t}$ until the CDW melts at times $\sim\lambda^{2}$. The entropy profile eventually loses memory of the initial charge distribution, but only on a parametrically longer time scale $\sim\lambda^{2}(1+\log(\lambda))$.

Our surface growth model predictions were supported by a circuit averaged calculation of the annealed second Rényi entropy (which bounds below the von Neumann entropy). Furthermore we have demonstrated that the deterministic, non-noisy model of Eq. (IV.2) also exhibits the diffusive growth of entanglement for the case of a single domain wall, (albeit with some additional subtleties related to regions with very low density of excitations where transport becomes ballistic). Finally, motivated by recent experiments Lukin et al. (2018), we investigated the dynamics of the so-called number entanglement, and argued that it grows logarithmically in time in systems with diffusive transport. Our argument relied on the assumption that the charge variance of extensive regions tends to grow as $\propto\sqrt{t}$, an assumption based on general hydrodynamic considerations and verified by a random circuit calculation.

Our work suggests further avenues for research. It would be interesting to see whether our heuristically motivated equation of motion for the entanglement can be put on firmer ground in the random circuit setting, perhaps using the methods of Ref. Zhou and Nahum, 2018. Furthermore, while a generalization of our model to higher dimensions, through the minimal cut/membrane picture, appears natural, establishing it more rigorously, and investigating its properties remains an open problem. Another interesting avenue is understanding entanglement growth in situations where there is an interplay between particle/spin transport and energy conservation. For example, the dynamics of charge domain walls in U(1)-symmetric Hamiltonian systems, can be more complex, due to such interplay constraining the domain wall motion. Investigating entanglement dynamics in this case, for generic non-integrable systems is an interesting open question.

Lastly, it would be interesting to pin down the connection between the minimal cut formalism and holographic calculations Mezei (2018); do the diffusive characteristics of entropy growth show up in analogous holographic quench calculations?

Acknowledgements.

We thank Herbert Spohn, Adam Nahum and Juah Garrahan for useful discussions. FP acknowledges the support of the DFG Research Unit FOR 1807 through grants no. PO 1370/2- 1, TRR80, the Nanosystems Initiative Munich (NIM) by the German Excellence Initiative, the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2111-390814868 and the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement no. 771537). CvK is supported by a Birmingham Fellowship.

Appendix A Mapping average purity to classical partition function

As noted in the main text, the simplest entanglement measure which we can calculate tor our random circuit models is the annealed average of the second Rényi entropy, which is given by the logarithm of the average purity $\mathcal{P}$,

[TABLE]

Since this is a quantity that only involves four instances of the time evolution operator $U$, the averaging over a single gate of the circuit can be done in a relatively straightforward manner, resulting in a 2D tensor network representation of the average purity. The boundary conditions of this partition function at times [math] and $t$ depend on the initial state $\rho(0)$ and on the choice of the subsystem $A$, respectively. This representation allows us to efficiently evaluate $S_{2}^{(a)}(t)$, by contracting the 2D tensor network via the time-evolving block decimation algorithm Vidal (2003), resulting in the data shown in Sec. IV.1 and below in App. B.1. The same mapping has been developed in Ref. Rakovszky et al., 2018 (generalizing the calculation of Ref. Nahum et al., 2017 to the U(1)-symmetric case) to calculate other quantities, and we summarize it briefly below.

The basic building block of the 2D partition function is a 4-leg tensor, which we obtain by averaging over four copies of a single unitary gate $U$. The operator $U$ acts on a Hilbert space $\mathcal{H}^{2}$ (i.e. the Hilbert space of two neighboring sites of the chain), which can be decomposed into sectors with fixed total charges as $\mathcal{H}^{2}=\bigoplus_{Q}\mathcal{H}^{2}_{Q}$, each of which has a dimension $d_{Q}$. The effective tensor associated to this gate then has the following form Rakovszky et al. (2018):

[TABLE]

where the states $|{{\mathcal{I},\mathcal{J}}}\rangle$ on the four-copy Hilbert space are defined as

[TABLE]

These states belong to four copies of the total (two-site) Hilbert space. In order to contract the tensor network one then needs to split them up to states living on four copies of a single-site Hilbert space, which in principle can be done in a many different ways, depending on the choice of local basis.

Each unitary gate therefore maps to such a 4-leg tensor, which then need to be contracted, as defined by the geometry of the circuit (e.g. the ‘brick wall’ geometry of Fig. 2(c)). The final step in calculating the average purity is to specify the boundary conditions. We note that the purity of subsystem $A$ can be rewritten as Hastings et al. (2010); Abanin and Demler (2012); Hayden et al. (2016); Nahum et al. (2018a); Rakovszky et al. (2019)

[TABLE]

where $\rho\otimes\rho$ is two copies of the density matrix of the whole system, and $\mathcal{S}_{A}$ is an operator on this doubled Hilbert space that swaps the two copies inside subsystem $A$ and acts as the identity on its complement. The trace on the right hand side should be understood as a sum over a complete set of states in the two-copy Hilbert space. The boundary conditions of the 2D tensor network representation therefore correspond to $\rho(0)\otimes\rho(0)$ on the boundary at time [math] and by $\mathcal{S}_{A}$ on the boundary at time $t$. In particular, denoting by $\mathcal{H}^{1}$ the single site Hilbert space, the latter reads

[TABLE]

where

[TABLE]

The other boundary condition depends on the initial state $\rho(0)=|{{\psi(0)}}\rangle\langle{{\psi(0)}}|$. For product states in the local charge basis, relevant in the present context, the two on-site states that appear are $|{{0000}}\rangle$ and $|{{1111}}\rangle$.

A.1 Equation of motion for well defined local charge

Apart from providing an efficient way of numerically evaluating $S_{2}^{(a})(t)$ in the random circuit, utilized in Sec. IV.1, the above calculation can also be used to derive an update rule in the case when a 2-site unitary is acting on a bond with fixed amount of charge, similarly to what has been done for the Hartley entropy in Sec. III.1.2. Here we show that this leads to an update rule similar to Eq. (6), but with a modified constant term, reflecting the differences between the von Neumann and second Rényi entropies.

For the purposes of understanding the effect of a single gate, it is enough to look at the two-site version of the swap operator, i.e. $\mathcal{I}\otimes\mathcal{K}$. Applying the operator (13) to this state we get two terms, which read

[TABLE]

where we have defined

[TABLE]

in terms of the dimension of the charge $a$ sector, $\eta_{a}\equiv\dim\mathcal{H}^{1}_{Q}$, and the total single site dimension $\eta\equiv\sum_{a}\eta_{a}$.

To arrive at the equation of motion we are going to assume that the local charge on the two sites in question is exactly $\bar{Q}$, so that

[TABLE]

where $P_{Q}$ is the operator that projects onto the charge $Q$ sector of the two-site Hilbert space. Therefore, when evaluated on a density matrix that satisfied the above relation, the sum of the two terms in Eq. (A.1) can be replaced with

[TABLE]

The right hand side is therefore proportional to the sum of the ‘purity superoperators’ at the two neighboring bonds. The last remaining step is to evaluate the prefactor. In the large-$N$ limit it simplifies to $\frac{\eta f_{\bar{Q}\bar{Q}}}{d^{2}_{\bar{Q}}}$. A short calculation reveals that this quantity is just the one-site purity associated with the maximally entangled two site density matrix with charge exactly $\bar{Q}$, namely $\rho=P_{\bar{Q}}/d_{\bar{Q}}$.

Using the assumption of well defined local charge, along with the large-$N$ limit, we therefore arrive at a closed equation of motion for the Purity $\mathcal{P}(x,t)$ under the effect of a local 2-site gate at position $x$, which reads

[TABLE]

Here $s_{2}(\bar{Q})\equiv\ln\left(\frac{\eta f_{\bar{Q}\bar{Q}}}{d^{2}_{\bar{Q}}}\right)$ is the density of second Rényi entropy, associated to the local charge density on the bond $x,x+1$. Taking the logarithm of Eq. (21) (with base $2^{N}$) and once again making use of the large-$N$ limit results in a very simple update rule for the annealed average second Rényi entropy under the effect of a 2-site gate:

[TABLE]

This has the same form as the update rule (6) we used in our effective surface growth model for the von Neumann entropy, except for the fact that the local density of von Neumann entropy has been replaced by the local density of $S_{2}$. This update rule is consistent with our numerical results in Sec. IV.1, where we observed that the profile of $S_{2}^{(a)}$ starting from a domain wall at long times behaves as if it obeyed sub-additivity. This is despite the fact that the update rule (22) clearly misses the important effect of charge fluctuations for homogenous states, uncovered in Ref. Rakovszky et al., 2019.

Appendix B Further numerical results on spin-$\frac{1}{2}$ chains

In this appendix we gather some further data on inhomogenous quenches, both in the random circuit and in the Floquet model, to complement the results shown in the main text.

B.1 Local equilibration in the spin-$\frac{1}{2}$ random circuit

Our general argument, putting an upper bound proportional to $\sqrt{t}$ on the growth of von Neumann entropy for a domain wall, relied on the notion of local equilibration. Here we check numerically that this assumption indeed holds in the spin-$\frac{1}{2}$ local random circuit we considered in the main text. This consists of two parts. First of all, in order for local equilibration to $\overline{\langle\hat{Q}_{x}(t)\rangle}$ to make sense, the circuit-to-circuit fluctuations in $\langle\hat{Q}_{x}(t)\rangle$ around this average should become small. This means that we can assign a local charge to each site, which should be approximately the same for all realizations of the circuit, making the right hand side of Eq. (8) well-defined. Second, we need the one-site density matrices to be close to the Gibbs state defined by this local charge density. As we show in Fig. 10 both of these requirements are satisfied at sufficiently long times.

To check the first requirement, we compute numerically the statistical variance of the local charge in the random circuit, i.e. the variance between circuit realizations of the quantum expectation value $\langle\hat{Q}_{x}\rangle$. This can be computed as a 2D partition function, similarly to the calculation of $S_{2}^{(a)}$ outlined in App. A. We find that the variance decays in time, as shown by Fig. 10(b), indicating that at long times it is meaningful to consider locally equilibrated states depending on a circuit-independent average local charge.

To check if the on-site density matrix is indeed locally equilibrated, we compare its annealed average second Rényi entropy with the entropy density expected for an eqilibrium state with the same charge density. Note that since the on-site density matrix only has two parameters, one of which is fixed by the charge density, the equivalence of Rényi entropies is sufficient to establish local equilibrium. We indeed find that the local Rényi entropy tends to its equilibrium value at long times; this is shown in Fig. 10(c). In both cases we show results for different bond dimensions $\chi$ used in the evaluation of these quantities as 2D tensor networks to indicate where the results are converged.

B.2 Entanglement profile in the Floquet model

As stated in Sec. IV.2 of the main text, while we do observe the expected sub-ballistic growth for the half-chain entanglement of a domain wall initial state in the Floquet spin chain model, we find that the spatial width of the entanglement profile spreads faster than $\sqrt{t}$. We interpreted this by noting that at in the tails of the melting domain wall, far on the left (right) from its center, the conserved spin density is close to minimal (maximal). In these regions, a weakly-interacting quasi-particle picture should therefore be approptiate, implying ballistic spreading. This is supported by considering the spatial profile of the charge, which shows these two distinct regimes: tails that spread linearly and a middle region which spreads out more slowly, as shown in the upper panels of Fig. 11

Based on the above interpretation it is expected that the diffusive (spatial) spreading of the entanglement profile should be more pronounced for a generalized domain wall, like the ones considered above for random circuits, where the average charge density on the two sides is $0<n_{L}<n_{R}<1$. Since such states have a finite density of particles everywhere, diffusion should hold (at long times) even far from the origin. Although for such states we are more limited in the times we can approach numerically, we indeed find a better scaling collapse with $\sqrt{t}$ both for the spin and for the entanglement profile. This is shown in the lower panels of Fig. 11.