Graph-based methods are popular in dependency parsing for decades. Recently, <cit.> propose a headed span-based method. Both of them score all possible trees and globally find the highest-scoring tree. In this paper, we combine these two kinds of methods, designing several dynamic programming algorithms for joint inference. Experiments show the effectiveness of our proposed methods[Our code is publicly available at <https://github.com/sustcsonglin/span-based-dependency-parsing>.].

§ INTRODUCTION

Dependency parsing is an important task in natural language processing. It has numerous applications in downstream NLP tasks.

There are many methods to tackle projective dependency parsing.
In this paper, we focus on two kinds of global methods: graph-based methods and headed span-based methods. They both score all possible trees and globally find the highest scoring tree. The difference is how they score dependency trees. In the simplest first-order graph-based methods <cit.>, they decompose the score of a dependency tree into the score of dependency arcs. In second-order graph-based methods <cit.>, they additionally score the sibling factors. There are more other higher-order graph-based methods <cit.>.  In contrast, in headed span-based methods <cit.>, they decompose the score into the score of headed spans. In a projective tree, the subtree rooted at each word covers a span in the surface order, the headword-span pair is referred to as headed span. Figure <ref> shows an example projective parse tree and all its headed spans.
Graph-based methods take word pairs as basic units for modeling while headed span-based methods directly model subtrees by scoring headed spans, so they have a very different view on the dependency trees, and combining them is intuitively helpful.

dg edges
[the [the, name=x1] ]
[child, name=x2] ]
[reads, name=x3]
[a [a, name=x4] ]
[book, name=x5] ] ]
ıin 2, ..., 5
(boundaryȷ) [below =1ex of x1] at  ( $ (x\j.east)!0.5!(x\i.west) $ ) ;

(boundary0) [left=4ex of boundary1] [];
(boundary5) [right=4ex of boundary4] [];
(boundary0) –  (boundary1) node [midway, below, blue] the;
(boundary3) – (boundary4) node [midway, below, blue] a;
($(boundary0)+(0,-4ex)$) – ($(boundary2)+(0,-4ex)$) node [midway, below, blue] child;
($(boundary3)+(0,-4ex)$) – ($(boundary5)+(0,-4ex)$) node [midway, below, blue] book;
($(boundary0)+(0,-8ex)$) – ($(boundary5)+(0,-8ex)$) node [midway, below, blue] reads;

An example projective dependency parse tree with all of its headed spans.

We design several dynamic programming algorithms for joint inference.
We first design an $O(n^4)$ algorithm adapted from the Eisner-Satta algorithm <cit.> to combine first-order graph-based and headed span-based methods using the hook trick. Then we design two $O(n^3)$ algorithms adapted from the Eisner algorithm <cit.> and its second-order extension <cit.> to combine first/second-order graph-based methods and headed span-based methods using the head-splitting trick: we decompose the score of headed spans into two terms assuming that the score of left span boundary is independent to that of the right span boundary for each headword. We conduct extensive experiments on PTB, CTB, and UD, showing the effectiveness of our proposed methods.

§ MODEL

§.§ Scoring

Given input sentences $x_1, ..., x_n$, we add <bos> and <eos> at the $x_0$ and $x_{n+1}$. We use the last layer of BERT <cit.> and use mean-pooling (i.e., averaging all subwords embeddings) to obtain the word-level embeddings $e_i$[For some dataset requiring the use of gold POS tags, we additionally concatenate the POS tag embedding to obtain $e_i$]. Then we feed $e_0, ..., e_{n+1}$ into a three-layer BiLSTM <cit.> network to get $c_0, ..., c_{n+1}$, where $c_i = [f_i; b_i]$, $f_i$ and $b_i$ are the forward and backward hidden state of the last BiLSTM layer at position $i$ respectively. We use $h_k=[f_k, b_{k+1}]$ to represent the $k$th boundary and use $e_{i,j}=h_j-h_i$ to represent span $(i,j)$. Then we feed $c_i, c_k$ into a deep biaffine function <cit.> to score $s^{\text{arc}}_{i,j}$ for arc $x_i \rightarrow x_j$; feed $c_i, c_k, c_j$ into a deep triaffine function <cit.> to score $s^{\text{sib}}_{i,j,k}$ where $x_i \rightarrow \{x_j,x_k\}$ with $k<j<i$ or $i<j<k$; feed $e_{i, j}, c_k$ to a deep biaffine function to score $s^{\text{span}}_{i,j,k}$ for headed span $(i, j, k)$ where $x_k$ is the headword of $(i,j)$; feed $c_i, h_j$ to a deep biaffine function to score $s^{\text{left}}_{i,j}$ where $j$ is the left boundary of the headed span of $x_i$; feed $c_i, h_k$ to a deep biaffine function to score $s^{\text{right}}_{i,j}$ where $k$ is the right boundary of the headed span of $x_i$.

§.§ Learning

We decompose the training loss $L$ into $L_{\text{parse}} + L_{\text{label}}$. For $L_{\text{parse}}$, we use the max-margin loss  <cit.>:

\begin{align*}
L_{\text{parse}} = \max (0, \max _{y^{\prime} \neq y}(s(y) + \Delta(y^{\prime}, y) - s(y))
\end{align*}

where $\Delta$ measures the similarity between the incorrect tree and gold tree $y$. Here we let $\Delta$ to be the Hamming loss on the basic units we model (arcs, headed spans, or half of the headed spans depending on the setting). We use the same label loss $L_{\text{label}}$ as <cit.>.

§ PARSING

We use the parsing-as-deduction framework <cit.> to describe the parsing algorithms of our proposed models:  (<ref>),  (<ref>),  (<ref>).

§.§ $O(n^4)$ modified Eisner-Satta algorithm

In this case, the score of a dependency tree $y$ is defined as:

\[
s(y) = \sum_{(x_i \rightarrow x_j) \in y} s^{\text{arc}}_{i,j} + \sum_{(l_i, r_i, x_i) \in y} s^{\text{span}}_{l_i, r_i, i}
\]

We design a dynamic programming algorithm adapted from the Eisner-Satta algorithm, which uses the hook trick to accelerate bilexicalized PCFG parsing.  The axiom items are $\tikz[baseline=-10pt]{\triangle[0.6][0.6]{i-1}{i}{i}{\phantom{0}}}$ with initial score 0 and the deduction rules are listed in Figure <ref>.  The original Eisner-Satta algorithm does not distinguish the “finished” span and “unfinished” span. An “unfinished” span can absorb a child span to form a larger span, while in a “finished” span, the headword has already generated all its children, so it cannot expand anymore and it is exactly the headed span for the given headword. By explicitly distinguishing the “unfinished“ span and “finished“ span, we can incorporate $s^{span}$ via the newly introduced rule .  We then modify the rule  and  accordingly as only the “finished” span can be attached.

§.§ $O(n^3)$ modifed Eisner algorithm

We decompose $s^{\text{span}}_{l_i, r_i, i}$ into two terms:

\[
s(y) = \sum_{(x_i \rightarrow x_j) \in y} s^{\text{arc}}_{i,j} + \sum_{(l_i, r_i, x_i) \in y} (s^{\text{left}}_{i, l_i} + s^{\text{right}}_{i, r_i})
\]

so that we can decrease the inference time complexity from $O(n^4)$ to $O(n^3)$ by modifying the Eisner algorithm.  The axiom items are  $\tikz[baseline=-10pt]{\righttrianglee[0.6][0.4]{i}{i+1}{}}$ and $\tikz[baseline=-10pt]{\lefttrianglee[0.6][0.4]{i}{i+1}{}}$
with initial score 0 and the deduction rules are shown at the first two rows of Figure <ref>.  Similar to the discussion in the previous section,  the origin Eisner algorithm does not distinguish the “finished” complete span and “unfinished” complete span. A “unfinished” complete span can absorb another complete span to form a larger incomplete span, while a “finished” complete span has no more child in the given direction and thus cannot expand anymore. We introduce new rules  and  to incorporate the left or right span boundary score respectively, and other rules are adjusted accordingly.

$\inferrule{\tikz[]{\triangle[0.4][0.4]{i}{h}{k}{s_1}} \quad
\tikz[]{\triangleleft[0.4][0.4][1.2]{k+1}{}{j}{s_2}}   } {
\tikz[]{\triangle[0.6][0.6]{i}{h}{j}{s_1 + s_2}}

\colorbox{red!20}{
\tikz[]{\triangleclose[0.4][0.4]{i}{c}{j}{s_1}} }} {
\tikz[]{\triangleleft[0.4][0.4][1.2]{i}{}{j}{s_1 + s^{\text{arc}}_{h,c}}

\tikz[]{\triangle[0.6][0.6]{i}{h}{j}{s_1}} } {
\colorbox{red!20}{
\tikz[]{\triangleclose[0.6][0.6]{i}{h}{j}{s_1+s_{h, i-1, j}^{\text{span}}

$\inferrule{\tikz[]{\triangleright[0.4][0.4][1.2]{i}{}{k-1}{s_1}} \quad
\tikz[]{\triangle[0.4][0.4]{k}{h}{j}{s_2}}   } {
\tikz[]{\triangle[0.6][0.6]{i}{h}{j}{s_1 + s_2}}

\colorbox{red!20}{
\tikz[]{\triangleclose[0.4][0.4]{i}{c}{j}{s_1}} }
} {
\tikz[]{\triangleright[0.4][0.4][1.2]{i}{}{j}{s_1 + s_{h,c}^{\text{arc}}}

Deduction rules for our modified Eisner-Satta algorithm <cit.>. Our modifications highlighted in red. All deduction items are annotated with their scores.

\tikz[baseline=-15pt]{\trapezoidd[0.6][0.6][0.3]{i}{k}{s_1}} \quad
\colorbox{red!20}{
\tikz[baseline=-10pt]{\rightclosetrianglee[0.5][0.4]{k}{j}{s_2}}
} {
\tikz[baseline=-10pt]{\righttrianglee[0.8][0.6]{i}{j}{s_1+s_2}}

\tikz[baseline=-9pt]{\righttrianglee[0.5][0.5]{i}{k}{s_1}}
\colorbox{red!20}{
\tikz[baseline=-9pt]{\leftclosetrianglee[0.5][0.5]{k+1}{j}{s_2}}
} {
\tikz[baseline=-10pt]{\trapezoidd[0.6][0.6][0.3]{i}{j}{s_1+s_2+s^{\text{arc}}_{i, j}}}

\tikz[baseline=-10pt]{\righttrianglee[0.8][0.6]{i}{j}{s_1}}
} {
\colorbox{red!20}{
\tikz[baseline=-10pt]{\rightclosetrianglee[0.8][0.6]{i}{j}{s_1 + s^{\text{right}}_{i, j}}}}

\colorbox{red!20}{
\tikz[baseline=-9pt]{\leftclosetrianglee[0.5][0.4]{j}{k}{s_1}}}\quad
\tikz[baseline=-9pt]{\trapezoiddleft[0.6][0.4][-0.3]{k}{i}{s_2}}
\tikz[baseline=-10pt]{\lefttrianglee[0.8][0.6]{j}{i}{s_1+s_2}}

\colorbox{red!20}{
\tikz[baseline=-9pt]{\rightclosetrianglee[0.5][0.5]{j}{k-1}{s_1}}}
\tikz[baseline=-9pt]{\lefttrianglee[0.5][0.5]{k}{i}{s_2}}
} {
\tikz[baseline=-6pt]{\trapezoiddleft[0.6][0.3][-0.3]{j}{i}{s_1+s_2+s^{\text{arc}}_{i,j}}}

\tikz[baseline=-10pt]{\lefttrianglee[0.8][0.6]{j}{i}{s_1}}
\colorbox{red!20}{
\tikz[baseline=-10pt]{\leftclosetrianglee[0.8][0.6]{j}{i}{s_1 + s^{\text{left}}_{i, j-1}}}

\colorbox{red!20}{
\tikz[baseline=-9pt]{\rightclosetrianglee[0.5][0.5]{j}{k-1}{s_1}}} \quad
\colorbox{red!20}{
\tikz[baseline=-9pt]{\leftclosetrianglee[0.5][0.5]{k}{i}{s_2}}
} {
\tikz[baseline=-6pt]{\square[1][0.5]{j}{i}{s_1+s_2}}

\tikz[baseline=-6pt]{\square[0.6][0.3]{j}{k}{s_1}} \quad
\tikz[baseline=-6pt]{\trapezoiddleft[0.6][0.3][-0.3]{k}{i}{s_2}}
} {
\tikz[baseline=-6pt]{\trapezoiddleftt[0.6][0.3][-0.3]{j}{i}{
s_1+s_2+s^{\text{arc}}_{i,j}}{+ s^{\text{sib}}_{i, k, j}}}}$

\tikz[baseline=-15pt]{\trapezoidd[0.6][0.6][0.3]{i}{k}{s_1}} \quad \tikz[baseline=-6pt]{\square[0.6][0.3]{k}{j}{s_1}}
} {
\tikz[baseline=-6pt]{\trapezoiddd[0.6][0.6][0.3]{i}{j}{
s_1+s_2+s^{\text{arc}}_{i,j}}{+ s^{\text{sib}}_{i, k, j}}}}

Deduction rules for our modified Eisner algorithm <cit.> (first two rows) and its second-order extension <cit.> (all rows). Our modifications highlighted in red. All deduction items are annotated with their scores.

§.§ $O(n^3)$ modified second-order Eisner algorithm

We further enhance the model with sibling information:

\begin{align*}
s(y) = &  \sum_{(x_i \rightarrow x_j) \in y} s^{\text{arc}}_{i,j} + \sum_{(x_i \rightarrow \{x_j, x_k\})\in y} s^{\text{sib}}_{i,j,k}\\
&	  + \sum_{(l_i, r_i, x_i) \in y} (s^{\text{left}}_{i, l_i} + s^{\text{right}}_{i, r_i})
\end{align*}

Similarly, we modify the second-order extension of the Eisner algorithm <cit.> by distinguishing “unfinished” and “finished” complete spans.  The additional deductive rules for second-order parsing are shown in the last row of Figure <ref> and the length of the “unfinished” complete span is forced to be 1 in the rule  and .

bg    	 ca     cs    de    en    es    fr     it   nl   no   ro   ru    Avg

MFVI2O     91.30     93.60     92.09     82.00     90.75     92.62     89.32     93.66     91.21     91.74     86.40     92.61      90.61

$\textit{Biaffine+MM}^{\dagger}$     90.30     94.49     92.65     85.98     91.13    93.78     91.77    94.72    91.04    94.21    87.24     94.53    91.82

Span     91.10     94.46    92.57    85.87     91.32     93.84     91.69     94.78     91.65     94.28     87.48    94.45     91.96

Biaffine+2O+MM      91.58     94.48      92.69     85.72    91.28      93.80     91.89     94.88      91.30     94.23      87.55     94.55     92.00

1O+Span     91.44    94.54    92.68    85.75    91.23    93.84    91.67     94.97     91.81    94.35     87.17    94.49    91.99

1O+Span+Headsplit    91.46    94.53    92.63    85.78    91.25    93.77     91.91    94.88    91.59     94.18    87.45    94.47    91.99

2O+Span+Headsplit     91.82     94.58    92.59    85.65    91.28     93.86    91.80    94.75    91.50     94.40     87.71    94.51     92.04

Labeled Attachment Score (LAS) on twelves languages in UD 2.2. We use ISO 639-1 codes to
represent languages.  $\dagger$ means reported by <cit.>. MFVI2O: <cit.>. Span: <cit.>.

2cPTB     2cCTB

UAS     LAS     UAS     LAS

2c+$\text{XLNet}_{\text{large}}$   2c+$\text{BERT}_{\text{base}}$

HPSG$^\flat$     97.20     95.72     -     -

HPSG+LAL$^\flat$     97.42     96.26     94.56     89.28

2c+$\text{BERT}_{\text{base}}$   2c+$\text{BERT}_{\text{base}}$

RNGTr     96.66     95.01     92.98     91.18

2c+$\text{BERT}_{\text{large}}$   2c+$\text{BERT}_{\text{base}}$

MFVI2O     96.91     95.34     92.55     91.69

HierPtr     97.01     95.48     92.65     91.47

$\textit{Biaffine+MM}^{\dagger}$     97.22     95.71     93.18     92.10

Span     97.24     95.73     93.33     92.30

Biaffine+2O+MM     97.28     95.73     93.42     92.34

1O+Span     97.26     95.68     93.56     92.49

1O+Span+HeadSplit     97.30      95.77     93.46    92.42

2O+Span+HeadSplit     97.23     95.69     93.57     92.47

Results for different model on PTB and CTB. $^\flat$: They use additional constituency tree data. $^{\dagger}$ means reported by <cit.>. HPSG: <cit.>; HPSG+LAL: <cit.>; RNGTr: <cit.>; HierPtr: <cit.>.

§ EXPERIMENTS

§.§ Setup

We conduct experiments on in Penn Treebank (PTB) 3.0 <cit.>, Chinese Treebank (CTB) 5.1 <cit.> and 12 languages on Universal Dependencies (UD) 2.2. We use the same data process, evaluation methods, and hyperparameters as <cit.> for fair comparison and we refer readers to their paper for details due to the limit of space. We set the hidden size of the Triaffine function as 300 additionally. The reported values are averaged over three runs with different random seeds.

§.§ Main result

Table <ref> and <ref> show the results on UD, PTB and CTB respectively. We additionally reimplement  by replacing the TreeCRF loss of <cit.> with max-margin loss for fair comparison.

We draw following observations:  (1) Second-order information is still helpful even with powerful encoders (i.e., BERT).  outperforms  in almost all cases. (2) Combining graph-based and span-based methods is effective in almost all cases. Both  and  beat  and  on PTB, CTB and UD.  outperforms  on CTB and UD.   (3) Decomposing the headed span scores is effective.  has a comparable performance to  while manages to decrease the complexity from $O(n^4)$ to $O(n^3)$. We speculate that powerful encoders mitigate the issue of independent scoring.

§ CONCLUSION

In this paper we have presented several dynamic programming algorithms to combine graph-based and headed span-based methods. Experiments show the effectiveness of our proposed methods.


# Supplemental Material for Noise Induced Universal Diffusive Transport in
Fermionic Chains

Christopher M. Langlett  Department of Physics & Astronomy, Texas A&M
University, College Station,Texas 77843, USA  Shenglong Xu  Department of
Physics & Astronomy, Texas A&M University, College Station,Texas 77843, USA

###### Abstract

This supplemental material’s purpose is to clarify the conclusions drawn in
the main text. We first give an overview of calculating the two-point
correlation function utilizing free fermion techniques. Next, we derive the
operator equation of motion that leads to the effective Liouvillian. In the
third section, we explicitly derive the differential equation for the operator
amplitudes. Solving the diffusion constant for different noisy and static
combinations constitutes the remaining sections. The final section outlines
the numerical methods for simulating the differential equations.

###### Contents

  1. S1 Fermion Operator Identities.
  2. S2 Free-Fermion Techniques.
    1. S2.1 Efficient Free-Fermion Simulations.
  3. S3 Operator Equation of Motion.
  4. S4 Details in Calculating the Coefficient Differential Equation.
    1. S4.1 Coherent Term from Static Hopping
    2. S4.2 Incoherent Term from Noisy Hopping.
    3. S4.3 Operator Equation of Motion.
    4. S4.4 Mapping to $\ell$ and $\mathcal{R}$.
    5. S4.5 Incoherent Term from Noisy Potential.
    6. S4.6 Mapping to $\ell$ and $\mathcal{R}$
  5. S5 Diffusion Constant for Static Hopping with Bond Noise.
  6. S6 Diffusion Constant for Static Hopping and Onsite Noise.
  7. S7 Diffusion Constant for Static Hopping with Bond and Onsite Noise.
  8. S8 Diffusion Constant for Linear Potential with Bond Noise.
    1. S8.0.1 Approximate solution to Bound State Energy.
  9. S9 Diffusion Constant for Linear Potential with Onsite Noise.
  10. S10 Diffusion Constant for General Operator Equation of Motion.
  11. S11 Numerical Methods for Determining the Coefficients.

##  S1 Fermion Operator Identities.

In this appendix, we provide helpful fermion operator identities used to
derive the equations of motion. We denote the spinless fermionic creation and
annihilation operators by $(c^{\dagger}_{\alpha},c_{\alpha})$, where $\alpha$
indicates the site index. These operators obey the algebra,

$\displaystyle\begin{split}\\{c_{\alpha},c_{\beta}\\}=\\{c^{\dagger}_{\alpha},c^{\dagger}_{\beta}\\}=0\\\
\\{c^{\dagger}_{\alpha},c_{\beta}\\}=\delta_{\alpha,\beta}.\end{split}$ (S1)

Further, defining the density operator
$n_{\alpha}=c^{\dagger}_{\alpha}c_{\alpha}$, we obtain,

$\displaystyle[n_{\alpha},c_{\beta}]$
$\displaystyle=-\delta_{\alpha,\beta}c_{\beta},\quad[n_{\alpha},c^{\dagger}_{\beta}]=\delta_{\alpha,\beta}c^{\dagger}_{\beta}$
$\displaystyle\\{n_{\alpha},c^{\dagger}_{\beta}\\}$
$\displaystyle=\delta_{\alpha,\beta}c_{\alpha}^{\dagger}+2c_{\beta}^{\dagger}n_{\alpha},\quad\\{n_{\alpha},c_{\beta}\\}=\delta_{\alpha,\beta}c_{\alpha}+2n_{\alpha}c_{\beta}$
(S2)

These relations were utilized in deriving the differential equations for the
coefficients in the operator expansion.

##  S2 Free-Fermion Techniques.

We now calculate the density-density correlation function using standard
methods. In the clean limit, $\Gamma_{x,y}=0$, with $J_{x,x}=0$, the exchange
couplings are uniform, and the Hamiltonian describes a one-dimensional free-
fermion. At this point, the Hamiltonian regains spatial and time translation
invariance, thereby in Fourier space, the fermionic operators become,

$\displaystyle\begin{split}c_{x}&=\frac{1}{\sqrt{N}}\sum_{k}e^{ik\cdot
x}c_{k}\\\ c^{\dagger}_{x}&=\frac{1}{\sqrt{N}}\sum_{k}e^{-ik\cdot
x}c^{\dagger}_{k}\\\ \end{split}$ (S3)

here $k$ is the Brillouin zone momentum, and the Hamiltonian,

$H=\sum_{k}\mathcal{E}_{k}c_{k}^{\dagger}c_{k}$ (S4)

with the band dispersion, $\mathcal{E}_{k}=2\cos(k)$ (lattice spacing $a=1$).
In the momentum basis, time-evolution of the operators
$c^{\dagger}_{k}~{}(c_{k})$ is simple, taking the form,

$\displaystyle\begin{split}c_{k}(t)&=e^{-i\mathcal{E}_{k}t}c_{k}\\\
c^{\dagger}_{k}(t)&=e^{i\mathcal{E}_{k}t}c^{\dagger}_{k}.\\\ \end{split}$ (S5)

In real space, we have that $c_{x}(t)=\sum_{y}u(x-y,t)c_{y}$. The function,
$u(x-y,t)$, is the Fourier transform of the phases $e^{-i\mathcal{E}_{k}t}$,
and controls how the correlation function grows in time. Specifically, the
envelope function takes the form,

$c_{x}(t)=\sum_{y}\underbrace{\bigg{[}\frac{1}{N}\sum_{k}e^{i(x-y)k}e^{-i\mathcal{E}_{k}t}\bigg{]}}_{\text{{\hbox{u(x-y,
t)}}}}c_{y}.$ (S6)

We now study the universal properties of this function. In the thermodynamic
limit, the mode functions take the form,

$\displaystyle\begin{split}u(x-y,t)&=\int_{-\pi}^{\pi}\frac{dk}{2\pi}e^{i(x-y)k}e^{-2i\cos(k)t}\\\
&=i^{(x-y)}\mathcal{J}_{x-y}(2t).\end{split}$ (S7)

In the previous equation, $\mathcal{J}_{\alpha}(2t)$ is a Bessel function of
the first kind of order $\alpha$. Utilizing the above result in combination
with Eq. (S6), the time-evolved operators become

$c_{x}(t)=\sum_{y}i^{(x-y)}\mathcal{J}_{x-y}(2t)c_{y}.$ (S8)

Transport is determined by determining the late-time behavior of the density-
density correlation function. First, write the time-evolved density operator,

$n_{x}(t)=\sum_{\alpha,\beta}(i)^{\alpha-\beta}\mathcal{J}_{x-\alpha}(2t)\mathcal{J}_{x-\beta}(2t)c^{\dagger}_{\alpha}c_{\beta}.$
(S9)

The correlation function from the main-text can be rewritten in the form,

$C_{xy}(t)=\frac{1}{2^{N}}\text{tr}(n_{x}(t)n_{y}(0))-\frac{1}{4}.$ (S10)

By inserting Eq. (S9) into the previous equation, we find,

$\displaystyle\begin{split}C_{x,y}(t)&=\frac{1}{2^{N}}\sum_{\alpha,\beta,\gamma,\delta}\left[i^{\alpha-\beta+\gamma-\delta}\mathcal{J}_{x-\alpha}(2t)\mathcal{J}_{x-\beta}(2t)\mathcal{J}_{y-\gamma}(0)\mathcal{J}_{y-\delta}(0)\right]\text{tr}(c^{\dagger}_{\alpha}c_{\beta}c^{\dagger}_{\gamma}c_{\delta})-\frac{1}{4}\\\
&=\frac{1}{2^{N}}\sum_{\alpha,\beta}\left[i^{\alpha-\beta}\mathcal{J}_{x-\alpha}(2t)\mathcal{J}_{x-\beta}(2t)\right]\text{tr}(c^{\dagger}_{\alpha}c_{\beta}n_{y})-\frac{1}{4}.\end{split}$
(S11)

we utilized the relation $\mathcal{J}_{\alpha-\beta}(0)=\delta_{\alpha,\beta}$
when moving between lines. In the previous equation, $C_{x,y}=0$ when
$\alpha\neq\beta$, in which the sum breaks into two separate cases: (i)
$\alpha\neq y$ and (ii) $\alpha=y$. Finally, using the identity
$\sum_{\alpha=-\infty}^{\infty}\mathcal{J}_{x-\alpha}^{2}(z)=1$, we evaluate
the correlation function to be,

$C_{x,y}(t)=\frac{1}{4}\mathcal{J}^{2}_{x-y}(2t).$ (S12)

To compare our analytical result to simulations, we utilize the non-
interacting structure of the model to simulate large systems to late time,
enabling reliable extraction hydrodynamic behavior. The following section
focuses on efficiently implementing free-fermion numerics to compute the
correlation function.

###  S2.1 Efficient Free-Fermion Simulations.

The Hamiltonian in the many-body basis is stored as a $2^{N}\times 2^{N}$
matrix, severely limiting the possible system size. However, utilizing the
non-interacting structure, we can simulate the dynamics in the single particle
basis; the Hamiltonian is a $N\times N$ matrix. While the system size is no
longer a barrier, here, because the model is Brownian, it must be diagonalized
at each time step, severely restricting the simulation time. Moreover, for the
system to remain uncorrelated in time, $dt$ must be small (at least less than
$0.1$). We now outline how to compute $C_{x,y}(t)$ numerically on the single-
particle basis. We begin with the general Hamiltonian,

$H=\sum_{ij}c_{i}^{\dagger}M_{ij}c_{j}.$ (S13)

Here, the matrix $M$ is a $N\times N$ real symmetric matrix that is
diagonalized as $M=VDV^{T}$ where $V$ are the eigenvectors satisfying
$VV^{T}=\mathbbm{1}$ and $D=\text{diag}(\lambda_{1},\ldots,\lambda_{N})$ is a
diagonal matrix. The Hamiltonian now becomes,

$\displaystyle H$
$\displaystyle=\sum_{ij}\sum_{kl}c_{i}^{\dagger}V_{ik}\lambda_{k}V_{jk}c_{j}$
$\displaystyle=\sum_{k}\lambda_{k}d_{k}^{\dagger}d_{k}.$ (S14)

We have defined the normal mode operators, $d_{k}=\sum_{j}V_{jk}c_{j}$ and
$d^{\dagger}_{k}=\sum_{i}c_{i}^{\dagger}V_{ik}$ which inherit the standard
fermionic anti-commutation relations. We now define the vectors,
$\hat{d}^{\dagger}=[d_{1}^{\dagger},\ldots,d_{N}^{\dagger}]=\hat{c}^{\dagger}V$
and $\hat{d}=[d_{1},\ldots,d_{N}]^{T}=V^{T}\hat{c}$. From the Heisenberg
evolution of the normal mode operators, we determine the original fermion
operators evolve as,

$\displaystyle\hat{c}^{\dagger}(t)=\hat{c}^{\dagger}\left(Ve^{iDt}V^{T}\right)$
$\displaystyle\hat{c}(t)=\left(Ve^{iDt}V^{T}\right)\hat{c}$ (S15)

The correlation function on the same site becomes,

$\displaystyle C_{x,x}(t)$
$\displaystyle=\frac{1}{2^{N}}\text{tr}\left(n_{x}(t)n_{x}\right)-\frac{1}{4}$
$\displaystyle=\frac{1}{2^{N}}\left(\sum_{\sigma,\rho}\sum_{\mu,\nu}(V_{\sigma\rho}e^{i\lambda_{\rho}t}V_{\rho
x}^{T})(V_{x\mu}e^{-i\lambda_{\mu}t}V_{\mu\nu}^{T})\text{tr}(c^{\dagger}_{\sigma}c_{\nu}c^{\dagger}_{x}c_{x})\right)-\frac{1}{4}$
$\displaystyle=\frac{1}{2^{N}}\left(\frac{2^{N}}{2}\sum_{\rho,\mu}V_{x\rho}V_{\rho
x}^{T}V_{x\mu}V_{\mu
x}^{T}e^{i\omega_{\rho\mu}t}+\frac{2^{N}}{4}\sum_{\rho,\sigma,\mu}V_{\sigma\rho}V_{\rho
x}^{T}V_{x\mu}V_{\mu\sigma}^{T}e^{i\omega_{\rho\mu}t}\right)-\frac{1}{4}$
$\displaystyle=\frac{1}{2}a_{x,x}(t)a^{*}_{x,x}(t)+\frac{1}{4}\sum_{\sigma,\sigma\neq
x}a_{\sigma,x}(t)a^{*}_{x,\sigma}(t)-\frac{1}{4}.$ (S16)

Where above we have defined the time-dependent amplitudes,
$a_{m,n}(t)=\sum_{k}V_{m,k}e^{i\lambda_{k}t}V^{T}_{k,n}$. Moving from the
first to the second line, we utilized that
$\text{tr}(c^{\dagger}_{\sigma}c_{\nu}c^{\dagger}_{x}c_{x})=0$ if
$\sigma\neq\nu$. In the third line we consider the two cases: $\sigma=x$ and
$\sigma\neq x$, where we used the relations $\text{tr}(n_{\sigma})=2^{N}/2$
and $\text{tr}(n_{\sigma}n_{x})=2^{N}/4$ and defined the energy difference,
$\omega_{\rho\mu}=\lambda_{\rho}-\lambda_{\mu}$. Efficient time evolution is
performed as follows:

  1. Step 1:

At each $dt$ draw $N-1$ random numbers from the Gaussian distribution with a
mean $\mu=0$ and standard deviation $\sigma=\left(\Gamma/dt\right)^{1/2}$.

  2. Step 2:

Diagonalize the $N\times N$ Hamiltonian.

  3. Step 3:

Compute the correlation function $C_{xx}(dt)$ from the eigenvalues and
eigenvectors.

  4. Step 4:

Repeat the above steps for each time-step.

  5. Step 5:

Repeat the entire procedure a designated number of disorder realizations and
average.

##  S3 Operator Equation of Motion.

After noise averaging, we now outline how to derive the effective Liouvillian
that governs the operator dynamics. We consider the Brownian couplings,
$\Gamma_{x,y}(t)$, which are time-dependent real Gaussian numbers with a mean
and variance given by,

$\displaystyle\begin{split}\mathbb{E}[\Gamma_{i,j}(t)]&=0\\\
\mathbb{E}[\Gamma_{i,j}(t)\Gamma_{l,m}(t^{\prime})]&=\Gamma_{x,y}\delta_{i,l}\delta_{j,m}\delta(t-t^{\prime}).\end{split}$
(S17)

Where the operation $\mathbb{E}[\cdot]$ denotes the Gaussian average
couplings, $\Gamma_{x,y}$ sets the energy scale of the noise, and
$\delta(t-t^{\prime})$ implies that the couplings are correlated at a single
instance in time.

We expand the flow of the Heisenberg operator up to order $dt^{2}$ because the
variance is of order $dt^{-1}$. Specifically, consider the Heisenberg
evolution of an operator $\mathcal{O}$,

$\displaystyle\begin{split}\mathcal{O}(t+dt)&=\mathcal{U}_{dt}\mathcal{O}(t)\mathcal{U}^{\dagger}_{dt}\\\
&=\mathcal{O}(t)+i[H_{t},\mathcal{O}]dt+\mathcal{L}[\mathcal{O}]dt^{2}\end{split}$
(S18)

With the unitary operator, $\mathcal{U}_{dt}=e^{-iH_{t}dt}$ and the effective
Lindblad,
$\mathcal{L}[\mathcal{O}]=H^{\dagger}_{t}\mathcal{O}H_{t}-\frac{1}{2}\\{H^{\dagger}_{t}H_{t},\mathcal{O}\\}$.
We also emphasize that the above convention differs from standard transport
theories in condensed matter that study the Green function, where the operator
equation is put into a Schrodinger equation-like form. Inserting the
Hamiltonian into the operator equation, we get,

$\displaystyle\mathcal{O}(t+dt)-\mathcal{O}(t)=\sum_{x,y}i\bigg{[}J_{x,y}[c^{\dagger}_{x}c_{y},\mathcal{O}]+\Gamma_{x,y}(t)[c^{\dagger}_{x}c_{y},\mathcal{O}]\bigg{]}dt+\mathcal{L}_{x,y}[\mathcal{O}]dt^{2}$
(S19)

The Lindblad term breaks into terms with: $J_{x,y}J_{x,y}$,
$J_{x,y}\Gamma_{x,y}(t)$, $\Gamma_{x,y}(t)\Gamma_{x,y}(t)$, moreover, only the
first and final term survive noise averaging. We now average over disorder
realizations,

$\displaystyle d\bar{\mathcal{O}_{t}}$
$\displaystyle=\sum_{x,y}i\bigg{[}J_{x,y}[c^{\dagger}_{x}c_{y},\bar{\mathcal{O}_{t}}]+\mathbb{E}[\Gamma_{x,y}(t)][c^{\dagger}_{x}c_{y},d\bar{\mathcal{O}_{t}}]\bigg{]}dt+\bigg{[}\mathbb{E}[\Gamma_{x,y}(t)\Gamma_{x,y}(t)]\mathcal{L}_{x,y}[\bar{\mathcal{O}_{t}}]+J_{x,y}J_{x,y}\mathcal{L}_{x,y}[\bar{\mathcal{O}_{t}}]\bigg{]}dt^{2}$
$\displaystyle=\sum_{x,y}iJ_{x,y}[c^{\dagger}_{x}c_{y},\bar{\mathcal{O}_{t}}]dt+\Gamma_{x,y}\mathcal{L}_{x,y}[\bar{\mathcal{O}_{t}}]dt+J_{x,y}J_{x,y}\mathcal{L}_{x,y}[\bar{\mathcal{O}_{t}}]dt^{2}.$
(S20)

Only keeping terms that are of order $dt$, leads to the noise-averaged
operator equation,

$\displaystyle\frac{d\bar{\mathcal{O}_{t}}}{dt}=\sum_{x,y}\bigg{[}iJ_{x,y}[c^{\dagger}_{x}c_{y},\bar{\mathcal{O}_{t}}]+\Gamma_{x,y}\mathcal{L}_{x,y}[\bar{\mathcal{O}_{t}}]\bigg{]}$
(S21)

Where we have defined $\mathbb{E}(d\mathcal{O})=d\bar{\mathcal{O}}$ and the
Lindblad has the form,
$\mathcal{L}_{x,y}[\bar{\mathcal{O}}]=L^{\dagger}_{x,y}\bar{\mathcal{O}}L_{x,y}-\frac{1}{2}\\{L_{x,y}^{\dagger}L_{x,y},\bar{\mathcal{O}}\\}$.
Generically, the above equation applies for long-range hopping, however, we
consider either onsite or nearest-neighbor interactions. Specifically, when
the noise is coupled to nearest-neighbor hopping ($y=x+1$) then the jump
operator is $L_{x,x+1}=c^{\dagger}_{x}c_{x+1}+c^{\dagger}_{x+1}c_{x}$,
whereas, when the noise is coupled to the density $L_{x,x}=n_{x}$. The most
general operator equation of Eq. (S21) we consider is written as follows,

$\frac{d\bar{\mathcal{O}_{t}}}{dt}=\sum_{x}\bigg{[}iJ[L_{x,x+1},\bar{\mathcal{O}_{t}}]-i\gamma
x[L_{x,x},\bar{\mathcal{O}_{t}}]+\mathcal{V}\mathcal{L}_{x,x}[\bar{\mathcal{O}_{t}}]+\Gamma\mathcal{L}_{x,x+1}[\bar{\mathcal{O}_{t}}]\bigg{]}.$
(S22)

The above equation follows from setting the diagonal (off-diagonal) elements
of the static coupling matrix to $J_{x,x}=-\gamma x$ (we focus on a linear
potential) and $J_{x,x+1}=J$, while the diagonal (off-diagonal) elements of
the noisy coupling matrix to $\Gamma_{x,x}=\mathcal{V}$ and
$\Gamma_{x,x+1}=\Gamma$.

##  S4 Details in Calculating the Coefficient Differential Equation.

In this section, we derive the differential equation for the operator
amplitudes by inserting $n_{x}(t)$ into Eq. (S21).

###  S4.1 Coherent Term from Static Hopping

We begin with the time averaged operator equation of motion with only the
static and noisy hopping contribution,

$\frac{d\bar{n}_{t}}{dt}=\sum_{x}\bigg{[}iJ[L_{x,x+1},\bar{n}_{t}]+\Gamma\mathcal{L}_{x,x+1}[\bar{n}_{t}]\bigg{]}.$
(S23)

In this section, we derive the term from the unitary dynamics,

$\displaystyle[c^{\dagger}_{x}c_{x+1},\bar{n}_{t}]$
$\displaystyle=\sum_{x}\sum_{mn}A_{m,n}\bigg{[}c^{\dagger}_{x}[c_{x+1},c^{\dagger}_{m}]c_{n}+[c^{\dagger}_{x},c^{\dagger}_{m}]c_{x+1}c_{n}+c^{\dagger}_{m}c^{\dagger}_{x}[c_{x+1},c_{n}]+c^{\dagger}_{m}[c^{\dagger}_{x},c_{n}]c_{x+1}\bigg{]}$
$\displaystyle=\sum_{x}\sum_{mn}A_{m,n}\bigg{[}c^{\dagger}_{x}\left(2c_{x+1}c^{\dagger}_{m}-\delta_{x+1,m}\right)c_{n}+2c^{\dagger}_{x}c^{\dagger}_{m}c_{x+1}c_{n}+2c^{\dagger}_{m}c^{\dagger}_{x}c_{x+1}c_{n}+c^{\dagger}_{m}\left(2c_{x}c^{\dagger}_{n}-\delta_{x,n}\right)c_{x+1}\bigg{]}$
$\displaystyle=\sum_{x}\sum_{m,n}A_{m,n}\bigg{[}2c^{\dagger}_{x}c_{x+1}c^{\dagger}_{m}c_{n}-\delta_{x+1,m}c^{\dagger}_{x}c_{n}+2c^{\dagger}_{m}c^{\dagger}_{x}c_{n}c_{x+1}-\delta_{x,n}c^{\dagger}_{m}c_{x+1}\bigg{]}$
$\displaystyle=\sum_{x}\sum_{m,n}A_{m,n}\bigg{[}\delta_{x+1,m}c^{\dagger}_{x}c_{n}-\delta_{x,n}c^{\dagger}_{m}c_{x+1}\bigg{]}.$
(S24)

In the first line, we utilized a commutator identity between four operators
and then used that $[A,B]=2AB-\\{A,B\\}$ where we then applied the anti-
commutator relation for fermion operators. For the second line, the product of
four fermion operators cancels by anti-commuting the creation operators and
then expand in the third line. Because the model is non-interacting the final
result cannot have more than a product of two fermion operators, thus we
simplify further using the anti-commutation relation
$\\{c^{\dagger}_{m},c_{x+1}\\}$ and end with the final equation. The second
commutator is found by a similar calculation or by setting $x+1\rightarrow x$
and $x\rightarrow x+1$ in the above equation,

$\displaystyle[c^{\dagger}_{x+1}c_{x},\bar{n}_{t}]$
$\displaystyle=\sum_{x}\sum_{mn}A_{m,n}\bigg{[}\delta_{x,m}c^{\dagger}_{x+1}c_{n}-\delta_{x+1,n}c^{\dagger}_{m}c_{x}\bigg{]}.$
(S25)

Moreover, we find the final equation after resolving the delta functions is,

$iJ\sum_{x}[L_{x,x+1},\bar{n}_{t}]=iJ\sum_{m,n}A_{m,n}\bigg{[}c^{\dagger}_{m-1}c_{n}-c^{\dagger}_{m}c_{n+1}+c^{\dagger}_{m+1}c_{n}-c^{\dagger}_{m}c_{n-1}\bigg{]}$
(S26)

Next, we want to shift the indices so that each operator is equivalently
$c^{\dagger}_{m}c_{n}$ which allow us to equate sides and get a differential
equation for $A_{m,n}$.

$\displaystyle A_{m,n}c^{\dagger}_{m-1}c_{n}$ $\displaystyle\rightarrow
A_{m+1,n}c^{\dagger}_{m}c_{n},\quad A_{m,n}c^{\dagger}_{m}c_{n+1}\rightarrow
A_{m,n-1}c^{\dagger}_{m}c_{n}$ $\displaystyle A_{m,n}c^{\dagger}_{m+1}c_{n}$
$\displaystyle\rightarrow A_{m-1,n}c^{\dagger}_{m}c_{n},\quad
A_{m,n}c^{\dagger}_{m}c_{n-1}\rightarrow A_{m,n+1}c^{\dagger}_{m}c_{n}$ (S27)

The final equation due to unitary dynamics is,

$iJ\sum_{x}[L_{x,x+1},\bar{n}_{t}]=iJ\sum_{mn}\bigg{[}A_{m+1,n}-A_{m,n-1}+A_{m-1,n}-A_{m,n+1}\bigg{]}c^{\dagger}_{m}c_{n}$
(S28)

In the next section, we derive the contribution from the effective Lindblad
term.

###  S4.2 Incoherent Term from Noisy Hopping.

In this section, we derive the incoherent piece that arises from the variance
of the couplings,

$\sum_{x}\mathcal{L}_{x,x}\big{[}\bar{n}_{t}\big{]}=\sum_{x}\sum_{m,n}A_{m,n}\bigg{[}L^{\dagger}_{x,x+1}c^{\dagger}_{m}c_{n}L_{x,x+1}-\frac{1}{2}\\{L^{\dagger}_{x,x+1}L_{x,x+1},c^{\dagger}_{m}c_{n}\\}\bigg{]}.$
(S29)

where $L_{x,x+1}=c^{\dagger}_{x}c_{x+1}+c^{\dagger}_{x+1}c_{x}$. After some
algebra and using fermion identities, we get the final equation,

$\displaystyle\frac{d\bar{n}_{t}}{dt}$
$\displaystyle=\sum_{mn}\bigg{[}-2A_{m,n}c^{\dagger}_{m}c_{n}+A_{m,n}\delta_{m,n-1}c^{\dagger}_{n}c_{n-1}+A_{m,n}\delta_{m,n+1}c^{\dagger}_{n}c_{n+1}+A_{m,n}\delta_{m,n}\left(n_{m+1}+n_{m-1}\right)\bigg{]}$
(S30)

Next, we want to shift the indices so that each operator is aligned with
$c^{\dagger}_{m}c_{n}$ which allow us to equate sides and get a differential
equation for $A_{m,n}$.

$\displaystyle A_{m,n}\delta_{m,n}c^{\dagger}_{m-1}c_{m-1}$
$\displaystyle\rightarrow A_{m+1,n+1}\delta_{m,n}c^{\dagger}_{m}c_{n},\quad
A_{m,n}\delta_{m,n}c^{\dagger}_{m+1}c_{m+1}\rightarrow
A_{m-1,n-1}\delta_{m,n}c^{\dagger}_{m}c_{n}$ $\displaystyle
A_{m,n}\delta_{m,n+1}c^{\dagger}_{n}c_{n+1}$ $\displaystyle\rightarrow
A_{m+1,n-1}\delta_{m+1,n}c^{\dagger}_{m}c_{n},\quad
A_{m,n}\delta_{m,n-1}c^{\dagger}_{n}c_{n-1}\rightarrow
A_{m-1,n+1}\delta_{m-1,n}c^{\dagger}_{m}c_{n}$ (S31)

The operator equation of motion becomes,

$\displaystyle\frac{d\bar{n}_{t}}{dt}$
$\displaystyle=\sum_{m,n}\bigg{[}-2A_{m,n}+\delta_{m+1,n}A_{m+1,n-1}+\delta_{m-1,n}A_{m-1,n+1}+\delta_{m,n}A_{m+1,n+1}+\delta_{m,n}A_{m-1,n-1}\bigg{]}c^{\dagger}_{m}c_{n}.$
(S32)

Moreover, we find that the contribution from the incoherent piece is as
follows,

$\frac{dA_{m,n}}{dt}=\bigg{[}-2A_{m,n}+\delta_{m+1,n}A_{m+1,n-1}+\delta_{m-1,n}A_{m-1,n+1}+\delta_{m,n}A_{m+1,n+1}+\delta_{m,n}A_{m-1,n-1}\bigg{]}$
(S33)

In the following section we will perform the mapping to the fictitious lattice
defined by the operator length and center-of-mass coordinates.

###  S4.3 Operator Equation of Motion.

Combining, the equation from the previous two sections, we arrive at the full
operator equation,

$\displaystyle\partial_{t}A_{m,n}=iJ\bigg{[}A_{m+1,n}-A_{m,n+1}+A_{m-1,n}-A_{m,n-1}\bigg{]}$
$\displaystyle+\Gamma\bigg{[}-2A_{m,n}+\delta_{m+1,n}A_{m+1,n-1}+\delta_{m-1,n}A_{m-1,n+1}+\delta_{m,n}A_{m-1,n-1}+\delta_{m,n}A_{m+1,n+1}\bigg{]}.$
(S34)

When a static potential is included in the coherent dynamics, the additional
term is $i\big{[}V_{m}-V_{n}\big{]}A_{m,n}$ which has no affect in the strong
noise limit.

###  S4.4 Mapping to $\ell$ and $\mathcal{R}$.

Spatial translation invariance is partially restored in a different coordinate
system defined by the variables $\ell=n-m$ and $\mathcal{R}=n+m$. Here $\ell$
represents the length of the operator and $\mathcal{R}$ the center-of-mass of
the operator. The coefficients under this mapping become,

$\displaystyle
A_{m+1,n-1}\overset{(n-1-m-1=\ell-2,n-1+m+1=\mathcal{R})}{\longrightarrow}A_{\ell-2,\mathcal{R}},\quad$
$\displaystyle
A_{m-1,n+1}\overset{(n+1-m+1=\ell+2,n+1+m-1=\mathcal{R})}{\longrightarrow}A_{\ell+2,\mathcal{R}}$
$\displaystyle
A_{m+1,n+1}\overset{(n+1-m-1=\ell,n+1+m+1=\mathcal{R}+2)}{\longrightarrow}A_{\ell,\mathcal{R}+2},\quad$
$\displaystyle
A_{m-1,n-1}\overset{(n-1-m+1=\ell,n-1+m-1=\mathcal{R}-2)}{\longrightarrow}A_{\ell,\mathcal{R}-2}$
$\displaystyle\delta_{m+1,n}\overset{(\ell=1)}{\longrightarrow}\delta_{\ell,1},\quad$
$\displaystyle\delta_{m-1,n}\overset{\ell=-1)}{\longrightarrow}\delta_{\ell,-1}$
(S35)

Applying this transformation to the equation above,

$\displaystyle\partial_{t}A_{\ell,\mathcal{R}}$
$\displaystyle=iJ\bigg{[}A_{\ell-1,\mathcal{R}+1}-A_{\ell+1,\mathcal{R}+1}+A_{\ell+1,\mathcal{R}-1}-A_{\ell-1,\mathcal{R}-1}\bigg{]}$
$\displaystyle+\Gamma\bigg{[}-2A_{\ell,\mathcal{R}}+\delta_{\ell,1}A_{\ell-2,\mathcal{R}}+\delta_{\ell,-1}A_{\ell+2,\mathcal{R}}+\delta_{\ell,0}\left(A_{\ell,\mathcal{R}+2}+A_{\ell,\mathcal{R}-2}\right)\bigg{]}$
(S36)

Now the above equation is simplified further by looking at the terms
$\delta_{\ell,-1}A_{\ell+2,\mathcal{R}}$ and
$\delta_{\ell,1}A_{\ell-2,\mathcal{R}}$ which are nonzero if $\ell=\pm 1$ with
the coefficients $A_{\mp 1}$. A similar argument can be made for
$\delta_{\ell,0}\left(A_{\ell,\mathcal{R}+2}+A_{\ell,\mathcal{R}-2}\right)$.
Moreover, the final equation is,

$\displaystyle\partial_{t}A_{\ell,\mathcal{R}}$
$\displaystyle=iJ\bigg{[}A_{\ell-1,\mathcal{R}+1}-A_{\ell+1,\mathcal{R}+1}+A_{\ell+1,\mathcal{R}-1}-A_{\ell-1,\mathcal{R}-1}\bigg{]}$
$\displaystyle+\Gamma\bigg{[}-2A_{\ell,\mathcal{R}}+\delta_{\ell,1}A_{-1,\mathcal{R}}+\delta_{\ell,-1}A_{1,\mathcal{R}}+\delta_{\ell,0}\left(A_{0,\mathcal{R}+2}+A_{0,\mathcal{R}-2}\right)\bigg{]}$
(S37)

Now because the delta functions are independent of $\mathcal{R}$ we can
perform a discrete Fourier transformation in this coordinate, namely,
$A_{\ell,\mathcal{R}}=\sum e^{ik\mathcal{R}}A_{\ell,k}$. The Fourier
transformed equation becomes,

$\displaystyle\partial_{t}A_{\ell,k}$
$\displaystyle=2J\sin(k)\bigg{[}A_{\ell+1,k}-A_{\ell-1,k}\bigg{]}+\Gamma\bigg{[}-2A_{\ell,k}+\delta_{\ell,1}A_{-1,k}+\delta_{\ell,-1}A_{1,k}+2\cos(2k)\delta_{\ell,0}A_{0,k}\bigg{]}$
(S38)

In the next section, we derive the contribution from a noisy potential,
specifically, $\Gamma_{x,x}\neq 0$.

###  S4.5 Incoherent Term from Noisy Potential.

In this section, we derive the incoherent piece that arises from the variance
of the couplings,

$\sum_{x}\mathcal{L}_{x}\big{[}\bar{n}_{t}\big{]}=\sum_{x}\sum_{m,n}A_{m,n}\bigg{[}n_{x}c^{\dagger}_{m}c_{n}n_{x}-\frac{1}{2}\\{n_{x},c^{\dagger}_{m}c_{n}\\}\bigg{]}.$
(S39)

As in the previous section, the final result must be a product of two
fermionic operators. We first calculate the second term in the above equation,

$\displaystyle\\{n_{x},c^{\dagger}_{m}c_{n}\\}$
$\displaystyle=c^{\dagger}_{x}c_{x}c^{\dagger}_{m}c_{n}+c^{\dagger}_{m}c_{n}c^{\dagger}_{x}c_{x}$
$\displaystyle=c^{\dagger}_{x}\left(\delta_{m,x}-c^{\dagger}_{m}c_{x}\right)c_{n}+c^{\dagger}_{m}\left(\delta_{x,n}-c^{\dagger}_{x}c_{n}\right)c_{x}$
$\displaystyle=\delta_{m,x}c^{\dagger}_{x}c_{n}+c^{\dagger}_{m}c^{\dagger}_{x}c_{x}c_{n}+c^{\dagger}_{m}c^{\dagger}_{x}c_{x}c_{n}+\delta_{n,x}c^{\dagger}_{m}c_{x}$
$\displaystyle=2c^{\dagger}_{m}c^{\dagger}_{x}c_{x}c_{n}+\delta_{m,x}c^{\dagger}_{x}c_{n}+\delta_{n,x}c^{\dagger}_{m}c_{x}.$
(S40)

We put the fermion operators in normal order to simplify all these
expressions. Therefore, in the first equation above, we use the anti-
commutator to move the creation operators to the left and then do algebra on
the subsequent lines. Combing everything, we have that,

$-\frac{1}{2}\\{n_{x},\bar{n}_{t}\\}=-\sum_{x}\sum_{m,n}A_{m,n}c^{\dagger}_{m}c^{\dagger}_{x}c_{x}c_{n}-\sum_{mn}A_{m,n}c^{\dagger}_{m}c_{n}.$
(S41)

We resolved the delta functions in the second term to eliminate the sum over
x. Note that if $A_{m,n}=\delta_{x,m}\delta_{x,n}$ we get that the above
equation is $-n_{x}$ which we expect from
$-\frac{1}{2}\\{n_{x},n_{x}\\}=-n_{x}$. Another check is if
$A_{m,n}=\delta_{y,m}\delta_{y,n}$ then we find $-n_{y}n_{x}$ as expected. We
now compute the first term, which is a product of six fermion operators,

$\displaystyle c^{\dagger}_{x}c_{x}c^{\dagger}_{m}c_{n}c^{\dagger}_{x}c_{x}$
$\displaystyle=c^{\dagger}_{x}\left(\delta_{mx}-c^{\dagger}_{m}c_{x}\right)c_{n}c^{\dagger}_{x}c_{x}$
$\displaystyle=\delta_{m,x}c^{\dagger}_{x}c_{n}c^{\dagger}_{x}c_{x}-c^{\dagger}_{x}c^{\dagger}_{m}c_{x}c_{n}c^{\dagger}_{x}c_{x}$
$\displaystyle=\delta_{m,x}c^{\dagger}_{x}c_{n}c^{\dagger}_{x}c_{x}-\delta_{n,x}c^{\dagger}_{x}c^{\dagger}_{m}c_{x}c_{x}+c^{\dagger}_{x}c^{\dagger}_{m}c_{x}c^{\dagger}_{x}c_{n}c_{x}$
$\displaystyle=\delta_{m,x}c^{\dagger}_{x}c_{n}c^{\dagger}_{x}c_{x}+c^{\dagger}_{x}c^{\dagger}_{m}c_{n}c_{x}$
(S42)

In total we have that,

$\sum_{x}\sum_{m,n}A_{m,n}c^{\dagger}_{x}c_{x}c^{\dagger}_{m}c_{n}c^{\dagger}_{x}c_{x}=\sum_{x}\sum_{mn}A_{m,n}\bigg{[}\delta_{m,x}c^{\dagger}_{x}c_{n}c^{\dagger}_{x}c_{x}+c^{\dagger}_{x}c^{\dagger}_{m}c_{n}c_{x}\bigg{]}$
(S43)

We now combine the previous two results,

$\displaystyle\sum_{x}\mathcal{L}_{x}\big{[}\bar{n}_{t}\big{]}$
$\displaystyle=\sum_{x}\sum_{m,n}A_{m,n}c^{\dagger}_{x}c_{x}c^{\dagger}_{m}c_{n}c^{\dagger}_{x}c_{x}-\frac{1}{2}A_{m,n}\\{n_{x},c^{\dagger}_{m}c_{n}\\}$
$\displaystyle=\sum_{x}\sum_{m,n}A_{m,n}\bigg{[}\delta_{m,x}c^{\dagger}_{x}c_{n}c^{\dagger}_{x}c_{x}+c^{\dagger}_{x}c^{\dagger}_{m}c_{n}c_{x}\bigg{]}-\sum_{x}\sum_{m,n}A_{m,n}c^{\dagger}_{m}c^{\dagger}_{x}c_{x}c_{n}-\sum_{m,n}A_{m,n}c^{\dagger}_{m}c_{n}$
$\displaystyle=\sum_{x}\sum_{m,n}A_{m,n}\delta_{m,x}c^{\dagger}_{x}c_{n}c^{\dagger}_{x}c_{x}-\sum_{m,n}A_{m,n}c^{\dagger}_{m}c_{n}$
$\displaystyle=\sum_{m,n}A_{m,n}c^{\dagger}_{m}c_{n}c^{\dagger}_{m}c_{m}-\sum_{m,n}A_{m,n}c^{\dagger}_{m}c_{n}$
(S44)

There are a few checks we can perform on the above equation. First, if
$A_{m,n}=\delta_{x,m}\delta_{x,n}$ or $A_{m,n}=\delta_{y,m}\delta_{y,n}$ the
above result vanishes. When the Lindblad acts on a number operator, i.e.,
$\mathcal{L}[c^{\dagger}_{\alpha}c_{\alpha}]$, it annihilates it. The final
simplification is the normal ordering of the first term above,

$\sum_{x}\mathcal{L}_{x}\big{[}\bar{n}_{t}\big{]}=\sum_{m,n}A_{m,n}\left(\delta_{n,m}c^{\dagger}_{m}c_{m}-c^{\dagger}_{m}c_{n}\right)$
(S45)

Notice above that if the operator is on the same site, i.e., $m=n$, then we
get zero, and if $m\neq n$, we get $-\bar{n}_{t}$. Thus we can rewrite the
above equation as,

$\sum_{x}\mathcal{L}_{x}\big{[}\bar{n}_{t}\big{]}=\sum_{mn}A_{m,n}\left(\delta_{nm}-1\right)c^{\dagger}_{m}c_{n}$
(S46)

Moreover, the final differential equation becomes,

$\frac{d\bar{n}_{t}}{dt}=\mathcal{V}\left(\delta_{nm}-1\right)A_{m,n}$ (S47)

###  S4.6 Mapping to $\ell$ and $\mathcal{R}$

In this section, we rewrite the final equation of the previous section in a
different coordinate system defined by the variables $\ell=n-m$ and
$\mathcal{R}=n+m$. Here $\ell$ represents the length of the operator, and
$\mathcal{R}$ is the center of mass of the operator. The coefficients under
this mapping become,

$\displaystyle
A_{m+1,n}\overset{(n-m-1=\ell-1,n+m+1=\mathcal{R}+1)}{\longrightarrow}A_{\ell-1,\mathcal{R}+1},\quad$
$\displaystyle
A_{m,n-1}\overset{(n-1-m=\ell-1,n-1+m=\mathcal{R}-1)}{\longrightarrow}A_{\ell-1,\mathcal{R}-1}$
$\displaystyle
A_{m-1,n}\overset{(n-m+1=\ell+1,n+m-1=\mathcal{R}-1)}{\longrightarrow}A_{\ell+1,\mathcal{R}-1},\quad$
$\displaystyle
A_{m,n+1}\overset{(n+1-m=\ell+1,n+1-m=\mathcal{R}+1)}{\longrightarrow}A_{\ell+1,\mathcal{R}+1}$
$\displaystyle
A_{m,n}\overset{(n-m=\ell,n+m=\mathcal{R})}{\longrightarrow}A_{\ell,\mathcal{R}},\quad$
$\displaystyle\delta_{m,n}\overset{(n-m=0)}{\longrightarrow}\delta_{\ell,0}.$
(S48)

Moreover, the differential equation becomes,

$\frac{dA_{\ell,\mathcal{R}}}{dt}=\mathcal{V}A_{\ell,\mathcal{R}}\left(\delta_{\ell,0}-1\right)$
(S49)

The advantage of moving into this coordinate system is that the delta function
is independent of $\mathcal{R}$; thus, we can perform a discrete Fourier
transformation $A_{\ell,\mathcal{R}}=\sum_{k}e^{ik\mathcal{R}}A_{\ell,k}$
where $k$ is the momentum of the Brillouin zone. The Fourier transformed
equation becomes,

$\frac{dA_{\ell,k}}{dt}=\mathcal{V}A_{\ell,k}\left(\delta_{\ell,0}-1\right).$
(S50)

With all these pieces in hand, the differential equation from Eq. (S22) takes
the form,

$\frac{dA_{\ell}}{dt}=t_{k}\big{[}A_{\ell+1}-A_{\ell-1}\big{]}+\Gamma\bigg{[}-2A_{\ell,k}+\delta_{\ell,1}A_{-1,k}+\delta_{\ell,-1}A_{1,k}+2\cos(2k)\delta_{\ell,0}A_{0,k}\bigg{]}+i\gamma\ell
A_{\ell}+\mathcal{V}\big{[}\delta_{\ell,0}-1\big{]}A_{\ell}.$ (S51)

The following sections will focus on various limits of the above equation and
solve for the diffusion constant exactly in the small momentum limit.

##  S5 Diffusion Constant for Static Hopping with Bond Noise.

Figure S1: Static and Noisy Hopping Characteristics. (a) Noisy-hopping
operator dynamics. (b) Bound state energy for the case of static and noisy
hopping, i.e., $\mathcal{V}=0$ and $\gamma=0$. The solid curves are the
analytical result, Eq. (S5) in the small momentum limit, while the black
crosses are from diagonalizing the eigenvalue equations, Eq. (S52). (c)-(f)
Real part of the eigenvalue spectrum of Eq. (S52) for different noisy and
static hopping strengths. Regardless of the noise strength, there is always a
bound state at small momentum. As momentum increases, the bound state vanishes
in a continuum of bound states at $-2\Gamma$ represented by the red line.
Increasing $k$ further reveals that the bound state reemerges at an
exceptional point. Before the exceptional point, the cubic equation admits two
physical solutions with real energies, which then coalesce and become a
conjugate pair of complex energies. Interestingly, near the strong drive limit
there is a reentrant-like point [see (d)] where the degeneracy is re-broken,
and the purely real energies emerge once more. Parameters: (a) $N=400$,
$dt=0.05$, $\ell_{\text{max}}=N/2$, (b) $N=400$, and (c)-(f) $N=400$,
$\ell_{\text{max}}=N/2$.

In this section, we give an explicit solution to the bound state energy in the
case of finite noise with no noisy or static potential, i.e., $\mathcal{V}=0$
and $\gamma=0$ in Eq. (S51). We begin with the eigenvalue equation,

$\mathcal{E}_{q}A_{\ell}=t_{k}\bigg{[}A_{\ell+1}-A_{\ell-1}\bigg{]}+\Gamma\bigg{[}-2A_{\ell}+\delta_{\ell,1}A_{-1}+\delta_{\ell,-1}A_{1}+2\cos(2k)\delta_{\ell,0}A_{0}\bigg{]}.$
(S52)

We have dropped the subscript $k$ since we are solving these for a fixed
momentum and $t_{k}=2J\sin(k)$ is the effective hopping strength. For finite
noise, the above equation is non-hermitian and effectively describes the
dynamics of a one-dimensional hopping model. The above equation splits into
four distinct eigenvalue equations,

$\displaystyle\mathcal{E}_{q}A_{0}$
$\displaystyle=t_{k}\left(A_{1}-A_{-1}\right)-4\Gamma\sin^{2}(k)A_{0},\quad\mathcal{E}_{q}A_{1}=t_{k}\left(A_{2}-A_{0}\right)-2\Gamma
A_{1}+\Gamma A_{-1}$ $\displaystyle\mathcal{E}_{q}A_{-1}$
$\displaystyle=t_{k}\left(A_{0}-A_{-2}\right)-2\Gamma A_{-1}+\Gamma
A_{1},\quad\mathcal{E}_{q}A_{\ell}=t_{k}\left(A_{\ell+1}-A_{\ell-1}\right)-2\Gamma
A_{\ell},$ (S53)

The final equation describes the bulk dynamics and sets the energy of the
problem, while the first three equations are boundary conditions that
constrain the allowed values that $q$ can take. For $J=0$, the above equations
decouple. The eigenvalue equation for $A_{0}$ maps to the standard diffusion
equation, and the entire spectrum forms horizontal lines at
$\mathcal{E}=-1,-2,-3$ where the fourth value is the diffusive mode $-4k^{2}$
visualized in Fig. (S1)(c). On the other hand, for $\Gamma=0$, the model is
free with ballistic transport. We now attempt to solve for the diffusion
constant for generic noise strength, $\Gamma$. First, we solve the bulk
equation with the solution, $A_{\ell}=Ae^{q\ell}+Be^{-q\ell+i\pi\ell}$,

$\mathcal{E}_{q}=4J\sin(k)\sinh(q)-2\Gamma$ (S54)

Generically, $q$ depends on momentum $k$ which we determine through the
equation $A_{0}$. First, if $A_{1}$ and $A_{-1}$ are known, then the entire
set of coefficients is determined through the above equations; thus we use the
ansatz,

$A_{\ell}=\begin{cases}A_{-1}e^{q(1+\ell)}&\text{if}\ \ell\leq-1\\\ \\\
-A_{1}e^{q(1-\ell)+i\pi\ell}&\text{if}\ \ell\geq 1\end{cases}$ (S55)

The three remaining equations solve for the variables $A_{1}/A_{-1}$,
$A_{0}/A_{-1}$, and $q$. The middle two equations solve for the ratios, and
the equation for $A_{0}$ restricts the values $q$ takes. First, rewrite the
middle two equations,

$\displaystyle\bigg{[}\mathcal{E}_{q}+2Je^{-q}\sin(k)+2\Gamma\bigg{]}\frac{A_{1}}{A_{-1}}+2J\sin(k)\frac{A_{0}}{A_{-1}}=\Gamma$
$\displaystyle\bigg{[}\mathcal{E}_{q}+2Je^{-q}\sin(k)+2\Gamma\bigg{]}-2J\sin(k)\frac{A_{0}}{A_{-1}}=\Gamma\frac{A_{1}}{A_{-1}}$
(S56)

Solving these equations yields,

$\frac{A_{1}}{A_{-1}}=-1,\quad\frac{A_{0}}{A_{-1}}=e^{q}+\frac{\Gamma}{2J\sin(k)}$
(S57)

Here $A_{1}$ and $A_{-1}$ are determined through normalizing $A_{\ell}$. The
value of $A_{0}$ is determined through the second equation. The ratios above,
combined with the final eigenvalue equation for $A_{0}$ leads to the
constraint,

$\bigg{[}\mathcal{E}_{q}+4\Gamma\sin^{2}(k)\bigg{]}\bigg{[}e^{q}+\frac{\Gamma}{2J\sin(k)}\bigg{]}=-4J\sin(k)$
(S58)

The above equation is a cubic equation for $q$, which can be solved exactly.
At small momentum, all three solutions to the cubic equation are physical and
admit purely real energies; one corresponds to the largest real eigenvalue,
the bound state energy. We solve the cubic equation and expand the bound state
energy around $k=0$ giving $q$,

$e^{q}=-\frac{\Gamma}{Jk}-\bigg{[}\frac{J}{3\Gamma}+\frac{11\Gamma}{6J}\bigg{]}k$
(S59)

Inserting the above solution into the equation for the energy and keeping the
lowest order in $k$, we find,

$\displaystyle\mathcal{E}_{q}$
$\displaystyle=2kJ\left(e^{q}-e^{-q}\right)-2\Gamma\approx-\bigg{[}\frac{8}{3}\frac{J^{2}}{\Gamma}-4\Gamma\bigg{]}k^{2}$
$\displaystyle\mathcal{D}$
$\displaystyle=\bigg{[}\frac{8}{3}\frac{J^{2}}{\Gamma}-4\Gamma\bigg{]}$ (S60)

Moreover, regardless of the strength of the noisy hopping, the late-time
hydrodynamic tail remains diffusive. In Fig. S1(c)-(f), we diagonalize Eq.
(S5) for different noise and static strengths and plot the real part of the
spectrum. For all strengths, the diffusive mode is present at small $k$, which
we fit with the analytical result Eq. (S5) [see yellow curve]. As momentum
increases, the bound state undergoes a phase transition by entering a
continuum of scattering states at $-2\Gamma$, indicated by the red line in
Fig. S1(c)-(f). Interestingly as $k$ increases further, the bound state
reemerges at an exceptional like point. At this point, the two physical
solutions to the cubic equation, which give purely real energies, collide and
become a pair of conjugate energies. Interestingly, there is a higher-order
exceptional-like point near the strong noise limit where the solutions
transition back to purely real non-degenerate energies; see Fig. S1(d). We
also emphasize that an alternative solution to the above eigenvalue equations
is to assume $A_{0}=0$, in which case you can determine analytically the form
of the branch beginning at $-\Gamma$ in Fig. S1(c)-(f).

##  S6 Diffusion Constant for Static Hopping and Onsite Noise.

Figure S2: Static Hopping and Noisy Potential Characteristics. (a) Operator
dynamics for static hopping with a noisy potential. (b) The bound state energy
for the noisy potential with static hopping. The dotted orange line is from
diagonalizing the eigenvalue equation with the other curves representing
different stages in approximating the low momentum scaling. (c) Real part
eigenvalue spectrum for the noisy potential, where there is a single bound
state corresponding to the diffusive mode with energy given by Eq. (S6). (d)
The coefficients $A_{\ell}$ from Eq. (S66) for different $k$ compared to
diagonalizing of Eq. (S61). Parameters: (a) $N=400$, $dt=0.05$,
$\ell_{\text{max}}=N/2$, (b) $N=4000$, $\mathcal{V}/J=1$, and (c) $N=4000$ (d)
$N=4000$, $\mathcal{V}/J=1$.

In this section, we solve Eq. (S51) with $\Gamma=\gamma=0$. We begin with the
eigenvalue equation,

$\mathcal{E}_{q}A_{\ell,k}=t_{k}\left(A_{\ell+1,k}-A_{\ell-1,k}\right)+\mathcal{V}\left(\delta_{\ell,0}-1\right)A_{\ell,k}$
(S61)

We now split the above equation into two different equations,

$\displaystyle\mathcal{E}_{q}A_{0,k}$
$\displaystyle=t_{k}\left(A_{1,k}-A_{-1,k}\right)$
$\displaystyle\mathcal{E}_{q}A_{\ell,k}$
$\displaystyle=t_{k}\left(A_{\ell+1,k}-A_{\ell-1,k}\right)-\mathcal{V}A_{\ell,k},$
(S62)

where $\ell$ represents the entire line without zero. Using the second
equation, we can determine the energy by assuming the ansatz,
$A_{\ell}=Ae^{q\ell}+Be^{-q\ell+i\pi\ell}$ leading to the relation,

$\displaystyle(\mathcal{V}+\mathcal{E}_{q})\left(Ae^{q\ell}+Be^{-q\ell+i\pi\ell}\right)$
$\displaystyle=2J\sin(k)\bigg{[}Ae^{q\ell}(e^{-q}-e^{q})+Be^{-q\ell+i\pi\ell}(e^{q+i\pi}-e^{-q+i\pi})\bigg{]}$
$\displaystyle=4J\sin(k)\sinh(q)\left(Ae^{q\ell}+Be^{-q\ell+i\pi\ell}\right).$
(S63)

Moreover, the energy is determined as

$\mathcal{E}_{q}=4J\sin(k)\sinh(q)-\mathcal{V}$ (S64)

The next step is to determine the values $q$ can take. First, we write the
solution to the second equation in both regions,

$A_{\ell}=\begin{cases}Ae^{q\ell}+Be^{-q\ell+i\pi\ell}&\text{if}\ \ell<0\\\
Ce^{q\ell}+De^{-q\ell+i\pi\ell}&\text{if}\ \ell>0\end{cases}$ (S65)

In the region $\ell\leq 0$ the second term is unbounded as
$\ell\rightarrow-\infty$ thereby we set $B=0$. Similarly, for $\ell\geq 0$ the
first term is unbounded as $\ell\rightarrow\infty$ and we therefore also set
$C=0$. A further constraint is made at $\ell=0$ which sets $A=D$ and the
general solution becomes,

$A_{\ell}=\begin{cases}Ae^{q\ell}&\text{if}\ \ell\leq 0\\\
Ae^{-q\ell+i\pi\ell}&\text{if}\ \ell\geq 0.\end{cases}$ (S66)

With this solution, we can determine $A_{1}$ and $A_{-1}$ and inset them into
the first equation above for $A_{0}$ which leads to the constraint on $q$,

$\displaystyle\mathcal{E}_{q}A$
$\displaystyle=2J\sin(k)\big{[}-Ae^{-q}-Ae^{-q}\big{]}$
$\displaystyle=-4J\sin(k)e^{-q}$ $\displaystyle=4J\sin(k)\sinh(q)-\mathcal{V}$

From this, we find that $q$ is

$q=\cosh^{-1}\bigg{[}\frac{\mathcal{V}}{4J\sin(k)}\bigg{]}$ (S68)

If now assume to be near the bottom of the band such that $\sin(k)\approx k$
then we find the energy and diffusion constant as follows,

$\displaystyle\mathcal{E}_{q}$
$\displaystyle=4J\sin(k)\bigg{[}\frac{\sqrt{\mathcal{V}^{2}-16J^{2}\sin^{2}(k)}}{4J\sin(k)}\bigg{]}-\mathcal{V}\approx\frac{-8J^{2}k^{2}}{\mathcal{V}}$
$\displaystyle\mathcal{D}$ $\displaystyle=\frac{8J^{2}}{\mathcal{V}}$ (S69)

Where because $k$ is small we binomially expand the square root and find that
the bound state energy scales as $k^{2}$ which indicates diffusion at late
times.

As $k$ increases, $q$ decreases, becoming zero at a critical $k$ given by
$4J\sin k=\mathcal{V}$ for $\mathcal{V}<4J$, at which the energy equals
-$\mathcal{V}$ and the bound state enters the continuum of scattering states
and vanishes. This transition manifests in the long-time behavior of the time-
dependent correlation function by displaying increased oscillations before
crossing over to diffusion; see Fig. S2(a) for numerical simulations of Eq.
(S61) for different $\mathcal{V}/J$. In Fig. S2(d), we compare our analytical
result for $A_{\ell}$ to diagonalization of the eigenvalue equations, when
$k=0.252$, $q\approx 0$ and increasing $k$ anymore cause the solution to
crossover to a scattering state. We also emphasize that the diffusion constant
vanishes as $\mathcal{V}\rightarrow\infty$ indicating emergent localization.

##  S7 Diffusion Constant for Static Hopping with Bond and Onsite Noise.

In this section, we give an explicit solution to the bound state energy in the
case of finite noise with no static potential, i.e., $\gamma=0$ in Eq. (S51).
We begin with the eigenvalue equation,

$\mathcal{E}_{q}A_{\ell}=t_{k}\bigg{[}A_{\ell+1}-A_{\ell-1}\bigg{]}+\Gamma\bigg{[}-2A_{\ell}+\delta_{\ell,1}A_{-1}+\delta_{\ell,-1}A_{1}+2\cos(2k)\delta_{\ell,0}A_{0}\bigg{]}+\mathcal{V}\bigg{[}\delta_{\ell,0}-1\bigg{]}A_{\ell}.$
(S70)

We have dropped the subscript $k$ since we are solving these for a fixed
momentum and $t_{k}=2J\sin(k)$ is the effective hopping strength. For finite
noise, the above equation is non-hermitian and effectively describes the
dynamics of a one-dimensional hopping model. The above equation splits into
four distinct eigenvalue equations,

$\displaystyle\mathcal{E}_{q}A_{0}$
$\displaystyle=t_{k}\left(A_{1}-A_{-1}\right)-4\Gamma\sin^{2}(k)A_{0},\quad\mathcal{E}_{q}A_{1}=t_{k}\left(A_{2}-A_{0}\right)-\big{[}2\Gamma+\mathcal{V}\big{]}A_{1}+\Gamma
A_{-1}$ $\displaystyle\mathcal{E}_{q}A_{-1}$
$\displaystyle=t_{k}\left(A_{0}-A_{-2}\right)-\big{[}2\Gamma+\mathcal{V}\big{]}A_{-1}+\Gamma
A_{1},\quad\mathcal{E}_{q}A_{\ell}=t_{k}\left(A_{\ell+1}-A_{\ell-1}\right)-\big{[}2\Gamma+\mathcal{V}\big{]}A_{\ell}.$
(S71)

The final equation describes the bulk dynamics and sets the energy of the
problem, while the first three equations are boundary conditions that
constrain the allowed values that $q$ can take. The three remaining equations
solve for the variables $A_{1}/A_{-1}$, $A_{0}/A_{-1}$, and $q$. Using the
same ansatz as the previous two sections, leads to the energy,

$\mathcal{E}_{q}=4J\sin(k)\sinh(q)-\big{[}\mathcal{V}+2\Gamma\big{]}$ (S72)

The middle two equations solve for the ratios, and the equation for $A_{0}$
restricts the values $q$ takes. First, rewrite the middle two equations,

$\displaystyle\bigg{[}\mathcal{E}_{q}+2Je^{-q}\sin(k)+2\Gamma+\mathcal{V}\bigg{]}\frac{A_{1}}{A_{-1}}+2J\sin(k)\frac{A_{0}}{A_{-1}}=\Gamma$
(S73)
$\displaystyle\bigg{[}\mathcal{E}_{q}+2Je^{-q}\sin(k)+2\Gamma+\mathcal{V}\bigg{]}-2J\sin(k)\frac{A_{0}}{A_{-1}}=\Gamma\frac{A_{1}}{A_{-1}}$
(S74)

Solving these equations yields the same consistency equations as above,

$\frac{A_{1}}{A_{-1}}=-1,\quad\frac{A_{0}}{A_{-1}}=e^{q}+\frac{\Gamma}{2J\sin(k)}$
(S75)

Here $A_{1}$ and $A_{-1}$ are determined through normalizing $A_{\ell}$. The
value of $A_{0}$ is determined through the second equation. The ratios above,
combined with the final eigenvalue equation for $A_{0}$ leads to the
constraint,

$\bigg{[}\mathcal{E}_{q}+4\Gamma\sin^{2}(k)\bigg{]}\bigg{[}e^{q}+\frac{\Gamma}{2J\sin(k)}\bigg{]}=-4J\sin(k)$
(S76)

The above equation is a cubic equation for $q$, which can be solved exactly.
We find that the bound state energy has the scaling,

$\displaystyle\mathcal{E}_{q}$
$\displaystyle=-4\bigg{[}\Gamma+\frac{2J^{2}}{\mathcal{V}+3\Gamma}\bigg{]}k^{2}$
$\displaystyle\mathcal{D}$
$\displaystyle=4\bigg{[}\Gamma+\frac{2J^{2}}{\mathcal{V}+3\Gamma}\bigg{]}$
(S77)

The above solution in the limit $\mathcal{V}=0$ becomes Eq. (S5) for the case
of static hopping with bond noise, while in the limit $\Gamma=0$ the bound
state energy becomes Eq. (S6) for the case of static hopping with onsite
noise. Notice, that the diffusion constant will generically have non-monotonic
behavior where when $\Gamma$ is increasing there will be a regime of
suppressed diffusion from the second term, while when $\Gamma$ continues to
increase the first term will dominate leading to monotonic growth to the
asymptotic value $4\Gamma$.

##  S8 Diffusion Constant for Linear Potential with Bond Noise.

Figure S3: Coefficient Ratios for Linear Potential. (a) Coefficient ratio
$A^{R}_{0}/A^{R}_{-1}$, which governs the bound state energy scaling, as a
function of momentum $k$ for different potential strengths and fixed
$\Gamma/J=1$. The black crosses are from numerical diagonalization of Eq.
(S78), while the curves are analytical results. (b) Same as (a) with fixed
potential strength $\gamma=4$ and different noise strengths. It is important
to note that in both (a) and (b) the final points at $k=0.2$ are
$A^{R}_{0}/A^{R}_{-1}>100$. (c) $A^{R}_{2}/A^{R}_{-1}$, (d)
$A^{R}_{-2}/A^{R}_{-1}$, (e) $A^{I}_{-2}/A^{R}_{-1}$, (f)
$A^{I}_{-2}/A^{R}_{1}$ for different potential strengths. To determine the
analytical form of the bound state energy Eq. (S8.0.1) we neglected all
$A^{R(I)}_{\pm 2}$ terms in the operator equations by arguing they are small
compared to $A^{R}_{0}/A^{R}_{-1}$. These plots provide a numerical
justification of this approximation, where clearly the ratios (c)-(f) are
significantly smaller over the range of momentum and give a negligent
contribution to the energy. Parameters: $N=400$

We now consider the case where we include a static potential in the operator
equation. Specifically, we consider a linear potential whereby the system
undergoes Stark localization, and transport is absent in the clean limit. The
eigenvalue equation has the form,

$\mathcal{E}_{q}A_{\ell}=2J\sin(k)\bigg{[}A_{\ell+1}-A_{\ell-1}\bigg{]}+\Gamma\bigg{[}-2A_{\ell}+\delta_{\ell,1}A_{-1}+\delta_{\ell,-1}A_{1}+2\cos(2k)\delta_{\ell,0}A_{0}\bigg{]}+i\gamma\ell
A_{\ell}.$ (S78)

We note that with this potential the coordinate $\mathcal{R}$ still has
spatially translation invariance, however the explicit dependence on $\ell$
prevents the ansatz used above. Moreover, we split Eq. (S78) into boundary and
bulk equations,

$\displaystyle\mathcal{E}_{q}A_{0}$
$\displaystyle=t_{k}\left(A_{1}-A_{-1}\right)-4\Gamma\sin^{2}(k)A_{0},\quad\mathcal{E}_{q}A_{1}=t_{k}\left(A_{2}-A_{0}\right)-2\Gamma
A_{1}+\Gamma A_{-1}+i\gamma A_{1}$ $\displaystyle\mathcal{E}_{q}A_{-1}$
$\displaystyle=t_{k}\left(A_{0}-A_{-2}\right)-2\Gamma A_{-1}+\Gamma
A_{1}-i\gamma
A_{-1},\quad\mathcal{E}_{q}A_{\ell}=t_{k}\left(A_{\ell+1}-A_{\ell-1}\right)-2\Gamma
A_{\ell}+i\gamma\ell A_{\ell}.$ (S79)

The bulk equation without the boundary conditions is anti-hermitian;
therefore, the energies are imaginary with a fixed real part, i.e.,
$\mathcal{E}_{q}=-2\Gamma+i\gamma q$ where
$q\in\\{-\ell_{\text{max}},\ell_{\text{max}}\\}$. The energy becomes purely
imaginary in the noise-free limit, forming a Wannier-Stark ladder leading to
non-thermal Bloch oscillations. Including the boundary conditions
$A_{0},A_{1},A_{-1}$ at finite noise has the effect of turning the set of
equations non-hermitian, causing an eigenstate to emerge from the ladder with
purely real energy, which is the bound state that governs the hydrodynamic
behavior. The bulk equation is the modified Bessel function recurrence
relation, solved with the following,

$A_{\ell}=\begin{cases}A\mathcal{I}_{\nu_{-}}(-2it_{k}/\gamma)&\text{if}\
\ell\leq-1\\\ \\\ B\mathcal{I}_{\nu_{+}}(-2it_{k}/\gamma)&\text{if}\ \ell\geq
1.\end{cases}$ (S80)

where $\mathcal{I}_{\nu_{\pm}}(iz)$ is the modified Bessel function of the
first kind of order
$\nu_{\pm}=\frac{i(\mathcal{E}_{q}+2\Gamma)}{\gamma}\pm\ell$. The difference
in the sign for the Bessel function order ensures that the solution decays to
zero as $\ell\rightarrow\pm\infty$; otherwise, it is generically a poorly
behaved function on the left half.

####  S8.0.1 Approximate solution to Bound State Energy.

Unlike the case without a potential, here, determining the analytical form of
the bound state energy is tricky because the solution to the recurrence
relation explicitly depends on the energy through the Bessel function order.
In the previous section, determining the energy came down to finding a
particular ratio of the coefficients. In this section, we construct a solution
to the bound state energy by making a few approximations that make the
calculation more tenable. To begin, we expand the three equations for $A_{0}$
and $A_{\pm 1}$ into six equations by breaking them into real and imaginary
parts, i.e., $A_{\ell}=A^{\mathbb{R}}_{\ell}+iA^{I}_{\ell}$. The following set
is as follows,

$\displaystyle\mathcal{E}_{q}A^{\mathbb{R}}_{0}$
$\displaystyle=t_{k}\left(A^{\mathbb{R}}_{1}-A^{\mathbb{R}}_{-1}\right)-4\Gamma\sin^{2}(k)A^{\mathbb{R}}_{0},\quad\mathcal{E}_{q}A^{\mathbb{R}}_{1}=t_{k}\left(A^{\mathbb{R}}_{2}-A^{\mathbb{R}}_{0}\right)-2\Gamma
A^{\mathbb{R}}_{1}+\Gamma A^{\mathbb{R}}_{-1}-\gamma A^{I}_{1}$
$\displaystyle\mathcal{E}_{q}A^{\mathbb{R}}_{-1}$
$\displaystyle=t_{k}\left(A^{\mathbb{R}}_{0}-A^{\mathbb{R}}_{-2}\right)-2\Gamma
A^{\mathbb{R}}_{-1}+\Gamma A^{\mathbb{R}}_{1}+\gamma
A^{I}_{-1},\quad\mathcal{E}_{q}A^{I}_{0}=t_{k}\left(A^{I}_{1}-A^{I}_{-1}\right)-4\Gamma\sin^{2}(k)A^{I}_{0}$
$\displaystyle\mathcal{E}_{q}A^{I}_{1}$
$\displaystyle=t_{k}\left(A^{I}_{2}-A^{I}_{0}\right)-2\Gamma A^{I}_{1}+\Gamma
A^{I}_{-1}+\gamma
A^{\mathbb{R}}_{1},\quad\mathcal{E}_{q}A^{I}_{-1}=t_{k}\left(A^{I}_{0}-A^{I}_{-2}\right)-2\Gamma
A^{I}_{-1}+\Gamma A^{I}_{1}-\gamma A^{\mathbb{R}}_{-1}.$ (S81)

From numerics we know the following features:
$A^{\mathbb{R}}_{1}/A^{\mathbb{R}}_{-1}=-1$, $A^{I}_{1}/A^{I}_{-1}=1$, and
$A^{I}_{0}=0$. Moreover, the bound state energy is solely determined by the
ratio $A^{\mathbb{R}}_{0}/A^{\mathbb{R}}_{-1}$,

$\bigg{[}\mathcal{E}_{q}+4\Gamma\sin^{2}(k)\bigg{]}\frac{A^{\mathbb{R}}_{0}}{A^{\mathbb{R}}_{-1}}=-2t_{k}$
(S82)

To solve the above set of equations, we first study the ratio
$A^{\mathbb{R}}_{0}/A^{\mathbb{R}}_{-1}$ as a function of momentum for
different $\gamma$ and $\theta$ in Fig. S3(a) and (b) where the functional
form roughly decays as $\sim 1/k$. We then numerically compare the ratios
including $A^{\mathbb{R}(I)}_{\pm 2}$ as a function of momentum in Fig.
S3(c)-(f) and find that they are small compared to
$A^{\mathbb{R}}_{0}/A^{\mathbb{R}}_{-1}$ which has a final value greater than
one-hundred. Moreover, in the equations above we set all
$A^{\mathbb{R}(I)}_{\pm 2}=0$. Next, we use the imaginary equations to
determine the ratio $A^{I}_{1}/A^{\mathbb{R}}_{-1}$ and
$A^{I}_{-1}/A^{\mathbb{R}}_{-1}$, which we find to be,

$\frac{A^{I}_{-1}}{A^{\mathbb{R}}_{-1}}=\frac{-\gamma}{\mathcal{E}_{q}+\Gamma}.$
(S83)

Inserting this relation into the equations for
$A^{\mathbb{R}}_{1}$,$A^{\mathbb{R}}_{-1}$ we find that,

$\frac{A^{\mathbb{R}}_{0}}{A^{\mathbb{R}}_{-1}}=\frac{\gamma^{2}+3\Gamma^{2}}{2J\Gamma\sin(k)}$
(S84)

Inserting into Eq. (S82) will give a cubic equation in $\mathcal{E}_{q}$ just
as in the previous section. Fig. S3(a) and (b) we plot the above analytical
result [solid curves] neglecting terms with $\mathcal{E}_{q}$ and find perfect
agreement with numerical diagonalization [black crosses]. With this
approximation, we find that the bound state energy has the scaling,

$\displaystyle\mathcal{E}_{q}$
$\displaystyle=-4\Gamma\sin^{2}(k)\bigg{[}\frac{2J^{2}+\gamma^{2}+3\Gamma^{2}}{\gamma^{2}+3\Gamma^{2}}\bigg{]}\approx-4\Gamma\bigg{[}\frac{2J^{2}+\gamma^{2}+3\Gamma^{2}}{\gamma^{2}+3\Gamma^{2}}\bigg{]}k^{2}$
$\displaystyle\mathcal{D}$
$\displaystyle=4\Gamma\bigg{[}\frac{2J^{2}+\gamma^{2}+3\Gamma^{2}}{\gamma^{2}+3\Gamma^{2}}\bigg{]}$
(S85)

Note, the above result preserves the feature that when $\Gamma=0$,
$\mathcal{E}_{q}=0$ because there is not a bound state with purely real
energy, but rather an entire spectrum of equally spaced bound states with
purely complex energies from the emergent Stark ladder. If instead we take the
limit, $\gamma=0$; we recover the bound state energy for static hopping with
bond noise. We point out that in this case, when momentum increases, the
diffusive bound state does not meet a continuum of scattering states but
rather other bound states from the localization. Interestingly, these bound
states do not appear in the correlation function, indicating that the ground
state we found is significantly more dominant.

##  S9 Diffusion Constant for Linear Potential with Onsite Noise.

In this section, we study a similar model as above except with the noise on
each site by coupling to the density rather than the hopping, i.e.,
$\Gamma=0$. Moreover, the boundary and bulk equations have the form,

$\displaystyle\mathcal{E}_{q}A_{0}$
$\displaystyle=t_{k}\left(A_{1}-A_{-1}\right),\quad\mathcal{E}_{q}A_{1}=t_{k}\left(A_{2}-A_{0}\right)+\big{[}i\gamma-\mathcal{V}\big{]}A_{1}$
$\displaystyle\mathcal{E}_{q}A_{-1}$
$\displaystyle=t_{k}\left(A_{0}-A_{-2}\right)-\big{[}i\gamma+\mathcal{V}\big{]}A_{-1},\quad\mathcal{E}_{q}A_{\ell}=t_{k}\left(A_{\ell+1}-A_{\ell-1}\right)+\big{[}i\gamma\ell-\mathcal{V}\big{]}A_{\ell}.$
(S86)

To elucidate the bound state energy scaling, we follow the same method in the
previous section where we first split the above equations into real and
imaginary parts and use the relations,
$A^{\mathbb{R}}_{1}/A^{\mathbb{R}}_{-1}=-1$, $A^{I}_{1}/A^{I}_{-1}=1$, and
$A^{I}_{0}=0$. Next, we use the imaginary equations to determine the ratio
$A^{I}_{1}/A^{\mathbb{R}}_{-1}$ and $A^{I}_{-1}/A^{\mathbb{R}}_{-1}$, which we
find to be,

$\frac{A^{I}_{-1}}{A^{\mathbb{R}}_{-1}}=\frac{-\gamma}{\mathcal{E}_{q}+\mathcal{V}}.$
(S87)

Inserting this relation into the equations for
$A^{\mathbb{R}}_{1}$,$A^{\mathbb{R}}_{-1}$ we find that,

$\frac{A^{\mathbb{R}}_{0}}{A^{\mathbb{R}}_{-1}}=\frac{\gamma^{2}+\mathcal{V}^{2}}{\mathcal{V}t_{k}}$
(S88)

where we have ignored terms with $\mathcal{E}_{q}$. With this approximation,
we find the bound state energy has the scaling,

$\displaystyle\mathcal{E}_{q}$
$\displaystyle=\bigg{[}\frac{-8J^{2}\mathcal{V}}{\gamma^{2}+\mathcal{V}^{2}}\bigg{]}k^{2}$
$\displaystyle\mathcal{D}$
$\displaystyle=\bigg{[}\frac{8J^{2}\mathcal{V}}{\gamma^{2}+\mathcal{V}^{2}}\bigg{]}$
(S89)

Note, the above result preserves the feature that when $\mathcal{V}=0$,
$\mathcal{E}_{q}=0$ because there is not a bound state with a purely real
energy, but rather an entire spectrum of bound states with purely complex
energies that are equally spaced from the emergent Stark ladder. Moreover, if
$\gamma=0$ then we recover the result,
$\mathcal{E}_{q}=\frac{-8J^{2}}{\mathcal{V}}k^{2}$ discussed in the main text
of static hopping with onsite noise.

##  S10 Diffusion Constant for General Operator Equation of Motion.

In this section, we present a derivation for the diffusion constant for Eq.
(S51) with all all terms contributing which is one of the primary results in
the main-text. Moreover, for this case the model effectively maps to a one-
dimensional non-hermitian hopping model in an imaginary linear potential with
delta function potential, the eigenvalue equation becomes,

$\displaystyle\mathcal{E}_{q}A_{0}$
$\displaystyle=t_{k}\left(A_{1}-A_{-1}\right)-4\Gamma\sin^{2}(k)A_{0},\quad\mathcal{E}_{q}A_{1}=t_{k}\left(A_{2}-A_{0}\right)-2\Gamma
A_{1}+\Gamma A_{-1}+i\gamma A_{1}-\mathcal{V}A_{1}$
$\displaystyle\mathcal{E}_{q}A_{-1}$
$\displaystyle=t_{k}\left(A_{0}-A_{-2}\right)-2\Gamma A_{-1}+\Gamma
A_{1}-i\gamma
A_{-1}-\mathcal{V}A_{-1},\quad\mathcal{E}_{q}A_{\ell}=t_{k}\left(A_{\ell+1}-A_{\ell-1}\right)-2\Gamma)A_{\ell}+i\gamma\ell
A_{\ell}-\mathcal{V}A_{\ell}.$ (S90)

Just like in the previous two sections, we split the above eigenvalue
equations into their respective real and imaginary parts and utilize the
relations $A^{\mathbb{R}}_{1}/A^{\mathbb{R}}_{-1}=-1$,
$A^{I}_{1}/A^{I}_{-1}=1$, and $A^{I}_{0}=0$ to simplify the equations. Using
the imaginary equations for $A_{\pm}$, we determine the following ratio,

$\frac{A^{I}_{-1}}{A^{\mathbb{R}}_{-1}}=\frac{-\gamma}{\mathcal{E}_{q}+\mathcal{V}+\Gamma}.$
(S91)

This result then leads to ratio,

$\frac{A^{\mathbb{R}}_{0}}{A^{\mathbb{R}}_{-1}}=\frac{\gamma^{2}+(3\Gamma+V)(\Gamma+V)}{2J\sin(k)(\Gamma+V)}$
(S92)

where we have ignored terms with $\mathcal{E}_{q}$. With this approximation,
we find the bound state energy has the scaling,

$\displaystyle\mathcal{E}_{q}$
$\displaystyle=-8\bigg{[}\frac{\Gamma}{2}+\frac{J^{2}(\mathcal{V}+\Gamma)}{\gamma^{2}+(3\Gamma+V)(\Gamma+V)}\bigg{]}k^{2}$
$\displaystyle\mathcal{D}$
$\displaystyle=8\bigg{[}\frac{\Gamma}{2}+\frac{J^{2}(\mathcal{V}+\Gamma)}{\gamma^{2}+(3\Gamma+V)(\Gamma+V)}\bigg{]}$
(S93)

##  S11 Numerical Methods for Determining the Coefficients.

We now outline how in the main text we numerically determined the two-point
correlation function from the coefficient, $A_{\ell,k}$. First, we put an
upper limit on the operator length labelled $\ell_{\text{max}}$ such that,

$\vec{A}=\begin{pmatrix}A_{-\ell_{\text{max}},k},\ldots,&A_{0,k},\ldots,&A_{\ell_{\text{max}},k}\end{pmatrix}.$
(S94)

The above vector obeys the matrix differential equation,

$\frac{d\vec{A}(t)}{dt}=M\vec{A}(t)$ (S95)

where the matrix $M$ is determined from Eq. (10) in the main text which is
diagonalized as $M=VDV^{-1}$ where $V$ is set of eigenvectors and
$D=\text{diag}(\lambda_{1},\ldots,\lambda_{2\ell_{\text{max}}+1})$. Upon
making the above change of basis the differential equation above becomes,

$\frac{d\vec{\Sigma}(t)}{dt}=D\vec{\Sigma}(t).$ (S96)

Here $\vec{\Sigma}=V^{-1}\vec{A}$ and the equation is easily integrated to
give, $\vec{\Sigma}(t)=e^{Dt}\vec{\Sigma}(0)$. Mapping back to the original
basis, the coefficients are found to be,

$\vec{A}(t)=\left(Ve^{Dt}V^{-1}\right)\vec{A}(0).$ (S97)

We then extract the element $[\vec{A}(t)]_{0,k}$ for each momentum and then
perform the following sum,

$C_{x,y}(t)=\frac{1}{8\pi}\sum_{k}[\vec{A}(t)]_{0,k}e^{i(x-y)k}.$ (S98)