ch2_1Drandompaper.tex

\def\maxnt{\mathrm{Max}^{N}_t}
\def\envnt{\mathrm{Env}^{N}_t}
\def\snt{\mathrm{Sam}^{N}_t}

% \makeatletter
% \makeatother

\chapter{Anomalous Fluctuations of Extremes in Many-Particle Diffusion}
\label{ch2_1Drandom}
%\author{Jacob, me, etc}
%\affiliation{$^*$Department of Physics and Materials Science Institute, University of Oregon, Eugene, Oregon 97403, USA \\ $^\dagger$Department of Physics and Astronomy University of Pennsylvania, Philadelphia, Pennsylvania 19104, USA}

From J.B. Hass, A.N. Carroll-Godfrey, I. Corwin, E.I. Corwin, ``Anomalous fluctuations of extremes in many-particle diffusion'' , Vol. 107, Issue 2, pp. L022101. DOI: 10.1103/PhysRevE.107.L022101 (2023).

\section{Abstract}
In many-particle diffusions, particles that move the furthest and fastest can play an outsized role in physical phenomena. A theoretical understanding of the behavior of such extreme particles is nascent. A classical model, in the spirit of Einstein's  treatment of single-particle diffusion, has each particle taking independent homogeneous random walks. This, however, neglects the fact that all particles diffuse in a common and often inhomogeneous environment that can affect their motion. A more sophisticated model treats this common environment as a space-time random biasing field which influences each particle's independent motion. While the bulk (or typical particle) behavior of these two models has been found to match to high degree, recent theoretical work of Barraquand, Corwin and Le Doussal on a one-dimensional exactly solvable version of this random environment model suggests that the extreme behavior is quite different between the two models. We transform these asymptotic (in system size and time) results into physically applicable predictions. Using high precision numerical simulations we reconcile different asymptotic phases in a manner that matches numerics down to realistic system sizes, amenable to experimental confirmation. We characterize the behavior of extreme diffusion in the random environment model by the presence of a new phase with anomalous fluctuations related to the Kardar-Parisi-Zhang universality class and equation.

\section{Introduction}
Our world is fueled by outliers. Information in signals is carried by the leading edge \cite{saxton_single-particle_1997,pinto_physics_2000, hofling_anomalous_2013, ghosh_non-universal_2015, manzo_review_2015,iyer-biswas_first_2015,metzler_non-brownian_2016,  shen_single_2017}.
A viral or bacterial infection is spread by the first few pathogens to enter a host and the first host to enter a new region \cite{hufnagel_forecast_2004,kao_virus_2006}.
A species is evolved by the fittest mutations \cite{ewens_mathematical_2012, metzler_first-passage_2014,waxman_diffusion_2017}. Scientific revolution is sparked by the first new idea. In all of these contexts the precipitating action is driven by the extremes among a great number of agents (varying from $N\sim 10^2$ to $N\sim 10^{60}$ depending on the context) evolving in a complex but shared environment. How does the nature of the shared environment affect these outlier behaviors? Conversely, can we infer the nature of the shared environment from the behavior of these outliers?  Despite their obvious importance, these overarching questions are still unanswered.

The classical model for many-particle diffusion as independent homogeneous random walks provides an easily calculable solution, but entirely neglects the effects of the shared and likely inhomogeneous environment. This model is the basis for diffusion coefficients~\cite{einstein_uber_1905,einstein_zur_1906, einstein_theoretische_1907}, which succinctly describe the behavior of typical particles in a many-particle diffusion. A more sophisticated model treats the shared environment as a space-time random biasing field with short-range space-time correlations. Each particle thus articulates independent random walks subject to forcing by the common biasing field. While this refined model does not affect typical particle diffusion behavior \cite{rassoul-agha_quenched_2013}, it drastically impacts the behavior of extreme particles. In this work, we provide predictions for the behavior of extreme particles moving in a random and inhomogenous environment.
We find that the variance in the position of the extreme particle is a robust and sensitive measurement of the nature of the environment and show how this variance can be understood as the sum of two contributions: the randomness present in the environment, and the sampling of random walks in that environment.  We show that by subtracting out the variance due to sampling we can produce direct measurements of the environment, inaccessible from measurements of the motion of a typical particle or of the bulk. This residual environmental variance is characterized by a novel power law that we demonstrate holds even when the number of particles is as small as a few hundred.

\begin{figure}[h]
\begin{center}
\includegraphics[width=.9\columnwidth]{Figures/Figure1.pdf}
\caption{A system of $N=10^5$ particles evolving in a given random environment. The heat map records the site occupancy density. We also plot in green the asymptotic theory mean location for the maximum particle location. Around this is a shaded region with a width of two standard deviations based on the asymptotic theory variance. This region generally contains the extreme-most particle over time. The zoomed-in inset shows the spatial locations of $N=10^2$ particles over time. Color indicates the bias (red is biased down and blue is biased up) and is chosen independently at each space-time box. The location of particles within each box is chosen for ease of visualization.}
  \label{fig:BCModel}
\end{center}
\end{figure}

\section{Background}
Building on observations by Brown~\cite{brown_xxvii_1828, brown_xxiv_1829} from 1827,  Einstein~\cite{einstein_uber_1905,einstein_zur_1906, einstein_theoretische_1907} (along with Langevin~\cite{langevin_sur_1908}, Sutherland~\cite{sutherland_lii_1893,sutherland_lxxv_1905} and Smoluchowski~\cite{von_smoluchowski_zur_1906, von_smoluchowski_notiz_1915}) proposed a theory of diffusion based on modeling particles by independent random walks with variance controlled by a diffusion coefficient intrinsic to the particle/environment pair. Soon after, Perrin experimentally verified Einstein's statistical predictions \cite{perrin_mouvement_1909, perrin_mouvement_1910}.

Probing the effectiveness and limitations of Einstein's diffusion model has remained a challenge. On short time scales, particle motion is ballistic, dominated by inertia~\cite{uhlenbeck_theory_1930,huang_direct_2011,hammond_direct_2017,lukic_direct_2005, franosch_resonances_2011, kheifets_observation_2014}.
Many physically relevant situations require the addition of new concepts to accurately model them.
Certain diffusive processes are better modeled by Levy flights \cite{wang_when_2012} or other types of anomalous diffusions~\cite{bouchaud_anomalous_1990, metzler_brownian_2019} instead of simple random walks. Other work has focused on active particles which inject energy into their environment~\cite{ramaswamy_mechanics_2010,kanazawa_loopy_2020}. Further, in environments which are slowly mixing, Einstein's theory may also break down due to the presence of quenched disorder \cite{zangi_frequency-dependent_2007, wang_when_2012}. Unlike the above deviations from the classical model, our approach is intended to describe generic many-particle diffusions.

The random walk in random environment (RWRE) model goes back to~\cite{chernov_replication_1967, temkin_one-dimensional_1972} (see also~\cite{havlin_diffusion_2002, bolthausen_ten_2002,sznitman_topics_2004, ofer_random_2004, hughes_random_1995}) and comes in two types -- long-range~\cite{kesten_limit_1975,sinai_limiting_1983,bouchaud_classical_1990, bouchaud_anomalous_1990, burlatsky_transient_1998, le_doussal_random_1999} and short-range~\cite{richardson_atmospheric_1926, hentschel_relative_1984, bouchaud_diffusion_1990,  chertkov_anomalous_1996, bernard_anomalous_1996, jullien_richardson_1999, balkovsky_intermittent_2001} temporally correlated environments. We focus here on the latter. In this context, typical RWRE particles behave like Brownian motion, matching the behavior from Einstein's model~\cite{rassoul-agha_almost_2005, deuschel_quenched_2016}.
The motion of atypical particles is controlled by large deviations of the RWRE's transition probability as first studied in ~\cite{balazs_random_2006}.

Barraquand and Corwin~\cite{barraquand_random-walk_2017} discovered the exactly solvable Beta RWRE discussed extensively below and uncovered a remarkable connection between its large deviations for times of order $\log(N)$ and the statistics of the Kardar-Parisi-Zhang (KPZ) universality class~\cite{corwin_kardar_2012, quastel_one-dimensional_2015}, namely the Gaussian Unitary Ensemble (GUE) Tracy-Widom distribution~\cite{tracy_level-spacing_1993}. Soon after~\cite{le_doussal_diffusion_2017} recognized that a phase transition should occur in the $(\log(N))^2$ time frame while~\cite{barraquand_moderate_2020} discovered that in this frame the GUE Tracy-Widom distribution is replaced by the KPZ equation one-point distribution~\cite{kardar_dynamic_1986,sasamoto_one-dimensional_2010,calabrese_free-energy_2010,dotsenko_bethe_2010,amir_probability_2011}. See~\cite{sabot_random_2017,balazs_large_2019, barraquand_large_2020, brockington_bethe_2021, oviedo_second_2021,korotkikh_hidden_2022,krajenbrink_crossover_2022} for further developments. The recursion relation \eqref{eq:kolmogorov} for RWRE transition probabilities solves a discrete version of the multiplicative noise stochastic heat equation (mSHE)
\begin{equation}
\partial_t Z(x,t) = \frac{1}{2} \partial_x^2 Z(x,t) + \xi(x,t)Z(x,t)
\end{equation}
with $\xi$ space-time white noise. The logarithm of the mSHE $h(x,t)=\log Z(x,t)$ solves the KPZ equation
\begin{equation}\label{eq:KPZ}
\partial_t h(x,t) = \frac{1}{2} \partial_x^2 h(x,t) + \frac{1}{2}(\partial_x h(x,t))^2 + \xi(x,t).
\end{equation}
Hence, large deviations for RWREs, in particular beyond the solvable model and even in experimental settings, may relate to the KPZ equation and its universality class -- especially in light of the rich canon of work on KPZ universality in various contexts using theoretical~\cite{alberts_intermediate_2010, corwin_kardar_2012}, numerical~\cite{krug_amplitude_1992, prolhac_height_2011}, and experimental~\cite{halpin-healy_kpz_2015} methods. The KPZ connection is  quite useful since its statistics and power-laws are well studied.

\section{Models for diffusion}
Although physical diffusion is continuous in time and (typically) occurs in three-dimensional space, here we work with discrete models in one spatial dimension. The principal reason for this choice is that it is the setting for the exactly solvable  Beta RWRE~\cite{barraquand_random-walk_2017} (a continuous {\it sticky Brownian motion} limit of this model exists~\cite{barraquand_large_2020}) that will enable us to compare numerical results to exact theoretical predictions.
Beyond that, discretization is common for numerical simulations and higher dimensions are more challenging numerically due to anisotropy issues arising from the choice of lattice and due to the lack of exactly solvable models, c.f. \cite{le_doussal_diffusion_2017}. In real diffusion in a common environment, there will be length and time scales on which the environment decorrelates.  Our discrete model can be thought of as coarse-graining the environment in space and time onto a lattice and thus we do not expect discrete and continuous models to differ greatly for long-times and large-scales. Our model ignores any higher order interactions as we expect them to be less present in the behavior of extreme particles, for which the local density is necessarily low. Additionally, there are physical settings where particles take discrete states~\cite{feller_diffusion_1951, moran_random_1958} or evolve in quasi-one-dimensional spaces~\cite{pollard_gaseous_1948,ahmadi_diffusion_2017}.

We study the {\it Beta RWRE} introduced in~\cite{barraquand_random-walk_2017} (see Fig.~\ref{fig:BCModel}). We model the environment by a collection $\mathbf{B}= \big\{B(x,t):x\in \mathbb{Z},t\in \mathbb{Z}_{\geq 0}\big\}$ of independent identically distributed random variables all drawn from the uniform distribution on $[0,1]$. At time $t=0$ we start with $N$ particles all at site $0$. Given an instance of the environment $\mathbf{B}$ the particles proceed as follows. Each particle at $x$ and $t$ independently flips the same weighted coin which has probability $B(x,t)$ of heads (moving the particle to site $x+1$ at time $t+1$) and $1-B(x,t)$ of tails (moving to $x-1$ instead). Thus, while particles do not interact with each other, those at the same place and time are all influenced by the common environment.

This model is exactly solvable when $B(x,t)$ are distributed according to the Beta distribution, $Beta(\alpha, \beta)$~\cite{barraquand_random-walk_2017}. For simplicity, we focus on the special case $\alpha=\beta=1$ corresponding to the uniform distribution. The classical simple symmetric random walk (SSRW) model arises in the limit $\alpha=\beta\to \infty$ where all $B(x,t)\equiv 1/2$ and the environment is deterministic.

We focus on the behavior of the right-most particle at time $t$. We denote this by $\maxnt$, with $N$ the number of particles in the system. Two types of randomness affect $\maxnt$: that of the environment and that of sampling the random walks in that environment. The effect of the environment is via the transition probability $p_{\mathbf{B}}(x,t)$, the probability that a single random walker initially at $0$ will end up at $x$ at time $t$ for a given environment $\mathbf{B}$. This satisfies the recursion relationship, %(i.e., the Kolmogorov backward equation, or master equation),
%
\begin{align} \label{eq:kolmogorov}
 \begin{split}
  p_{\mathbf{B}}(x,t) = & p_{\mathbf{B}}(x-1,t-1)B(x-1,t-1) +\\
  & p_{\mathbf{B}}(x+1,t-1)\big(1-B(x+1,t-1)\big)
 \end{split}
\end{align}
%
with initial condition $p_{\mathbf{B}}(0,0) = 1$ and $p_{\mathbf{B}}(x \neq 0,0) = 0$. Since each random walker is independent, conditional on the environment, the distribution of the ensemble of $N$ walks is determined by $p_{\mathbf{B}}(x,t)$. Given the environment $\mathbf{B}$, the probability that a single random walker is at or above $x$ at time $t$ is given by the tail probability, $P_{\mathbf{B}}(x,t) = \sum_{y\geq x} p_{\mathbf{B}}(y,t)$. This and the independence of random walkers, conditional on the environment, imply
%
\begin{equation}\label{eq:maxpb}
\mathrm{Prob}_{\mathbf{B}}(\maxnt \leq x) = \big(1-P_{\mathbf{B}}(x,t)\big)^N,
\end{equation}
%
where the left-hand side is the probability, given the environment $\mathbf{B}$, that $\maxnt \leq x$.

We study how $\maxnt$ varies upon sampling a new environment and random walkers therein. Eq. \eqref{eq:maxpb} suggests that a good proxy for $\maxnt$ is the location $\envnt$ of the $1/N$-quantile of $P_{\mathbf{B}}(x,t)$, i.e., $\envnt$ equals the maximal $x$ such that $P_{\mathbf{B}}(x,t)>1/N$. Notice that $\envnt$ only accounts for the variation due to the environment. The variation due to sampling in that environment is denoted $\snt$ and defined by $\maxnt = \envnt + \snt$. We use the notation $\mean{\bullet}$ and $\var{\bullet}$ for the mean and variance of a quantity $\bullet$ (e.g. $\maxnt, \envnt,\snt$) averaged over both the environment and the sampling of random walkers in that environment.

\section{Numerical Methods}
We numerically simulate our models for system sizes varying from $N=10^2$ to $N=10^{300}$. We consider such large and physically unrealistic system sizes like $10^{300}$ in order to see how asymptotic theory applies for as wide a  range as possible of finite system sizes. We evolve the system for times from $t=0$ to $t= 5000 \log(N)$.  As explained below, $\log(N)$ and $(\log(N))^2$ set key timescales and our range of times ensure that for all choices of $N$, we encompass these scales. We simulate such large systems by tracking occupation variables instead of individual particle trajectories. In particular, if there are $N(x,t)$ particles at site $x$ at time $t$, then the number that move to site $x+1$ is binomially distributed with $N(x,t)$ samples and success probability $B(x,t)$ (the remainder move to site $x-1$). We sample these binomial distributions utilizing quadruple-precision floating point numbers and making approximations to the binomial distribution when dealing with sizes beyond our precision limits, as described in Section \ref{ch2_supp}. The rightmost particle location (identified by the maximal $x$ with $N(x,t)\geq 1$) at each time represents a sample of $\maxnt$. By repeatedly sampling new environments along with random walk occupation variables $N(x,t)$ therein we numerically measure $\var{\maxnt}$. To distinguish from the true value we denote this numerically measured variance by $\varnum{\maxnt}$ and plot it in Fig. \ref{fig:MaxVar}. In like fashion, we measure $\varnum{\envnt}$ for each sampled environment by using Eq. \ref{eq:kolmogorov} to compute $p_{\mathbf{B}}(x,t)$.  Fig. \ref{fig:QuantileVar} shows $\varnum{\envnt}$ as a function of time (see  Section \ref{ch2_supp} for $\meannum{\maxnt}$ and $\meannum{\envnt}$). The data presented in Fig. \ref{fig:MaxVar} and \ref{fig:QuantileVar} took approximately three weeks to run in parallel on 500 cores of the University of Oregon high performance computing cluster, Talapas.

\begin{figure}[h]
\begin{center}
 \includegraphics[width=0.9\columnwidth]{Figures/MaxVar.pdf}
 \caption{Plots of $\varnum{\maxnt}$ (solid) computed over 10000, 5000, 1000, 500 and 500 environments (respectively) and  $\varasy{\maxnt}$ (dashed) for $N=10^2, 10^{7}, 10^{24}, 10^{85}, 10^{300}$}
 \label{fig:MaxVar}
\end{center}
\end{figure}

\begin{figure}[h]
\begin{center}
	\includegraphics[width=0.9\columnwidth]{Figures/QuantileVar.pdf}
	\caption{Plots of $\varnum{\envnt}$ (transparent solid) computed over 500 environments, $\varasy{\envnt}$ (dashed), and $\varnum{\maxnt} - \varasy{\snt}$ (dark solid) smoothed in each $1/25^{th}$ of a decade for $N=10^2, 10^{7}, 10^{24}, 10^{85},10^{300}$. The three curves agree as shown in the zoomed-in inset.}
	\label{fig:QuantileVar}
\end{center}
\end{figure}

\section{Asymptotic Theory Results}
We describe asymptotic results on the behavior of $\maxnt$, $\envnt$ and $\snt$ as both $N$ and $t$ tend to infinity in different limits. Given a fixed relationship between $t$ and $\log(N)$ such as $t/\log(N)=\hat{t}$ or  $t/\log(N)^2=\hathat{t}$ for $\hat{t}$ or $\hathat{t}$ fixed, we write $f(N,t)\gg g(N,t)$ if $f(N,t)/g(N,t)$ tends to infinity as $N$ and $t$ do subject to their relationship. We use the notation $\varasy{\bullet}$ to denote the asymptotic theory formula for the variance of $\bullet$, interpolated back to finite $N$ and $t$. SSRW theory follows from Stirling's formula while asymptotic results for the RWRE rely on tools from quantum integrable systems~\cite{barraquand_random-walk_2017, barraquand_moderate_2020,krajenbrink_crossover_2022} and are derived first for $\envnt$ and then for $\maxnt$ and $\snt$.

\noindent\textbf{SSRW $\boldsymbol{\maxnt}$:}
For $t/\log(N)=\hat{t}$ with fixed $\hat{t} < (\log 2)^{-1}$, we have $N\gg 2^t$ and hence with very high probability every reachable site in the lattice at time $t$ is occupied, hence $\var{\maxnt}\approx 0$. When $\hat{t} > (\log 2)^{-1}$, we show  in  Section \ref{ch2_supp} that $\maxnt$ is asymptotically a Gumbel random variable. For $\hat{t}$ large, $\varasy{\maxnt} \approx \frac{\pi^2}{12} \frac{t}{\log(N)}$.

%the location of the maximal random walk at time $t$ will likely occur at a position significantly below $t$ with non-trivial fluctuations which can be calculated by combining \eqref{eq:maxpb} with Stirling's formula for binomial coefficient asymptotics. Doing this in  Section \ref{ch2_supp} we see that $\maxnt$ is asymptotically distributed as a Gumbel random variable. In particular, for $t\gg \log(N)$, $\varasy{\maxnt} \approx \frac{\pi^2}{12} \frac{t}{\log(N)}$.

\noindent\textbf{RWRE $\boldsymbol{\envnt}$:}
For $t/\log(N)=\hat{t}$ with fixed $\hat{t}<1$, $\var{\envnt}\approx 0$. To see this, note that $P_\mathbf{B}(t,t) = B_{0,0} \cdots B_{t-1,t-1}$. Taking logs and applying the central limit theorem shows that $\log \big(P_\mathbf{B}(t,t)\big) \approx -t  + t^{1/2} G$ for $G$ a standard Gaussian. This implies that $P_\mathbf{B}(t,t) \approx e^{-t}\gg 1/N$. Thus the RWRE stops saturating the lattice when $t= \log(N)$ plus a order $(\log(N))^{1/2}$ Gaussian fluctuation. For the SSRW this happens at time $\log_2(N)$ plus order one fluctuations.


$\var{\envnt}$ displays two asymptotic regimes. For fixed $t/\log(N)=\hat{t}>1$, $\var{\envnt}$ takes  asymptotic form,
\begin{equation}\label{eqlogN}
V_1(N, t):=\Big(\frac{\log(N)}{t}\Big)^{2/3} \sigma_{\chi}^2 \frac{2^{2/3}\big(1-\frac{\log(N)}{t}\big)^{4/3}}{1- \big(1- \frac{\log (N)}{t}\big)^2},
\end{equation}
where $\sigma_{\chi}^2\approx~0.813$ is the variance of the GUE Tracy-Widom distribution \cite{prahofer_universal_2000, tracy_level-spacing_1993}. As shown in  Section \ref{ch2_supp}, this follows from the result of \cite{barraquand_random-walk_2017}: For $v\in(0,1)$, $\log P_\mathbf{B}(vt,t)=-t I(v) + t^{1/3} \sigma(v)\chi_t$ where $I(v) = 1-\sqrt{1-v^2}$, $\sigma(v) = (2I(v)^2/(1-I(v)))^{1/3}$ and $\chi_t$ is random converging to the GUE Tracy-Widom distribution as $t$ goes to infinity.

For $t/(\log(N))^2\!=\!\hathat{t}$,  $\var{\envnt}$ takes asymptotic form
\begin{equation}\label{eqlogNsq}
V_2(N, t):=\frac{t}{2\log(N)} \cdot \var{h\Big(0,\frac{4 (\log(N))^2}{t}\Big)},
\end{equation}
where $h(0,s)$ is the height at $0$ and time $s$ of the {\it narrow wedge} solution to the KPZ equation \eqref{eq:KPZ}. As shown in  Section \ref{ch2_supp}, this follows from \cite{barraquand_moderate_2020}: For $v\in (0,\infty)$,
$
\log P_\mathbf{B}(vt^{3/4},t) \approx -\frac{v^2t^{1/2}}{2} -\log(t)/4+\log(v) - v^4/12 + h(0,v^4).
$

Interpolating between these regimes, and extrapolating past $(\log(N))^2$ (see also \cite{krajenbrink_crossover_2022}) we find two power-laws,
\begin{equation}\label{eq:varqnt}
\varasy{\envnt} \approx
    \begin{cases}
    \sigma_{\chi}^2 (\frac{\log(N)}{2})^{\frac{1}{3}} t^{\frac{1}{3}}& 1 \ll \frac{t}{\log(N)}\ll \log(N),\\
    \frac{1}{2}\pi^{\frac{1}{2}} t^{\frac{1}{2}}& \frac{t}{\log(N)}\gg \log(N).
    \end{cases}
\end{equation}
For finite $N$ and $t$ these regimes have a gentle crossover that we capture by setting
$
\varasy{\envnt} := I(N,t) V_1(N, t)+\big(1-I(N,t)\big) V_2(N, t)
$
where $I(N,t) := \frac{1}{2} \cdot \left(1-\text{erf}\left(\frac{t-(\log(N))^{3/2}}{(\log(N))^{4/3}}\right)\right)$ (with  $\text{erf}(x)= 2/\sqrt{\pi} \int_{0}^{x} e^{-t^2}dt$ the error function) interpolates from $1$ to $0$ over an interval of width $(\log(N))^{4/3}$ around $(\log(N))^{3/2}$.


\noindent\textbf{RWRE $\boldsymbol{\snt}$ and $\boldsymbol{\maxnt}$:}
We identify the additional contribution from sampling the many-particle diffusion given an environment.
Using Eq. \eqref{eq:maxpb} and Taylor expansion of the results of \cite{barraquand_random-walk_2017} and \cite{barraquand_moderate_2020} quoted above,  Section \ref{ch2_supp} shows that for $t/\log(N)=\hat{t}>1$ the sample fluctuation $\snt$ is of Gumbel type with variance
%
\begin{equation}\label{eq:varnst}
\varasy{\snt} = \frac{\pi^2}{6} \frac{\big(\frac{t}{\log(N)} -1\big)^2}{2\frac{t}{\log(N)} -1}\approx \frac{\pi^2}{12} \frac{t}{\log(N)}
\end{equation}
as $\hat{t}$ grows. This limit matches the behavior of the SSRW model. In Section \ref{ch2_supp} we also show that $\snt$ is asymptotically independent of $\envnt$, thus
\begin{equation}\label{eq:varadd}
\var{\maxnt} \approx \var{\envnt} + \var{\snt}.
\end{equation}

\begin{figure}[h]
\begin{center}
\includegraphics[width=0.9\columnwidth]{Figures/MaxQuantComp.pdf}
\caption{Variance of the maximal particle $\var{\maxnt}$ (red), environment $\var{\envnt}$ (blue), and sampling $\var{\snt}$ (green) for RWRE, and variance of the maximal particle  $\var{\maxnt}$ (purple) for SSRW, all for $N=10^7$. Dashed lines are $\varasy{\bullet}$, while solid lines are $\varnum{\bullet}$. $\varnum{\maxnt} - \varasy{\snt}$ (orange; smoothed as in Fig. \ref{fig:QuantileVar}) closely matches the environment curve (blue).}
\label{fig:RWREAssembly}
\end{center}
\end{figure}

\section{Comparison of Numerical and Theoretical Results}
Fig. \ref{fig:MaxVar} and \ref{fig:QuantileVar} show that the asymptotic theoretical predictions for $\var{\maxnt}$ and $\var{\envnt}$ are in excellent agreement with the numerical measurements.  Fig. \ref{fig:QuantileVar} further shows that we reliably recover $\var{\envnt}$ using $\varnum{\maxnt} - \varasy{\snt}$, as expected from Eq. \ref{eq:varadd}. Notably, while these results were derived for asymptotically large $\log(N)$ and $t$, they hold nearly perfectly down to $N=10^{2}$.
Fig. \ref{fig:QuantileVar} reveals that while we readily see the long-time $t^{1/2}$ power-law for $\var{\envnt}$ from Eq. \ref{eq:varqnt}, the $t^{1/3}$ power-law is elusive. Although the full characterization of the short-time regime is in excellent agreement with the numerical results, the $t^{1/3}$ power-law is difficult to capture since the transitional window of $\log(N)$ to $(\log(N))^2$ is too narrow for realistic sizes of $N$, even up to $N=10^{300}$. By measuring the long-time $t^{1/2}$ power law, we measure the short-time scaling behavior of the KPZ equation up to a prefactor using Eq. \ref{eqlogNsq}. Fig. \ref{fig:RWREAssembly} shows the tight matching of the asymptotic theory curves and numerically measured values for the variance of $\maxnt, \envnt$ and $\snt$ for a given value of $N=10^7$. Notice that for $t\approx \log(N)$ the asymptotic theory and numerical values for the variance of $\snt$ do not fit as well as for large $t$. This is likely a result of finite-size effects and quickly goes away at larger values of $t$ or when $N$ increases. The fit for $N=10^{300}$ in Fig. \ref{fig:MaxVar} and \ref{fig:QuantileVar} remains tight over the entire range of $t$.

\section{Conclusion}
The link between RWREs and KPZ universality with its wealth of theoretical, numerical and experimental evidence strongly suggests that aspects of the picture presented here will persist beyond discrete and solvable models, even to experiments. When $t$ is of order $\log(N)$, variances should be non-universal, depending in a difficult to determine way on the nature of the environment. By contrast, when $t\gg\log(N)$, we anticipate that the scaling exponents and functional forms we have identified for the variances of $\envnt,\snt$ and $\maxnt$ will be universal, as will the relation \eqref{eq:varadd}. The leading coefficients  in Eq. \eqref{eq:varqnt} should be non-universal and hold within them all of the accessible information about the correlation structure of the environment -- we call these {\it extreme diffusion coefficients}. Further theoretical study, such as for the general $\alpha,\beta$ Beta RWRE model, should provide a natural first test of this universal picture and an understanding of how the extreme diffusion coefficients relate to the microscopic environment. A continuum model that should provide an even wider testing-ground amenable to numerics involves particles $x_i(t)$ for $i=1,2,\ldots$ satisfying $dx_i(t) = F(x_i(t),t)dt + D(x_i(t),t) dB_i(t)$ where $F(x,t)$ and $D(x,t)$ are random forcing (as in~\cite{le_doussal_diffusion_2017}) and diffusivity (generalizing diffusing diffusivity, c.f.~\cite{chechkin_brownian_2017}) fields common to all particles while $B_i$ are  Brownian motions independent between different $i$. Changing the correlation structures of $F$ and $D$ will probe the transition between temporally mixing versus quenched environments, which should have very different behavior (c.f.~\cite{noskowicz_average_1988,le_doussal_first-passage_1989}) and warrants further study. Considering higher dimensions as in~\cite{le_doussal_diffusion_2017} may lead to further theory that better model real physical systems. A study of higher order cumulants may reveal other ways to probe the hidden environment, although they may be harder to observe numerically or experimentally.

In physical systems it is impossible to directly measure the environmental variance. However, an indirect measurement can be performed via the approach presented here by using $\var \envnt\approx \var \maxnt - \var \snt$. The sample variance $\var \snt$ is now computed using
$
 \var \snt = \frac{\pi^2 D}{6} \frac{t}{\log(N)}
$
where $D$ is the diffusion coefficient. One could repeatedly track the motion of the leading edge of diffusing particles in a system of colloids confined to a quasi-1D channel thereby directly measuring $\var \maxnt$ for system sizes ranging from $N ~ 10^2$ to $N ~ 10^{10}$. Further, one can also perform complementary measurements on the time of first passage of diffusing objects, which opens the door to experiments done on all manner of diffusing objects, including light or sound diffusing through a scattering medium, dye molecules in a fluid, or any other object whose first passage can be measured.
By measuring the environmental variance and extreme diffusion coefficient we will gain a new microscope through which to probe the hidden nature of the underlying environment in which the diffusion occurs. Our work should serve as a guide in the development and analysis of novel experimental measurements of the extreme behavior of many-particle diffusion. Two examples of such experimental work will be described in the following chapters.

\section{Acknowledgements}
We thank G. Barraquand and P. Le Doussal for discussions and S. Prolhac for providing numerics for $\var{h(0,s)}$.
This work was funded under the W.M. Keck Foundation Science and Engineering grant on ``Extreme Diffusion''.
I.C. also wishes to acknowledge ongoing support from the NSF through DMS:1811143 and DMS:1937254, the Simon Foundation through a Simons Fellowship in Mathematics (Grant No. 817655), and the Packard Foundation Fellowship for Science and Engineering. Much of this work was performed while I.C. held a Miller Visiting Professorship from the Miller Institute for Basic Research in Science, and while in residence at the Mathematical Sciences Research Institute in Berkeley, California (NSF Grant No. 1440140).
E.I.C. wishes to acknowledge ongoing support from the Simons Foundation for the collaboration Cracking the Glass Problem via award 454939. This work benefited from access to the University of Oregon high performance computing cluster, Talapas.