# Towards Brain-inspired System: Deep Recurrent Reinforcement Learning for
Simulated Self-driving Agent

Jieneng Chen1, Jingye Chen2, Ruiming Zhang1 and Xiaobin Hu3 *This work was
financially supported by the German Research Foundation (DFG) and the
Technical University of Munich (TUM) in the framework of the Open Access
Publishing Program. This research was also funded by the Chinese Ministry of
Education’s National University Student Innovation and Entrepreneurship
Training Program (2018).1Jieneng Chen and Ruiming Zhang are with College of
Electronics and Information Engineering, Tongji University, 4800 Cao An
Highway, Shanghai, China $\\{$chenjn<EMAIL_ADDRESS>Chen
is with Department of Computer Science, Fudan University, Shanghai, China
<EMAIL_ADDRESS>Hu is with Department of Computer Science,
Technical University of Munich, Munich, Germany<EMAIL_ADDRESS>
###### Abstract

An effective way to achieve intelligence is to simulate various intelligent
behaviors in the human brain. In recent years, bio-inspired learning methods
have emerged, and they are different from the classical mathematical
programming principle. In the perspective of brain inspiration, reinforcement
learning has gained additional interest in solving decision-making tasks as
increasing neuroscientific research demonstrates that significant links exist
between reinforcement learning and specific neural substrates. Because of the
tremendous research that focuses on human brains and reinforcement learning,
scientists have investigated how robots can autonomously tackle complex tasks
in the form of a self-driving agent control in a human-like way. In this
study, we propose an end-to-end architecture using novel deep-Q-network
architecture in conjunction with a recurrence to resolve the problem in the
field of simulated self-driving. The main contribution of this study is that
we trained the driving agent using a brain-inspired trial-and-error technique,
which was in line with the real world situation. Besides, there are three
innovations in the proposed learning network: raw screen outputs are the only
information which the driving agent can rely on, a weighted layer that
enhances the differences of the lengthy episode, and a modified replay
mechanism that overcomes the problem of sparsity and accelerates learning. The
proposed network was trained and tested under a third-partied OpenAI Gym
environment. After training for several episodes, the resulting driving agent
performed advanced behaviors in the given scene. We hope that in the future,
the proposed brain-inspired learning system would inspire practicable self-
driving control solutions.

##  I Introduction

Recently, research in brain science has gradually received the public’s
attention. Given the rapid progress in brain imaging technologies and in
molecular and cell biology, much progress has been made in understanding the
brain at the macroscopic and microscopic levels. Currently, the human brain is
the only truly general intelligent system that can cope with different
cognitive functions with extremely low energy consumption. Learning from
information processing mechanisms of the brain is clearly the key to building
stronger and efficient machine intelligence [1]. In recent years, some bio-
inspired intelligent methods have emerged, and they are definitely different
from the classical mathematical programming principle. Bio-inspired
intelligence has the advantages of strong robustness and efficient, good
distributed computing mechanism. Besides, it is easy to combine with other
methods.

The mammalian brain has multiple learning subsystems. Niv et al. [2]
categorized major learning components into four classes: the neocortex, the
hippocampal formation (explicit memory storage system), the cerebellum
(adaptive control system) and the basal ganglia (reinforcement learning ).
Among these learning components, reinforcement learning is particularly
attractive to researches. Nowadays, converging evidence links reinforcement
learning to specific neural substrates, thus assigning them to precise
computational roles. Most notably, much evidence suggests that the
neuromodulator known as dopamine provides basal ganglia target structures with
phasic signals that convey a reward prediction error which can influence
learning and action selection, particularly in stimulus-driven habitual
instrumental behaviors [3]. Hence, many efforts have been made to investigate
the capability of bio-inspired reinforcement learning by applying them to
artificial intelligence related tasks. [4, 5, 6, 7, 8, 9]

In recent years, deep reinforcement learning has contributed to many of the
spectacular success stories of artificial intelligence in recent years [10,
11]. After the initial success of the deep Q network (DQN) [12], a variety of
improved models have been published successively. Later on and based on the
former discoveries, Mnih et al. [7] proposed the Nature DQN in 2015 and
introduced the replay memory mechanism to break the strong correlations
between the samples. Mnih et al. [13] proposed a deep reinforcement learning
approach, in which the parameters of the deep network are updated by multiple
asynchronous copies of the agent in the environment. Van Hasselt et al. [14]
suggested the Double DQN to eliminate overestimation; they added a target Q
network independent from the the current Q network. It was shown to apply to
large-scale function approximation [14]. Newer techniques included deep
deterministic policy gradients and mapping an observation directly to action,
both of which could operate over continuous action spaces [15]. Schaul et al.
[16] suggested prioritized replay, adding prior to the replay memory to
relieve the sparse reward and slow converge the problem slowly [16]. In the
case of partially observable states, the recurrent neural network (RNN) and
long short-term memory (LSTM) have been proved to be effective in processing
sequence data [17]. hausknecht2015deep [18] replaced the last fully connected
layer in the network with an LSTM layer. Also, Foerster et al. [19] proposed
to use multi-agent to describe a state distributively. However, most previous
studies such as Atari games focused on the simple environment and action
space. Because of this limitation, there is an urgent need to further improve
the capability of deep reinforcement learning in a more challenging and
complex scenario such as simulated 3D driving control problem.

One deep-learning-based studied the simulated self-driving game [20]. However,
three problems existed in their implementation. First, they created a
handcrafted dataset. Obviously, one can never create this ideal benchmark
dataset that includes all the bad situations encountered by the driving agent
during training. At best, one can include the best behavior that the driving
agent should implement in each step. The driving agent was reported to perform
well when it had a good position in the driveway. However, the behavior
deteriorated rapidly when the driving agent deviated from the driveway. Such
behaviors indicated the dissimilar distribution and instability even though
correctional measures were taken on the dataset. Second, they trained and
supervised their network in a supervised way. As there are many possible
scenarios, manually tackling all possible cases using supervised learning
methods will likely yield a more simplistic policy [21]. Third, their
experiments were built on ideal conditions; for example, they assumed that the
brakes were ignored. In our experiments, we take the brakes into
consideration. Moreover, to support autonomous capabilities, a robotic driven
agent should adopt human driving negotiation skills when overtaking, halting,
braking, taking left and right turns, and pushing ahead in unstructured
roadways. It comes naturally that a trial-and-error way of study is more
suitable for this simulated self-driving game. Hence, the bio-inspired
reinforcement learning method in the study is a more suitable way for the
driving agent to learn how to make decisions.

In our study, we proposed a deep recurrent reinforcement learning network to
solve simulated self-driving problems. Rather than creating a handcrafted
dataset and training in a supervised way, we adopted a bio-inspired trail-and-
error technique for the driving agent to learn how to make decisions.
Furthermore, this paper provides three innovations. First, intermediate game
parameters were completely abandoned, and the driving agent relied on only raw
screen outputs. Second, a weighting layer was introduced in the network
architecture to strengthen the intermediate effect. Third, a simple but
effective experience recall mechanism was applied to deal with the sparse and
lengthy episode.

The rest of this study is organized as follows: Section II describes deep
Q-learning, recurrent reinforcement learning, network architecture and
implementation details. Section III verifies experimental results. The
conclusion of this study is drawn in Section IV .

##  II Methodology

Deep Q-learning is used to help AI agents operate in environments with
discrete actions spaces. Based on the knowledge of Deep Q-learning, we
proposed a modified DRQN model in order to infer the full state in partially
observable environments.

###  II-A Deep Q-learning

Reinforcement learning manages learning policies for an agent interacting in
an unknown environment. In each step, an agent observes the current states of
the environment, makes decisions according to a policy $\pi$, and observes a
reward signal $r_{t}$ [22]. Given the current states and a set of available
actions, the main aim of the DQN is to approximate the maximum sum of
discounted rewards. The value of a given action-state pair or the ‘Q-value’ is
given by:

$Q^{*}\left(s_{t},a_{t}\right)=maxR_{t+1}$ (1)

According to the Bellman equation, it gives the approximating form of Q-values
by combining the reward obtained with the current state-action pair and the
highest Q-value at the next state $S^{{}^{\prime}}$, and the best action
$a^{{}^{\prime}}$:

$Q(s,a)=Q(s,a)+\gamma*\max_{a^{\prime}}Q\left(s^{\prime},a^{\prime}\right)$
(2)

We often use the form involving a iterative process:

$Q(s,a)=Q(s,a)+\alpha\left(r+\gamma*\max_{a^{\prime}}Q\left(s^{\prime},a^{\prime}\right)-Q(s,a)\right)$
(3)

Actions are typically chosen following a $\varepsilon$-greedy exploration
policy. The epsilon value $\varepsilon$ was valued between 0.0 to 1.0. In
order to encourage the agent to explore the environment, the $\varepsilon$ was
set to 1.0 at the beginning. During the training process, the value decayed
gradually as the experience accumulated. Then, the agent could use the
experience to achieve the task.

$P(a|s)=\frac{\varepsilon}{n(A)}+1-\varepsilon*\operatorname{argmax}_{a}Q(s,a)$
(4)

The network stands for the approximation of the expected Q-value $Q^{*}$,
which leads to the loss function:

$\operatorname{Loss}(\theta)=\sum\left(y_{t}-Q_{\theta_{t}}(s,a)\right)^{2}$
(5)

###  II-B Recurrent Reinforcement Learning

Figure 1: The modified DRQN model. The value function was divided into two
categories: the current value function $Q$ and target value function
$Q^{{}^{\prime}}$. The parameters in $Q$ were assigned to $Q^{{}^{\prime}}$
per N episodes. The state contained two elements: $o_{t}$ gained from current
environment and $h_{t-1}$ gained from former information. The agent performed
action $a$ using a specific policy, and the sequence
$\left(o_{t},a,r,o_{t}^{{}^{\prime}}\right)$ was stored in the replay memory
unit. We used a prioritized experience replay memory unit here. During
training, the sequence was randomly chosen from the replay memory unit. We
trained the network using gradient descent to make the current value function
$Q$ approach $Q^{{}^{\prime}}$ given a specific sequence. The loss function
was shown in Eq. 5.

For some special games which are three-dimensional and partially observable,
the DQN lacks the ability to solve the problem. In partially observable
environments, the agent only receives an observation $o_{t}$ of the current
environment, which is usually insufficient to infer the full state of the
system. The real state $s_{t}$ is the combination of the current observation
$o_{t}$ and an unfixed length of history states. Hence, we adopted the DRQN
model on top of the DQN to deal with such conditions (see Figure 1). The last
full connect layer was replaced by the LSTM in the DRQN model in order to
record former information. Figure 2 shows the sequential updates in the
recurrent network. An additional input $h_{t-1}$ stood for the previous
information is added to the recurrent model. The output of the LSTM
$z\left(o_{t},h_{t-1}\right)$ , which combined the current observation $o_{t}$
and the history information $h_{t-1}$, was used to approximate the Q-value
$Q\left(o_{t},h_{t-1},a_{t}\right)$. The history information was updated and
passed through the hidden state to the network in the next time step:

$h_{t}=LSTM\left(h_{t-1},o_{t}\right)$ (6)

Figure 2: Sequence updates in the recurrent network. Only the scores of the
actions taken in states 5, 6 and 7 will be updated. The first four states
provide a more accurate hidden state to the LSTM, while the last state
provides a target for state 7.

###  II-C Network Architecture

Figure 3: An illustration of the architecture of our model. The input image is
assigned to four convolution layers. The output of the convolution layers is
split into two streams. The first stream (bottom) flattens the output and
feeds it to an LSTM. The second one (top) flattens the output and feeds it to
a full connect layer. Then, we obtained an advantage stream and value steam
individually and multiplied them as the output. We stored the information in
prioritized experience memory unit. As was shown in Figure 1, the network was
trained using the DRQN loss function.

In the beginning, we used the baseline DRQN model carried out by the CMU in
VizDoom2016 [22]. However, we obtained unsatisfied results with the same
model. Hence, we have made three improvements in our modified model to make
things work better. The whole architecture is shown in Figure 3. First, the
network was built on top of the NVIDIAs autopilot model [23]. To reduce
overfitting, the original model was modified by adding several batch
normalization layers. We used a four-layer stronger CNN for feature
extraction. The input size was set to 320*240. The first convolutional layer
contained 32 kernels, with the a size of 8*8, and the stride was 4. The second
convolutional layer contained 64 kernels, with the a size of 4*4 and the
stride was 4. The third layer contained 128 kernels, with the a size of 3*3
and the stride was 1. The last convolutional layer contained 256 kernels with
a size of 3*3 and the stride is 1. The activate functions were all Relu, and
the size of pooling layers were all 2*2. Second, we abandoned the fully
connected layer before the LSTM layer in the original DRQN model and fed the
LSTM directly with the high-level feature. The number of units in LSTM was set
to 900. Third, the subsequent structure was divided into two groups for
different purposes. One was mapped to the set of possible actions, and the
other was a set of scalar values. The final action value function was acquired
using both of them. We will introduce their functions respectively.

We used two collateral layers rather than a single DRQN network to approximate
the value function. The auxiliary layer also had the same dimension as the
action space. The original intention was to balance the impact of history and
the current. The DRQN model enabled the driving agent to make a reasonable
action from a fully observing perspective. Nevertheless, we wanted to focus
our attention on the precise instantaneous change of the current view. The
auxiliary layer was mapped to [0, 1.0] using the softmax method, thus
suggesting the direction of correction for the raw approximations using DRQN.
Because of this intervention, the network would not only learn the best action
for a state but also understand the importance of taking actions. Imagine an
agent driving in a straight line and having an obstacle far ahead, an original
DRQN model could learn that it is time to make the driving agent move a bit to
avoid hitting the obstacle. The modified model will also have an insight into
when it is the best time to move, with the knowledge that the danger of the
obstacle increases as it gets closer. $V(s,a)$ was used to represent the
original output of the DRQN, and $A(s,a)$ was used to represent the
significance provided by the auxiliary layer. We used the formula to express
the final strengthen of the Q-value (see Figure 3):

$Q(s,a)=V(s,a)^{T}*A(s,a)$ (7)

The result was stored in the prioritized experience memory unit. During
training, the sequence in the memory unit was removed, and we used the Eq. 5
to calculate the loss.

###  II-D Implementation Details

Reinforcement learning consists of two basic concepts: action and reward.
Action is what an agent can do in each state. Given that the screen is the
input, a robot can take steps within a certain distance. An agent can take
finite actions. When a robot takes an action in a state, it receives a reward.
Here, the term reward is an abstract concept that describes the feedback from
the environment. A reward can be positive or negative. When the reward is
positive, it corresponds to our normal meaning of reward. However, when the
reward is negative, it corresponds to what we usually call punishment. We also
describe the training details such as the hyperparameters, input size
selection, frameskip, and prioritized replay.

####  II-D1 Action Space

The game had five actions, including Left, Right, Straight, Brake and
Backwards. The rocking range of the joystick reflects the numerical value of
the speed control, which is mapped into a region of -80 to 80. To keep the
model simple, we considered that only the turning control had a high
requirement for precision. The speed section was discretized into a set of [0,
20, 40, 80] for speed control. Another three actions were represented by 1/0
flag. Thus, each action was represented by a 5-dimensional vector shown below:

$\displaystyle\text{$actions$}$ $\displaystyle=[-40,0,1,0,0],\text{ $left$ }$
$\displaystyle=[40,0,1,0,0],\text{ $right$ }$
$\displaystyle=[0,0,0,1,0],\text{ $go$ $backwards$ }$
$\displaystyle=[0,0,1,0,0],\text{ $go$ $straight$ }$
$\displaystyle=[0,0,0,0,1],\text{ $brake$ }$

####  II-D2 Reward

Under most circumstances, the driving agent cannot explore a path with big
rewards initially. The driving agent often gets stuck somewhere in the halfway
and wait for the time to elapse before resetting. We have to make the rewards
of these cases variant in order to make these experience meaningful. Hence, we
set a series of landmarks along the track. A periodical reward was given to
the driving agent when each landmark was reached. The closer the distance
between the landmark and the destination, the bigger the phased reward was
given. We have considered formulating a more precise reward system to increase
density, such as giving the driving agent a slight punishment when it deviates
from the driveway. However, hidden game parameters were needed, and they were
against the original aim of this research.

####  II-D3 Hyperparameters

The network was trained using the RMSProp algorithm and the minibatch size was
40. The replay memory was set to contain 10000 recent frames. The learning
rate $\alpha$ was set, from 1.0 to 0.1, to follow a linear degradation and
then fixed at 0.1. The discount factor $\gamma$ was set to 0.9. The
exploration rate $\varepsilon$ for the $\varepsilon$-greedy policy was set to
follow a linear degradation, from 1.0 to 0.1, and then fixed at 0.1. The
exploration rate $\varepsilon$ was set to 0 when we evaluated the model. The
original screen outputs were three channel RGB images. They were first
transformed into grey scale images and then fed to the network.

####  II-D4 Input size

In the beginning, the original 640*480 screen resolution was resized to
160*120 to accelerate training. After several hours of trials, the driving
agent still got stuck in most cases and could not complete one lap. The
resulting rewards oscillated around the baseline reward for not finishing the
game in the limited number of steps, thus indicating the resolution was too
low to recognize. Then, the input size was resized to 320*240. We also reduced
the punishments to encourage positive rewards. After the observation of the
same length of time, the distribution of the resulting rewards cast off random
oscillation and started to turn positive over the baseline. Hence, we kept
this resolution as the input size of the network despite the memory
consumption as the system started to learn within the acceptable limit of
time.

####  II-D5 Frameskip

Obviously, it is not necessary to feed every single frame into the network
because there is only a weak change between adjacent frames. Because of this,
we used the frame-skip technique [24] as in most approaches. One frame was
taken out as the network input by every k + 1 frames. The same action was
repeated over the skipped frames. The higher k was set, the faster the
training speed. However, the control effect became imprecise at the same time.
Considering the balance of low computing resource consumption and smooth
control, we finally set a frame skip of k = 3 by relying on our experience.

####  II-D6 Prioritized replay

In most cases, there was little replay memory with high rewards, which would
be time-consuming with huge replay table and only sparse rewards. As
prioritized replay was a method that can make learning from experience replay
more efficient, we simply timed the important experiences in proportion to
their rewards and stored them into the replay memory. That way, they would
have a greater opportunity to be recalled [16].

##  III Experiments

The model was trained using three different tracks, which cover all the track
that Stanford used for comparison. An individual sets of weights was trained
separately for each model because each track has different terrain textures.
The rewards were low initially because it is equivalent to a random
exploration at the beginning of training and because the driving agent would
get stuck somewhere without making any significant progress. After about 1400
episodes of training, the driving agent finished the race. Under most
circumstances, the driving agent did not finish the race in given steps so the
reward was positive but not as high as receiving the final big reward. We set
this step limit because of a lack of a reset mechanism for dead situations
which was very useful in the early stage of training. In Stanford’s report,
they created a DNF flag to represent the driving agent getting stuck. In our
experiment, the driving agent had learned better policy and displayed better
behaviors, proving better robustness of the system. We also visualized the CNN
in order to validate the ability of the model.

###  III-A Experimental Environment

We chose the car racing game to carry out our simulated self-driving
experiment. In order to play the car racing game autonomously, we used a
third-partied OpenAI Gym environment wrapper for the game developed by
Bzier111https://github.com/bzier/gym-mupen64plus. The API provides direct
access to the game engine and allows us to run our scripts while playing the
game frame-by-frame. Through the API, we can easily obtain the game
information, either the screen output or the high-level intermediate game
parameters like its location in the minimap. Our model was verified to be able
to deal with the partially observable game environment efficiently.

Figure 4: There are reward trends for two different maps. The picture (a)
corresponded to Farm and the picture (b) corresponded to Raceway. After
training for 1750 episodes, we obtained the reward tendency. At about 400
episodes, the stability of the driving agent began to increase. After 1400
episodes, the reward stabilized at a high level.

###  III-B Rewards Analysis

For comparison, the model was trained and tested using the same tracks like
those used by Stanford in their supervised learning method. Each track has
different terrain textures and difficulty routes. Therefore, an individual set
of weights was trained separately for each model. The experiments were carried
out on a common configured portable laptop, and all models converged after
spending over 80 hours each.

###  III-C CNN Visualization

The quality of the high-level feature was regarded as the key measurement of
the learning process. The network could be evaluated from many aspects such as
the loss function and the validation accuracy in traditional supervised
learning. However, in unsupervised learning, we did not have this kind of
standard to provide a quantitative assessment of the network. Hence, we
visualized the output of the last layer (see Figure 5), with the aim of
ascertaining the kind of features that are captured by the network.The high-
level layers output seemed quite abstract through direct observation. Thus, we
visualized the high-level features through deconvolution. Deconvolution and
unpooling are the common methods for reconstructing entity.

Figure 5: We visualized the last layer and obtained a direct view (left).
Because the high-level layers output was quite abstract, we then visualized
the last layer using deconvolution and obtained the right picture. It seemed
to represent the wall element in the original graph.

###  III-D Result and Discussion

The test results are shown in Table I. The driving agent was tested on several
tracks in line with those of Stanford’s experiment such as Farm, Raceway, and
Mountain. We evaluated our model based on its achieved time. For each track,
we run 10 races in real-time and calculate the mean race times as the final
result. The Stanford results [20] and the Human results were borrowed from the
Stanford report. The human results were obtained by a real human participant
playing each track twice: The first time is to get used to the track, and the
second time is to record their time. The results were obtained by testing the
driving agent on each track ten times and then taking the average (excluding
the DNF cases).

TABLE I: Performance comparison.  Track |  Our model |  Stanford model [20] |  Human
---|---|---|---
Farm |  98.33 |  97.46 |  94.07
Raceway |  166 |  129.09, 1DNF |  125.30
Mountain |  213 |  138.37, 2DNF |  129.50

In sum, our results seemed to require a longer period to finish the race.
Through observation, we found that our agent’s driving path was not as smooth
as that provided by supervised learning results. We closed this course for two
reasons. First, we reduced the epsilon too fast. As shown in the Figure 4, the
epsilon value is rapidly adjusted to 0.1 as the rewards increase. This
indicates the experiences show more reliability. Then during the latter phase
of training, we depended mainly on experiences even though there were still
many better state-action sets to explore. Second, we discretized the actions
roughly. We used a set of [0; 20; 40; 80] as the options for speed. To be more
precise, we used the control variable of the joystick. We observed that the
numerical value of the joystick parameter did not maintain a linear
relationship with the real effect on the driving agent as shown on the screen.
Through observation, a speed of 20 m/s and 40 m/s both produced a tiny effect
while a speed of 80m/s would make a radical change. An accurate modeling
between the joystick parameter and the real effect would make the segmentation
more reasonable. However, compared with the experiment done by the Stanford
group, our experiments performed satisfactory even if the driving agent
deviated from the driveway. We considered the brake and trained the driving
agent using a trial-and-error method, which was more in line with the real
situation. Hence, the bio-inspired reinforcement learning method in the study
was a more suitable approach for the driving agent to make decisions.

It is worth mentioning that the driving agent is more sensitive to the green
grass than the yellow sand texture. The reason is that we used a fixed set of
weights to transform the color image into a gray-scale image before inputting,
which is in charge of this phenomenon. In the future, we will use the original
three channels of the RGB data to carry out the experiments to end this
conjecture and improve the results.

##  IV Conclusion

Brain-inspired learning has recently gained additional interests in solving
decision-making tasks. In this paper, we propose an effective brain-inspired
end-to-end learning method with the aim of controlling the simulated self-
driving agent. Our modified DRQN model has proven to manage plenty of error
states effectively, thus indicating that our trial-and-error method using deep
recurrent reinforcement learning could achieve better performance and
stability. By using the screen pixels as the only input of the system, our
method highly resembles the experience of human-beings solving a navigation
task from the first-person perspective. This resemblance makes this research
inspirational for real-world robotics applications. Hopefully, the proposed
brain-inspired learning system will inspire real-world self-driving control
solutions.

## References

  * [1] M.-m. Poo, J.-l. Du, N. Y. Ip, Z.-Q. Xiong, B. Xu, and T. Tan, “China brain project: basic neuroscience, brain diseases, and brain-inspired computing,” _Neuron_ , vol. 92, no. 3, pp. 591–596, 2016.
  * [2] Y. Niv, “Reinforcement learning in the brain,” _Journal of Mathematical Psychology_ , vol. 53, no. 3, pp. 139–154, 2009.
  * [3] F. Rivest, Y. Bengio, and J. Kalaska, “Brain inspired reinforcement learning,” in _Advances in neural information processing systems_ , 2005, pp. 1129–1136.
  * [4] M. McPartland and M. Gallagher, “Learning to be a bot: Reinforcement learning in shooter games.” in _AIIDE_ , 2008.
  * [5] J. Peters and S. Schaal, “Reinforcement learning of motor skills with policy gradients,” _Neural networks_ , vol. 21, no. 4, pp. 682–697, 2008.
  * [6] H. Duan, S. Shao, B. Su, and L. Zhang, “New development thoughts on the bio-inspired intelligence based control for unmanned combat aerial vehicle,” _Science China Technological Sciences_ , vol. 53, no. 8, pp. 2025–2031, 2010\.
  * [7] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski _et al._ , “Human-level control through deep reinforcement learning,” _Nature_ , vol. 518, no. 7540, p. 529, 2015.
  * [8] Y. Zhu, R. Mottaghi, E. Kolve, J. J. Lim, A. Gupta, L. Fei-Fei, and A. Farhadi, “Target-driven visual navigation in indoor scenes using deep reinforcement learning,” in _Robotics and Automation (ICRA), 2017 IEEE International Conference on_. IEEE, 2017, pp. 3357–3364.
  * [9] S. Gu, E. Holly, T. Lillicrap, and S. Levine, “Deep reinforcement learning for robotic manipulation with asynchronous off-policy updates,” in _Robotics and Automation (ICRA), 2017 IEEE International Conference on_. IEEE, 2017, pp. 3389–3396.
  * [10] J. Kober, J. A. Bagnell, and J. Peters, “Reinforcement learning in robotics: A survey,” _The International Journal of Robotics Research_ , vol. 32, no. 11, pp. 1238–1274, 2013.
  * [11] P. Henderson, R. Islam, P. Bachman, J. Pineau, D. Precup, and D. Meger, “Deep reinforcement learning that matters,” in _Thirty-Second AAAI Conference on Artificial Intelligence_ , 2018.
  * [12] V. Mnih, K. Kavukcuoglu, D. Silver, A. Graves, I. Antonoglou, D. Wierstra, and M. Riedmiller, “Playing atari with deep reinforcement learning,” _arXiv preprint arXiv:1312.5602_ , 2013.
  * [13] V. Mnih, A. P. Badia, M. Mirza, A. Graves, T. Lillicrap, T. Harley, D. Silver, and K. Kavukcuoglu, “Asynchronous methods for deep reinforcement learning,” in _International conference on machine learning_ , 2016, pp. 1928–1937.
  * [14] H. Van Hasselt, A. Guez, and D. Silver, “Deep reinforcement learning with double q-learning.” in _AAAI_ , vol. 2. Phoenix, AZ, 2016, p. 5.
  * [15] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” in _International Conference on Learning Representations (ICLR)_ , 2016.
  * [16] T. Schaul, J. Quan, I. Antonoglou, and D. Silver, “Prioritized experience replay,” in _International Conference on Learning Representations (ICLR)_ , 2016.
  * [17] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [18] M. Hausknecht and P. Stone, “Deep recurrent q-learning for partially observable mdps,” in _2015 AAAI Fall Symposium Series_ , 2015.
  * [19] J. N. Foerster, Y. M. Assael, N. de Freitas, and S. Whiteson, “Learning to communicate to solve riddles with deep distributed recurrent q-networks,” in _Advances in Neural Information Processing Systems_ , 2016, pp. 2137–2145.
  * [20] H. Ho, V. Ramesh, and E. T. Montano, “Neuralkart: A real-time mario kart 64 ai,” 2017.
  * [21] S. Shalev-Shwartz, S. Shammah, and A. Shashua, “Safe, multi-agent, reinforcement learning for autonomous driving,” _arXiv preprint arXiv:1610.03295_ , 2016.
  * [22] G. Lample and D. S. Chaplot, “Playing fps games with deep reinforcement learning.” in _AAAI_ , 2017, pp. 2140–2146.
  * [23] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in _Advances in neural information processing systems_ , 2012, pp. 1097–1105.
  * [24] M. G. Bellemare, Y. Naddaf, J. Veness, and M. Bowling, “The arcade learning environment: An evaluation platform for general agents,” _Journal of Artificial Intelligence Research_ , vol. 47, pp. 253–279, 2013.


# Equivalence of dynamics of disordered quantum ensembles and semi-infinite
lattices

Hallmann Óskar Gestsson<EMAIL_ADDRESS>Charlie Nation
<EMAIL_ADDRESS>Alexandra Olaya-Castro<EMAIL_ADDRESS>Department of
Physics and Astronomy, University College London, London WC1E 6BT, United
Kingdom

###### Abstract

We develop a formalism for mapping the exact dynamics of an ensemble of
disordered quantum systems onto the dynamics of a single particle propagating
along a semi-infinite lattice, with parameters determined by the probability
distribution of disorder realizations of the original heterogeneous quantum
ensemble. This mapping provides a geometric interpretation on the loss of
coherence when averaging over the ensemble and allows computation of the exact
dynamics of the entire disordered ensemble in a single simulation.
Alternatively, by exploiting the reverse map, one can obtain lattice dynamics
by averaging over realisations of disorder. The potential of this equivalence
is showcased with examples of the map in both directions: obtaining dephasing
of a qubit via mapping to a lattice model, and solving a simple lattice model
via taking an average over realizations of disorder of a unit cell.

Introduction.―A disordered quantum ensemble is composed of quantum systems
whose properties vary randomly from one member of the ensemble to another
according to some distribution. Predicting observables of such ensembles then
necessitates the introduction of stochastic elements into the Hamiltonian
description of a quantum system. Examples of disordered ensembles are found
across a wide range of fundamental and applied scenarios such as condensed
matter physics Anderson (1958); Gornyi  _et al._ (2005); Paredes  _et al._
(2005); Basko  _et al._ (2006); Bermudez  _et al._ (2010); Pal and Huse
(2010); Bardarson  _et al._ (2012); Breitweiser  _et al._ (2020), quantum
gravity and black hole physics Maldacena and Stanford (2016); Pollack  _et
al._ (2020); Peng (2021), quantum optics Schlawin  _et al._ (2012); Gilead
_et al._ (2015); Bai  _et al._ (2016); Lib and Bromberg (2022), statistical
physics Popescu  _et al._ (2006); Cotler  _et al._ (2017); Nation and Porras
(2018); Foini and Kurchan (2019); Nation and Porras (2020); Dabelow and
Reimann (2020), quantum chemistry Frost  _et al._ (2010); Gelin  _et al._
(2021) and quantum biology Olaya-Castro  _et al._ (2008); Kolli  _et al._
(2012); Stross  _et al._ (2016); Cohn  _et al._ (2022).

When one considers observable properties of an ensemble of disordered quantum
systems there will be a noticeable dephasing that is incurred, even if every
individual system follows a closed system (unitary) quantum dynamics:
information regarding coherences is lost in the averaging that is inextricably
carried out when an observable is measured. This is a manifestly classical
mechanism that is fundamentally distinct from those inducing incoherent
dynamics for a quantum system interacting with an environment. Owing to the
broad relevance of such phenomena in physical systems, a central issue of
interest thus surrounds the description of the ensemble averaged dynamics of
disordered quantum systems, and how this may be related to physically relevant
processes.

Recent works have established an effective master equation description of the
ensemble averaged state to analyze its dynamical evolution Gneiting  _et al._
(2016); Kropf  _et al._ (2016); Chen  _et al._ (2018); Gneiting (2020),
indicating that disordered quantum systems may be analyzed by applying the
frameworks developed within the theory of open quantum systems Breuer and
Petruccione (2007). A unitary description of disordered quantum ensembles has
been considered for discrete probability distributions Paredes  _et al._
(2005). Inspired by the chain-mapping technique Bulla  _et al._ (2005, 2008)
that has been applied to open quantum systems Prior  _et al._ (2010); Chin
_et al._ (2010), we consider an alternative and versatile approach where the
Hamiltonian of the entire ensemble continuum is transformed using a unitary
change of basis that leverages the properties of orthogonal polynomials Ismail
(2005).

In doing so, we show that the dynamical description of the disordered ensemble
is unitarily equivalent to a single semi-infinite lattice Teschl (2000) with
couplings that are local for physically relevant models of disorder. Our work
then allows gaining new physical insights on the mechanisms on the loss of
coherence while at the same time establishing new conceptual links between two
apparently dissimilar physical contexts: disordered ensembles and lattice
dynamics.

Figure 1: Equivalence of disordered quantum ensembles to lattice dynamics.

This lattice representation of the dynamics gives a geometric intuition to the
mechanism by which information on coherence is lost when an ensemble average
is carried out. Yielding exact dynamics of the entire continuum of
realizations, which may then be recovered by performing a partial trace over
the lattice degree of freedom, bypassing reliance on numerical quadrature
rules and sampling from the space of realizations. In practical terms, the
advantage of the proposed method stems from the simplicity of the Hamiltonian
representation and the accompanying initial wavefunction in the new basis. It
is furthermore straightforward to consider initial states of the ensemble that
are a function of the disorder, e.g. eigenstates of the system.

This work is arranged as follows. First, we introduce the pure state formalism
for disordered quantum ensembles, and the orthogonal polynomial mapping
procedure, showing the general equivalence between the ensemble average
dynamics and semi-infinite lattice models. To exemplify our approach, we then
discuss the case of linear disorder, and practical aspects of the mapping for
numerical application. We then provide examples of the exact equivalence in
realistic systems. Additional details are provided in the supplemental
material (SM).

Disordered Quantum Ensembles.―The ensemble of disordered quantum systems we
consider is composed of every possible realization of the disorder, which we
assume are not interacting with one another, meaning that each system
belonging to the ensemble evolves independent of all other systems in the
ensemble. Unitary dynamics of a specific realization of disorder is determined
by a Hamiltonian of the form

$\hat{H}_{\lambda}=\hat{H}_{0}+\hat{V}_{\lambda},$ (1)

where $\hat{H}_{0}$ and $\hat{V}_{\lambda}$ are components that are
independent and dependent on the disorder, respectively. The multi-index
$\lambda$ consists of $l$ unit-less independent random variables, such that we
have $\lambda=(\lambda_{1},\lambda_{2},\ldots,\lambda_{l})$ and a factorizing
joint probability distribution function
$p(\lambda)=\prod_{i=1}^{l}p^{(i)}(\lambda_{i})$.

The state space of each realization is spanned by a basis
$\\{\ket{n,\lambda}\\}_{n=1}^{N}$ which we use to expand Eq. (1) as

$\hat{H}_{\lambda}=\sum_{n,m=1}^{N}\left(\braket{n}{\hat{H}_{0}}{m}+f_{n,m}(\lambda)\right)\ket{n,\lambda}\bra{m,\lambda},$
(2)

where the $\braket{n}{\hat{H}_{0}}{m}$ elements are constant for all $\lambda$
so it is dropped as a label, whilst the
$f_{n,m}(\lambda)=\braket{n,\lambda}{\hat{V}_{\lambda}}{m,\lambda}$ elements
vary as a given function of the disorder. We furthermore expand the initial
state of each realization with respect to this basis as
$\ket{\psi_{0},\lambda}=\sum_{n=1}^{N}c_{n}(\lambda)\ket{n,\lambda}$, with
$c_{n}(\lambda)=\braket{n,\lambda}{\psi_{0},\lambda}$. The Hamiltonian of the
entire ensemble can be stated in terms of an integral over all possible
realizations as

$\hat{H}_{\text{Ens}}=\int\text{d}\lambda\hat{H}_{\lambda},$ (3)

with an initial state of the ensemble of the form

$\ket{\Psi_{0}}=\int\text{d}\lambda\sqrt{p(\lambda)}\ket{\psi_{0},\lambda}.$
(4)

Note that the distribution of disorder informs the initial state of the
ensemble (different distributions of disorder result in different initial
states) and that $\ket{\Psi_{0}}$ is normalized as
$\braket{\Psi_{0}}{\Psi_{0}}=\int\text{d}\lambda p(\lambda)=1$. This is a
proper physical choice for the initial state of the ensemble as the population
density of realizations is given by $p(\lambda)$, such that the probability
amplitude may be taken as $\sqrt{p(\lambda)}$. With $\hat{H}_{\text{Ens}}$ and
$\ket{\Psi_{0}}$ at hand we are now in position to determine exact dynamics of
the entire ensemble.

Recovering ensemble averaged dynamics.―We now express disorder averaged
operators as the partial trace over the disorder label $\lambda$,

$\bar{O}(t)=\int\text{d}\lambda\braket{\lambda}{O(t)}{\lambda},$ (5)

where $O(t)$ is a generic operator onto the Hilbert space of the entire
ensemble, and $\bar{O}(t)$ is its corresponding disorder average. In
particular we have that for average density matrix dynamics with its initial
state to be of the form given in Eq. (4) we have

$\displaystyle\bar{\rho}(t)$
$\displaystyle=\int\text{d}\lambda\braket{\lambda}{e^{-i\hat{H}_{\text{Ens}}}}{\Psi_{0}}\braket{\Psi_{0}}{e^{i\hat{H}_{\text{Ens}}}}{\lambda}$
(6) $\displaystyle=\int\text{d}\lambda
p(\lambda)\braket{\lambda}{e^{-i\hat{H}_{\lambda}t}}{\psi_{0},\lambda}\braket{\psi_{0},\lambda}{e^{i\hat{H}_{\lambda}t}}{\lambda},$
(7)

where we have made use of the fact that realizations are mutually independent
such that
$\braket{\lambda}{e^{-i\hat{H}_{\text{Ens}}}}{\psi_{0},\lambda^{\prime}}=\braket{\lambda}{e^{-i\hat{H}_{\lambda}t}}{\psi_{0},\lambda}\delta(\lambda-\lambda^{\prime})$.
In the special case where each initial state is independent of $\lambda$ we
can further simplify the expression for $\bar{\rho}(t)$ by making use of
$\braket{\lambda}{e^{-i\hat{H}_{\lambda}t}}{\psi_{0},\lambda}=e^{-i\hat{H}_{\lambda}t}\ket{\psi_{0}}$,
yielding

$\bar{\rho}(t)=\int\text{d}\lambda
p(\lambda)e^{-i\hat{H}_{\lambda}t}\ket{\psi_{0}}\bra{\psi_{0}}e^{i\hat{H}_{\lambda}t},$
(8)

such that we are able to recover from our pure state formalism the average
density matrix dynamics considered in Refs. Gneiting  _et al._ (2016); Kropf
_et al._ (2016); Chen  _et al._ (2018); Gneiting (2020).

Mapping ensemble dynamics to a semi-infinite lattice.―The typical strategy for
calculating average dynamics of a disordered quantum ensemble is to sample the
space of realizations, either randomly or according to numerical quadratures,
which can be parallelized in a straightforward manner as
$\hat{H}_{\text{Ens}}$ is block diagonal in $\lambda$. We consider a different
approach and perform a change of basis using a family of orthonormal
polynomials that transforms the integral representation of Eq. (3) into a
discrete semi-infinite lattice with nearest-neighbour hopping terms. This
procedure is equivalent to a Lanczos tridiagonalization of the ensemble
continuum Hamiltonian when the disorder parameter enters as a linear term into
the Hamiltonian de Vega  _et al._ (2015). Note this approach is also used in
the study of operator complexity growth Parker  _et al._ (2019); Mück and Yang
(2022).

We now introduce the set of polynomials,
$\\{\phi_{k}(\lambda)\\}_{k=0}^{\infty}$, that are mutually orthonormal with
respect to the measure $p(\lambda)\text{d}\lambda$ (see SM for a brief
introduction to orthogonal polynomials, or e.g. Ref. Ismail (2005)) and define
a new discrete basis,

$\ket{n,\lambda}=\sum_{k=0}^{\infty}\sqrt{p(\lambda)}\phi_{k}(\lambda)\ket{n,k}.$
(9)

The disorder has been assumed to consist of $l$ independent random variables
such that we have the factorization
$p(\lambda)=\prod_{i=1}^{l}p^{(i)}(\lambda_{i})$, which results in a
factorization of the polynomials as
$\phi_{K}(\lambda)=\prod_{i=1}^{l}\phi_{k_{i}}(\lambda_{i})$, where we have
introduced a discrete multi-index $K=(k_{1},k_{2},\ldots,k_{l})$. The
disorder-independent component of the ensemble Hamiltonian in the discrete
basis will reduce as

$\int\text{d}\lambda\ket{n,\lambda}\bra{m,\lambda}=\sum_{K}\ket{n,K}\bra{m,K}.$
(10)

As for the disorder-dependent terms we have

$\int\text{d}\lambda
f_{n,m}(\lambda)\ket{n,\lambda}\bra{m,\lambda}=\sum_{K,K^{\prime}}f_{n,m}^{(K,K^{\prime})}\ket{n,K}\bra{m,K^{\prime}},$
(11)

where $f_{n,m}^{(K,K^{\prime})}=\int\text{d}\lambda
p(\lambda)f_{n,m}(\lambda)\phi_{K}(\lambda)\phi_{K^{\prime}}(\lambda)$, such
that the Hamiltonian of the ensemble in the lattice basis is written as

$\hat{H}_{\text{Ens}}=\sum_{K,K^{\prime}}\left(\hat{H}_{0,K}\delta_{K,K^{\prime}}+\sum_{n,m=1}^{N}f_{n,m}^{(K,K^{\prime})}\ket{n,K}\bra{m,K^{\prime}}\right),$
(12)

with
$\hat{H}_{0,K}=\sum_{n,m=1}^{N}\braket{n}{\hat{H}_{0}}{m}\ket{n,K}\bra{m,K}$
and $\delta_{K,K^{\prime}}=\prod_{i=1}^{l}\delta_{k_{i},k^{\prime}_{i}}$. The
SM provides additional details on how one may arrive at Eq. (12). The disorder
independent components have been projected onto each node of the lattice,
whilst disorder dependent terms result in inter- and intra-node coupling
terms.

Therefore the dynamics of the ensemble of Hamiltonians is unitarily equivalent
to a single semi-infinite lattice whose dimension is equal to the number of
disorder parameters $l$. We may therefore bypass numerical quadrature
techniques that yield approximate solutions to the ensemble dynamics of Eq.
(3) and whose accuracy generally relies on the number of points sampled on the
integral interval, and directly compute the exact dynamics of the ensemble via
straightforward simulation of the unitary dynamics corresponding to Eq. (12).
Utility of this change of basis will however rely on both the initial state of
the ensemble and the form of the $f_{n,m}(\lambda)$ functions. Note that a
disordered Hamiltonian with many random variables will generally result in a
high dimensional lattice which can become numerically intractable, such that
randomly sampling the space of realizations would be the practical numerical
approach for computing average dynamics in such cases. We are able to
effectively simulate the exact ensemble dynamics corresponding to low
dimensional lattices by utilizing the sparsity of the Hamiltonian matrix
representation.

The initial state of the ensemble can be expanded in terms of the lattice
basis states as

$\ket{\Psi_{0}}=\sum_{n=1}^{N}\sum_{K}d_{n,K}\ket{n,K},$ (13)

with $d_{n,K}=\int\text{d}\lambda p(\lambda)\phi_{K}(\lambda)c_{n}(\lambda)$
that also serve as expansion coefficients for
$c_{n}(\lambda)=\sum_{K}d_{n,K}\phi_{K}(\lambda)$. A special case of the
initial state of the ensemble is when the wavefunction of each realization of
disorder does not change as a function of the disorder, i.e.
$c_{n}(\lambda)=c_{n}$ for all $\lambda$. We then simply have
$c_{n}(\lambda)=c_{n}\phi_{\underline{0}}(\lambda)$ which results in Eq. (13)
reducing to a single term,

$\ket{\Psi_{0}}=\sum_{n=1}^{N}c_{n}\ket{n,\underline{0}}.$ (14)

The initial state of the ensemble is found to be entirely localized to the
single node at the origin of the lattice due to the distributions of the
disorder and initial states being identical in this special case, which yields
a simple representation of the initial ensemble wavefunction in the lattice
basis. This fact motivates a straightforward truncation scheme to the
numerical simulation of the exact ensemble dynamics by introducing a cutoff to
the lattice at a distance $D$ from the origin such that the ensemble dynamics
up to a time $t<MD$ is accurately simulated, where $M$ is some positive
constant determined by the specifics of the disorder model.

Another form of the initial state of general interest is to consider an
ensemble of eigenstates of $\hat{H}_{\lambda}$, e.g. an ensemble of energy
ground states. As each realization is found in a stationary state this case
reduces to a case of spectral disorder where the ensemble distribution is
determined by the ground state energy distribution $p(E_{G}(\lambda))$, rather
than $p(\lambda)$. The dephasing due to the ensemble average can then be
exactly simulated by performing the change of basis, Eq. (9), with respect to
$p(E_{G}(\lambda))$.

Linear disorder.―A ubiquitous functional form of the disorder component is a
static linear function
$f_{n,m}(\lambda)=\sum_{i=1}^{l}c_{n,m}^{(i)}\lambda_{i}$, where
$c_{n,m}^{(i)}$ are arbitrary coefficients that depend on the specifics of the
model for disorder Anderson (1958); Pal and Huse (2010); Bardarson  _et al._
(2012); Olaya-Castro  _et al._ (2008); Kolli  _et al._ (2012); Gneiting  _et
al._ (2016); Kropf  _et al._ (2016); Chen  _et al._ (2018); Gneiting (2020);
Kropf (2020). For this case we find the ensemble dynamics to map onto a
lattice with nearest-neighbour interaction terms only. The expression for the
expansion coefficients $f_{n,m}^{(K,K^{\prime})}$ will simplify considerably
and be written in terms of polynomial recurrence coefficients as

$f_{n,m}^{(K,K^{\prime})}=\sum_{i=1}^{l}c_{n,m}^{(i)}(\alpha_{k_{i}}\delta_{K,K^{\prime}}+\sqrt{\beta_{k_{i}+1}}\delta_{K,K^{\prime}\pm
1_{i}}),$ (15)

where $1_{i}$ denotes a multi-index which has a single non-zero entry of unity
as its $i$-th component (see SM for derivation of Eq. (15)). The interaction
of the $l$-lattice Hamiltonian is then reduced to nearest-neighbours and we
obtain

$\displaystyle\hat{H}_{\text{Ens}}$
$\displaystyle=\sum_{K}\left(\hat{H}_{0,K}+\sum_{n,m=1}^{N}\sum_{i=1}^{l}c_{n,m}^{(i)}\alpha_{k_{i}}\ket{n,K}\bra{m,K}\right)$
(16)
$\displaystyle+\sum_{n,m=1}^{N}\sum_{i=1}^{l}c_{n,m}^{(i)}\sqrt{\beta_{k_{i}+1}}\ket{n,K}\bra{m,K+1_{i}}+\text{h.c.}$

It is by virtue of the three-term recurrence relation that every set of
mutually orthogonal polynomials satisfy that we may write
$f_{n,m}^{(K,K^{\prime})}$ in terms of their recurrence coefficients (see SM
for more details regarding the recurrence relation).

Figure 2: Dephasing dynamics for an ensemble of disordered qubits for
different distributions of the energy calculated via simulation of the
equivalent lattice model. The form of the distribution has strong influence on
how the ensemble will propagate in time. Here we have set $E=E_{1}-E_{0}$.

Disordered qubit ensemble dephasing.―Let us consider an ensemble of non-
interacting qubits with randomly sampled excited state energies over a
probability distribution $p(\lambda)$, such that each realization is described
by a Hamiltonian of the form

$\hat{H}_{\lambda}=E_{0}\ket{0,\lambda}\bra{0,\lambda}+(E_{1}+\lambda)\ket{1,\lambda}\bra{1,\lambda},$
(17)

where $E_{0}$ and $E_{1}$ are the central energies of the ground and excited
state, respectively, and $\lambda$ is the static disorder. The Hamiltonian of
the ensemble is then mapped onto a nearest-neighbour coupled chain form via
Eq. (16),

$\displaystyle\hat{H}_{\text{Ens}}$
$\displaystyle=\sum_{k=0}^{\infty}\left(E_{0}\ket{0,k}\bra{0,k}+(E_{1}+\alpha_{k})\ket{1,k}\bra{1,k}\right)$
$\displaystyle+\sum_{k=0}^{\infty}\sqrt{\beta_{k+1}}\ket{1,k}\bra{1,k+1}+\text{h.c.}$
(18)

Figure 2 demonstrates the dephasing dynamics for an initial state where each
realization is in the superposition
$\ket{\psi_{0},\lambda}=(\ket{0,\lambda}+\ket{1,\lambda})/\sqrt{2}$ such that
the initial wavefunction of the ensemble may be represented in the lattice
basis as $\ket{\Psi_{0}}=(\ket{0,0}+\ket{1,0})/\sqrt{2}$. Dephasing is shown
for four different probability distributions: Gaussian, Cauchy, Semicircle,
and Uniform. Notice that the semicircle and uniform distributions cases
exhibit coherence revivals, which is a consequence of the fact that they have
finite support Kropf  _et al._ (2016). For this particular case it has been
found that the dephasing dynamics is characterized completely by the
characteristic function of the random variable Kropf  _et al._ (2016) and this
can be seen by noting that each realization will have its coherence to be
proportional to $e^{i\lambda t}$, and that the expected value of $e^{i\lambda
t}$ is the definition of the characteristic function Lukacs (1972).

We can determine the exact error of our numerical approach for this example as
analytical solutions are known. We find that the error is on the order of
machine precision where we have considered the Gaussian, semicircle, and
uniform distributions, confirming that we are recovering exact dynamics. There
is however a finite error for the Cauchy distribution stemming from the energy
cutoff we introduced in order to obtain well defined recurrence coefficients.
See the SM for details on the parameters used and the closed form solutions we
compare to. We also demonstrate in the SM how one may recover the closed form
solution of the ensemble averaged dynamics by tracing over the lattice degree
of freedom.

Figure 3: Disorder averaged dynamics for site populations of a dimer sub-unit
of the light-harvesting complex LH2 Cogdell  _et al._ (2006). Novel non-
unitary dynamics are exhibited in this example whereby relaxation is
incomplete, meaning that the populations reach a new oscillatory state after
an initial relaxation. Parameters used are taken from Novoderezhkin  _et al._
(2003).

Disorder averaged dimer dynamics.―Let us consider an ensemble of non-
interacting dimers whose site energies are sampled from independent identical
Gaussian distributions, such that each realization is described by a
Hamiltonian of the form

$\displaystyle\hat{H}_{\lambda_{1},\lambda_{2}}$
$\displaystyle=(E_{1}+\lambda_{1})\ket{1,\lambda_{1},\lambda_{2}}\bra{1,\lambda_{1},\lambda_{2}}$
(19) $\displaystyle+\hskip
1.13791pt(E_{2}+\lambda_{2})\ket{2,\lambda_{1},\lambda_{2}}\bra{2,\lambda_{1},\lambda_{2}}$
$\displaystyle+V(\ket{1,\lambda_{1},\lambda_{2}}\bra{2,\lambda_{1},\lambda_{2}}+\text{h.c.}),$

where $E_{1}$ and $E_{2}$ are the two central energies of the site states, $V$
is the coupling, and $\lambda_{1}$, $\lambda_{2}$ are the static disorder. The
parameters we assume for this particular dimer model are
$E_{1}=12325\,\text{cm}^{-1}$, $E_{2}=12025\,\text{cm}^{-1}$,
$V=273\,\text{cm}^{-1}$, so as to simulate a dimerized pair of chromophores
Novoderezhkin  _et al._ (2003); Cogdell  _et al._ (2006). The standard
deviation we have chosen is $\sigma=200\,\text{cm}^{-1}$. The change of basis
we propose will now expand the ensemble into a two dimensional semi-infinite
lattice. Performing the trace over the lattice degrees of freedom will then
yield the disorder averaged density matrix in terms of the sites states, whose
populations we have plotted in Fig. 3. We assume that the initial state of the
ensemble is such that every instance is localized to site 1. Every instance of
the system evolves in time according to the Schrödinger equation, with the
population oscillating between sites in a unitary manner. This is not the case
for the disordered ensemble average of the dimer as Fig. 3 demonstrates. Here
we observe an initial relaxation of the average population dynamics, which is
however not a complete relaxation as the average dynamics eventually reach a
new oscillating stationary state.

Wigner semicircle distribution and constant couplings.―In the above text we
have established an equivalency between disordered quantum ensembles and semi-
infinite lattices. When we transform from the $\lambda$ basis and into the $k$
basis, we are effectively leveraging the spectral information we have at hand
($p(\lambda)$) in order to compute the ensemble averaged dynamics. Fully
utilizing the duality that we have established, one can go from a Hamiltonian
of a semi-infinite lattice with couplings that can be interpreted as
recurrence coefficients and (block-) diagonalize the Hamiltonian.

For example, the Wigner semicircle distribution gives rise to recurrence
coefficients that are constant and can therefore be used to generate semi-
infinite lattices that are translationally invariant for a given unit cell.
The corresponding set of orthogonal polynomials will be the Chebyshev
polynomials of the second kind, $U_{k}(\lambda)$, which can be made equivalent
to a sine wave basis expansion with a change of variables Boyd (1989). By
setting $\lambda=\cos\theta$ we have $p(\theta)=\frac{2}{\pi}\sin^{2}\theta$
and $U_{k}(\cos\theta)=\frac{\sin((k+1)\theta)}{\sin\theta}$ such that Eq. (9)
will become
$\ket{n,\theta}=\sqrt{\frac{2}{\pi}}\sum_{k=0}^{\infty}\sin((k+1)\theta)\ket{n,k}$.

In other words, the dynamics of a single particle propagating across a semi-
infinite constant coupling lattice is equivalent to dynamics of an ensemble
whose disorder follows a Wigner semicircle distribution. This stems from the
fact that the eigenstates of the lattice form a continuum that is indexed by a
quantum number which in turn plays the role of disorder in our ensemble
picture. To clarify, what we are highlighting here is not an analytical
technique to investigate lattice dynamics, but rather a new interpretation of
those dynamics.

Conclusions.―We have shown an equivalency between the dynamics of an ensemble
of disordered quantum systems and of a single semi-infinite lattice. This has
been achieved by identifying the pure state description of ensembles which may
then be unitarily transformed to a lattice description using a basis of
orthogonal polynomials. Terms representing disorder in the ensemble picture
will correspond to exchange (hopping) terms in the lattice picture. We
furthermore show that exact ensemble average dynamics may be recovered by
performing a partial trace over the lattice degrees of freedom, yielding non-
unitary evolution that is fundamentally distinct from open quantum systems.

We have demonstrated that the transformation to the lattice picture is exact
by simulating disorder-induced dephasing of a qubit and comparing to known
analytical solutions. We have also studied the average population dynamics for
a disordered dimer which exhibits novel non-unitary dynamics, demonstrating
the generality of our approach. Finally, we have discussed the case of semi-
infinite lattices with constant nearest-neighbour coupling and highlighted
their equivalency to ensembles of quantum systems with disorder that follows a
Wigner semicircle distribution.

Acknowledgements.― We thank the and the Gordon and Betty Moore Foundation
(Grant GBMF8820) for financial support and the Engineering and Physical
Sciences Research Council (EPSRC UK Grant No. EP/V049011/1). We would like to
thank D. Porras for insightful discussions.

Supplemental Material:
Equivalence of dynamics of disordered quantum ensembles and semi-infinite
lattices Hallmann Óskar Gestsson, Charlie Nation, and Alexandra Olaya-Castro

Department of Physics and Astronomy, University College London, London WC1E
6BT, United Kingdom
(Dated: )

##  I Primer on orthogonal polynomials

For a given probability distribution $p$ with finite moments that is defined
for the real domain $\mathcal{D}\subseteq\mathbb{R}$ we may define an inner
product

$\braket{f,g}=\int_{\mathcal{D}}f(x)g(x)p(x)\text{d}x$ (S1)

for arbitrary real-valued functions $f$, $g$ defined on $\mathcal{D}$. There
then exists a sequence of monic orthogonal polynomials
$\\{P_{n}(x)\\}_{n=0}^{\infty}$ that satisfies the following condition,

$\braket{P_{n},P_{m}}=\int_{\mathcal{D}}P_{n}(x)P_{m}(x)p(x)\text{d}x=\zeta_{n}\delta_{n,m},$
(S2)

with $\delta_{n,m}$ being the Kronecker delta, and positive constants
$\zeta_{n}$ where $\zeta_{0}=1$ Ismail (2005). The subscript of the $P_{n}$
elements is taken to correspond to the degree of the polynomial, e.g.
$P_{0}(x)=1$. The sequence forms a basis for real-valued functions on
$\mathcal{D}$ such that for any $f:\mathcal{D}\rightarrow\mathbb{R}$ there
exists a non-zero sequence of real-valued elements
$\\{c_{n}\\}_{n=0}^{\infty}$ such that

$f(x)=\sum_{n=0}^{\infty}c_{n}P_{n}(x),$ (S3)

where $c_{n}=\braket{f,P_{n}}/\zeta_{n}$.

We can analogously introduce a corresponding sequence of orthonormal
polynomials $\\{p_{n}(x)\\}_{n=0}^{\infty}$ whose elements are defined as
$p_{n}(x)=P_{n}(x)/\sqrt{\zeta_{n}}$ and will satisfy

$\braket{p_{n},p_{m}}=\delta_{n,m}.$ (S4)

###  I.1 Recurrence Relation

A powerful property of any sequence of orthogonal polynomials, that is at the
core of the mapping in the main text, is the three-term recurrence relation
they satisfy. For monic polynomials we have

$P_{n+1}(x)=(x-\alpha_{n})P_{n}(x)-\beta_{n}P_{n-1}(x),$ (S5)

with initial conditions: $P_{0}(x)=1$, $P_{1}(x)=x-\alpha_{0}$, and recurrence
coefficients $\alpha_{n}$, $\beta_{n}$, that are determined as

$\alpha_{n}=\braket{xP_{n},P_{n}}/\braket{P_{n},P_{n}},$ (S6)
$\beta_{n}=\braket{P_{n},P_{n}}/\braket{P_{n-1},P_{n-1}}.$ (S7)

We can show that Eq. (S5) holds by expanding $xP_{n}(x)$ in terms of monic
polynomials using Eq. (S3),

$xP_{n}(x)=\sum_{k=0}^{n+1}\frac{\braket{xP_{n},P_{k}}}{\zeta_{k}}P_{k}(x),$
(S8)

where we have noted that the expansion involves only polynomials of degree
$n+1$ or lower, as $xP_{n}(x)$ will be a polynomial of degree $n+1$. Now we
consider the coefficients appearing in the expansion and write

$\displaystyle\braket{xP_{n},P_{k}}$ $\displaystyle=\braket{P_{n},xP_{k}}$
(S9)
$\displaystyle=\sum_{m=0}^{k+1}\frac{\braket{P_{n},xP_{m}}}{\zeta_{m}}\braket{P_{n},P_{m}}$
(S10) $\displaystyle=\sum_{m=0}^{k+1}\braket{P_{n},xP_{m}}\delta_{n,m},$ (S11)

where in going from the first to second line we have expanded $xP_{k}(x)$
again in terms of the monic polynomials. We therefore have that
$\braket{xP_{n},P_{k}}=0$ for $n>k+1$ such that the right hand side of Eq.
(S8) reduces to three terms,

$xP_{n}(x)=d_{n+1}P_{n+1}(x)+d_{n}P_{n}(x)+d_{n-1}P_{n-1}(x).$ (S12)

The leading coefficient of $P_{n}(x)$ is unity by definition such that the
same will be the case for $xP_{n}(x)$, which in turn implies $d_{n+1}=1$. As
for $d_{n-1}$ we can write

$\displaystyle d_{n-1}\zeta_{n-1}$ $\displaystyle=\braket{xP_{n},P_{n-1}}$
(S13) $\displaystyle=\braket{P_{n},xP_{n-1}}$ (S14)
$\displaystyle=\braket{P_{n},P_{n}+d_{n-1}P_{n-1}+d_{n-2}P_{n-2}}$ (S15)
$\displaystyle=\braket{P_{n},P_{n}},$ (S16)

where we have applied Eq. (S12) when going from the second to third line.
Putting all of this together and rearranging the terms of Eq. (S12) we can
write

$P_{n+1}(x)=\left(x-\frac{\braket{xP_{n},P_{n}}}{\zeta_{n}}\right)P_{n}(x)-\frac{\zeta_{n}}{\zeta_{n-1}}P_{n-1}(x),$
(S17)

which is the three-term recurrence relation given in Eq. (S5).

##  II Transforming the integral representation of the ensemble to a lattice
representation

The Hamiltonian of the entire ensemble continuum can be represented in terms
of an integral over all elements of the ensemble as

$\hat{H}_{\text{Ens}}=\int\text{d}\lambda\hat{H}_{\lambda},$ (S18)

where $\lambda$ is a label that spans the elements. Rewriting Eq. (S18) in
terms of the $\\{\ket{n,\lambda}\\}_{n=0}^{N}$ basis yields

$\displaystyle\hat{H}_{\text{Ens}}=$
$\displaystyle\sum_{n,m=0}^{N}\braket{n}{\hat{H}_{0}}{m}\int\text{d}\lambda\ket{n,\lambda}\bra{m,\lambda}+$
(S19) $\displaystyle\sum_{n,m=0}^{N}\int\text{d}\lambda
f_{n,m}(\lambda)\ket{n,\lambda}\bra{m,\lambda}.$

In order to develop the integral terms in Eq. (S19) we make an unitary change
of basis using a sequence of polynomials $\\{p_{k}(\lambda)\\}_{n=0}^{\infty}$
that are mutually orthonormal with respect to a given probability distribution
$p$,

$\ket{n,\lambda}=\sum_{k=0}^{\infty}\sqrt{p(\lambda)}p_{k}(\lambda)\ket{n,k},$
(S20)

whose inverse is

$\ket{n,k}=\int\text{d}\lambda\sqrt{p(\lambda)}p_{k}(\lambda)\ket{n,\lambda}.$
(S21)

We can see that the change of basis is indeed unitary by considering

$\displaystyle\braket{n,k}{n,k^{\prime}}$
$\displaystyle=\iint\text{d}\lambda\text{d}\lambda^{\prime}\sqrt{p(\lambda)p(\lambda^{\prime})}p_{k}(\lambda)p_{k^{\prime}}(\lambda^{\prime})\braket{n,\lambda}{n,\lambda^{\prime}}$
(S22)
$\displaystyle=\iint\text{d}\lambda\text{d}\lambda^{\prime}\sqrt{p(\lambda)p(\lambda^{\prime})}p_{k}(\lambda)p_{k^{\prime}}(\lambda^{\prime})\delta(\lambda-\lambda^{\prime})$
(S23) $\displaystyle=\int\text{d}\lambda
p(\lambda)p_{k}(\lambda)p_{k^{\prime}}(\lambda)$ (S24)
$\displaystyle=\delta_{k,k^{\prime}}.$ (S25)

For terms of the form appearing in the first line of Eq. (S19) we will have

$\displaystyle\int\text{d}\lambda\ket{n,\lambda}\bra{m,\lambda}$
$\displaystyle=\sum_{k,k^{\prime}=0}^{\infty}\int\text{d}\lambda
p(\lambda)p_{k}(\lambda)p_{k^{\prime}}(\lambda)\ket{n,k}\bra{m,k^{\prime}}$
(S26)
$\displaystyle=\sum_{k,k^{\prime}=0}^{\infty}\braket{p_{k},p_{k^{\prime}}}\ket{n,k}\bra{m,k^{\prime}}$
(S27) $\displaystyle=\sum_{k=0}^{\infty}\ket{n,k}\bra{m,k},$ (S28)

where going from the first to second line we have noted that the integral is
the inner product that we had in Eq. (S4) and then from going from the second
to last line we applied the orthogonality of the polynomials. For terms of the
form appearing in the second line of Eq. (S19) we will have

$\displaystyle\int\text{d}\lambda
f_{n,m}(\lambda)\ket{n,\lambda}\bra{m,\lambda}$
$\displaystyle=\sum_{k,k^{\prime}=0}^{\infty}\int\text{d}\lambda
f_{n,m}(\lambda)p(\lambda)p_{k}(\lambda)p_{k^{\prime}}(\lambda)\ket{n,k}\bra{m,k^{\prime}}$
(S29)
$\displaystyle=\sum_{k,k^{\prime}=0}^{\infty}f_{n,m}^{(k,k^{\prime})}\ket{n,k}\bra{m,k^{\prime}}$
(S30)

where going from the first to second line we have made the definition
$f_{n,m}^{(k,k^{\prime})}=\braket{f_{n,m}p_{k},p_{k^{\prime}}}$. We can now
rewrite Eq. (S19) as a semi-infinite lattice model of the form

$\hat{H}_{\text{Ens}}=\sum_{k,k^{\prime}=0}^{\infty}\sum_{n,m=0}^{N}\left(\braket{n}{\hat{H}_{0}}{m}\delta_{k,k^{\prime}}+f_{n,m}^{(k,k^{\prime})}\right)\ket{n,k}\bra{m,k^{\prime}},$
(S31)

where the $\lambda$-independent component of the ensemble has been projected
onto each node of the lattice, whilst the $\lambda$-dependent component is
responsible for the $f_{n,m}^{(k,k^{\prime})}$ lattice coupling terms.

###  II.1 Linear disorder

At this point it is worth noting that we have transformed the integral
representation of $\hat{H}_{\text{Ens}}$ into a lattice model with respect to
an arbitrary probability distribution with well defined moments, and that the
utility of the lattice mapping will depend on whether or not our choice of $p$
will result in a simple form for the $f_{n,m}^{(k,k^{\prime})}$ coupling
constants. We now demonstrate how such a choice is simple to determine for
most models of disorder that are considered in the literature by exploiting
the recurrence relation that orthogonal polynomials satisfy.

Let us assume that $f$ is a linear function of disorder such that
$f(\lambda)=w\lambda$, where $w$ has units of energy and can be considered to
characterize the strength of the disorder contribution. We can then compute
the coupling coefficients of the lattice expansion using the three-term
recurrence relation,

$\displaystyle f^{(k,k^{\prime})}$ $\displaystyle=w\braket{\lambda
p_{k},p_{k^{\prime}}}$ (S32)
$\displaystyle=\frac{w}{\sqrt{\zeta_{k}\zeta_{k^{\prime}}}}\braket{\lambda
P_{k},P_{k^{\prime}}}$ (S33)
$\displaystyle=w\left(\sqrt{\beta_{k+1}}\delta_{k+1,k^{\prime}}+\alpha_{k}\delta_{k,k^{\prime}}+\sqrt{\beta_{k}}\delta_{k-1,k^{\prime}}\right).$
(S34)

We can see that the resulting interaction terms in the lattice will be of the
nearest-neighbour form, and the strength of this coupling is determined by the
square root of the recurrence coefficients, $\sqrt{\beta_{n}}$. We also have
terms corresponding to energy shifts that are given by $\alpha_{n}$. Finally
note that $\alpha_{n}=\alpha$ for all $n$ if the probability distribution is
even with respect to a point $\alpha$, i.e.
$p(\alpha-\lambda)=p(\alpha+\lambda)$. This can be seen by recalling their
definition given in Eq. (S6) and noting that the resulting polynomials will
have an even and odd parity with respect to the point $\alpha$.

###  II.2 Example: Two-level system ensemble

Let us consider an ensemble of non-interacting qubits that have their excited
state energies to be randomly sampled from a probability distribution
$p(\lambda)$, such that the dynamics of each realization is determined by a
Hamiltonian of the form

$\hat{H}_{\lambda}=E_{0}\ket{0,\lambda}\bra{0,\lambda}+(E_{1}+w\lambda)\ket{1,\lambda}\bra{1,\lambda},$
(S35)

where $E_{0}$ and $E_{1}$ are the central energies of the ground and excited
state, respectively, $\lambda$ is drawn from a random distribution
$p(\lambda)$, and $w$ characterizes the width of the static disorder. The
Hamiltonian of the ensemble is then mapped onto a nearest-neighbour coupled
chain form,

$\displaystyle\hat{H}_{\text{Ens}}$
$\displaystyle=\sum_{k=0}^{\infty}\left(E_{0}\ket{0,k}\bra{0,k}+(E_{1}+w\alpha_{k})\ket{1,k}\bra{1,k}\right)$
(S36)
$\displaystyle+w\sum_{k=0}^{\infty}\sqrt{\beta_{k+1}}\ket{1,k}\bra{1,k+1}+\text{h.c.}$

where $\alpha_{k}$ and $\beta_{k}$ are the recurrence coefficients stemming
from performing the lattice expansion with respect to the distribution of the
disorder, $p(\lambda)$. We take each realization to have an initial state of
the form $\ket{\psi_{0}}_{\lambda}=a\ket{0,\lambda}+b\ket{1,\lambda}$, with
$|a|^{2}+|b|^{2}=1$, which will then yield an initial state of the ensemble in
terms of the lattice basis that is of the form
$\ket{\Psi_{0}}=a\ket{0,0}+b\ket{1,0}$. We can now perform a partial trace
over all nodes of the lattice in order to compute the ensemble averaged
dynamics,

$\displaystyle\bar{\rho}(t)$
$\displaystyle=\sum_{k=0}^{\infty}\braket{k}{\rho(t)}{k}$ (S37)
$\displaystyle=\sum_{n,m=0,1}\left[\sum_{k=0}^{\infty}\braket{n,k}{\rho(t)}{m,k}\right]\ket{n}\bra{m},$
(S38)

where $\rho(t)=e^{-i\hat{H}_{\text{Ens}}t}\rho_{0}e^{i\hat{H}_{\text{Ens}}t}$
and
$\rho_{0}=|a|^{2}\ket{0,0}\bra{0,0}+|b|^{2}\ket{1,0}\bra{1,0}+ab^{*}\ket{0,0}\bra{1,0}+a^{*}b\ket{1,0}\bra{0,0}$.

We define functions $c_{n,m}(t)=\sum_{k=0}^{\infty}\braket{n,k}{\rho(t)}{m,k}$
for the sake of brevity and compute the averaged dynamics for this example. We
can immediately determine $c_{0,0}(t)=|a|^{2}$ by noting that $\ket{0,0}$ is a
stationary state of $\hat{H}_{\text{Ens}}$. We can reduce $c_{0,1}(t)$ to a
single term,

$\displaystyle c_{0,1}(t)$
$\displaystyle=ab^{*}\sum_{k=0}^{\infty}\braket{0,k}{e^{-i\hat{H}_{\text{Ens}}t}}{0,0}\braket{1,0}{e^{i\hat{H}_{\text{Ens}}t}}{1,k}$
(S39)
$\displaystyle=ab^{*}e^{-iE_{0}t}\braket{1,0}{e^{i\hat{H}_{\text{Ens}}t}}{1,0},$
(S40)

and then note that $\braket{1,0}{e^{i\hat{H}_{\text{Ens}}t}}{1,0}$ is the
characteristic function of the random variable $\Lambda=w\lambda+E_{1}$ Lukacs
(1972). We denote the characteristic function of a random variable $X$ as
$\varphi_{X}$ such that we may now write

$c_{0,1}(t)=ab^{*}e^{-i(E_{0}-E_{1})t}\varphi_{w\lambda}(t),$ (S41)

where we have made use of the fact that
$\varphi_{\Lambda}(t)=e^{iE_{1}t}\varphi_{w\lambda}(t)$. As
$c_{1,0}(t)=c_{0,1}^{*}(t)$, we now only need to determine $c_{1,1}(t)$, which
will reduce as

$c_{1,1}(t)=|b|^{2}\sum_{k=0}^{\infty}\left|\braket{1,k}{e^{-i\hat{H}_{\text{Ens}}t}}{1,0}\right|^{2}.$
(S42)

We might already expect $c_{1,1}(t)$ to be constant for all times as we know
that populations of the $\ket{0,\lambda}$, $\ket{1,\lambda}$ states are
conserved quantities. It is however not so obvious that this should be the
case when considering Eq. (S42).

Notice that we can write the $\ket{1,k}$ elements in terms of $\ket{1,0}$ in a
compact form if we allow for the corresponding orthonormal polynomials to take
operator arguments, i.e. we have

$\ket{1,k}=p_{k}(\hat{H}_{\text{Ens}})\ket{1,0}.$ (S43)

We then use this fact in conjunction to
$\hat{H}_{\text{Ens}}e^{-i\hat{H}_{\text{Ens}}t}=i\partial_{t}e^{-i\hat{H}_{\text{Ens}}t}$
in order to write

$\displaystyle\braket{1,k}{e^{-i\hat{H}_{\text{Ens}}t}}{1,0}$
$\displaystyle=\braket{1,0}{p_{k}(\hat{H}_{\text{Ens}})e^{-i\hat{H}_{\text{Ens}}t}}{1,0}$
(S44)
$\displaystyle=p_{k}(i\partial_{t})\braket{1,0}{e^{-i\hat{H}_{\text{Ens}}t}}{1,0}$
(S45) $\displaystyle=p_{k}(i\partial_{t})\varphi_{\Lambda}^{*}(t),$ (S46)

and in general we have

$\braket{1,k}{e^{i\hat{H}_{\text{Ens}}t}}{1,k^{\prime}}=p_{k}(-i\partial_{t})p_{k^{\prime}}(-i\partial_{t})\varphi_{\Lambda}(t).$
(S47)

We can now rewrite Eq. (S42) in terms of the characteristic function as

$c_{1,1}(t)=|b|^{2}\sum_{k=0}^{\infty}\left|p_{k}(i\partial_{t})\varphi_{\Lambda}^{*}(t)\right|^{2}.$
(S48)

In this form we can now readily consider $i\partial_{t}c_{1,1}(t)$ and apply
the three-term recurrence relation for $P_{k}$ in order to see that
$c_{1,1}(t)$ is indeed a constant. We have

$\displaystyle i\partial_{t}c_{1,1}(t)$
$\displaystyle=|b|^{2}\sum_{k=0}^{\infty}\frac{1}{\zeta_{k}}\left(\left[i\partial_{t}P_{k}(i\partial_{t})\varphi_{\Lambda}^{*}(t)\right]P_{k}(-i\partial_{t})\varphi_{\Lambda}(t)-P_{k}(i\partial_{t})\varphi_{\Lambda}^{*}(t)\left[-i\partial_{t}P_{k}(-i\partial_{t})\varphi_{\Lambda}(t)\right]\right)$
(S49)
$\displaystyle=|b|^{2}\sum_{k=0}^{\infty}\frac{1}{\zeta_{k}}\left(\left[P_{k+1}(i\partial_{t})+\alpha_{k}P_{k}(i\partial_{t})+\beta_{k}P_{k-1}(i\partial_{t})\right]\varphi_{\Lambda}^{*}(t)P_{k}(-i\partial_{t})\varphi_{\Lambda}(t)\right.$
$\displaystyle\left.\hskip 54.06006pt-
P_{k}(i\partial_{t})\varphi_{\Lambda}^{*}(t)\left[P_{k+1}(-i\partial_{t})+\alpha_{k}P_{k}(-i\partial_{t})+\beta_{k}P_{k-1}(-i\partial_{t})\right]\varphi_{\Lambda}(t)\right)$
(S50)
$\displaystyle=2|b|^{2}\sum_{k=0}^{\infty}\text{Im}\left(\frac{1}{\zeta_{k}}P_{k+1}(i\partial_{t})\varphi_{\Lambda}^{*}(t)P_{k}(-i\partial_{t})\varphi_{\Lambda}(t)-\frac{1}{\zeta_{k-1}}P_{k}(i\partial_{t})\varphi_{\Lambda}^{*}(t)P_{k-1}(-i\partial_{t})\varphi_{\Lambda}(t)\right)$
(S51)
$\displaystyle=2|b|^{2}\lim_{k\to\infty}\frac{1}{\zeta_{k}}\text{Im}\left(P_{k+1}(i\partial_{t})\varphi_{\Lambda}^{*}(t)P_{k}(-i\partial_{t})\varphi_{\Lambda}(t)\right),$
(S52)

whereby having collected like-terms we have noted that we have an infinite
series of the form $\sum_{k=0}^{\infty}a_{k}-a_{k-1}=\lim_{k\to\infty}a_{k}$
(with $a_{-1}=0$), so we are left with a single term of the form
$\frac{1}{\zeta_{k}}P_{k+1}(i\partial_{t})\varphi_{\Lambda}^{*}(t)P_{k}(-i\partial_{t})\varphi_{\Lambda}(t)=\sqrt{\beta_{k+1}}\braket{1,k+1}{e^{-i\hat{H}_{\text{Ens}}t}}{1,0}\braket{1,0}{e^{i\hat{H}_{\text{Ens}}t}}{1,k}$
which goes to zero as $k\to\infty$ for a local Hamiltonian at finite times.

Combining the fact that $\partial_{t}c_{1,1}(t)=0$ with the initial condition
$c_{1,1}(0)=|b|^{2}$ we have that $c_{1,1}(t)=|b|^{2}$. We therefore have the
ensemble averaged dynamics to be

$\displaystyle\bar{\rho}(t)=$
$\displaystyle|a|^{2}\ket{0}\bra{0}+|b|^{2}\ket{1}\bra{1}$ (S53)
$\displaystyle+ab^{*}e^{-i(E_{0}-E_{1})t}\varphi_{w\lambda}(t)\ket{0}\bra{1}+\text{h.c.}$

which is in agreement with what has been derived in Kropf  _et al._ (2016).

###  II.3 Example: Disorder dependent initial state

The pure state formalism is able to treat physical situations where the
initial state of each realization is correlated to the disorder. We illustrate
this point here with an example of an ensemble of systems,
$\hat{H}_{\lambda}$, where $\lambda$ follows a Wigner semicircle distribution
$p(\lambda)=\frac{2}{\pi}\sqrt{1-\lambda^{2}}$ such that $\lambda$ can take
values on the range $[-1,1]$. We shall consider the initial state of each
realization to be of the form

$\ket{\psi_{0},\lambda}=\lambda\ket{0,\lambda}+\sqrt{1-\lambda^{2}}\ket{1,\lambda}$
(S54)

The initial state of the ensemble is now written as

$\displaystyle\ket{\Psi_{0}}$
$\displaystyle=\int_{-1}^{1}\text{d}\lambda\sqrt{p(\lambda)}\left(\lambda\ket{0,\lambda}+\sqrt{1-\lambda^{2}}\ket{1,\lambda}\right)$
(S55)
$\displaystyle=\sum_{k=0}^{\infty}d_{0,k}\ket{0,k}+\sum_{k=0}^{\infty}d_{1,k}\ket{1,k},$
(S56)

with coefficients

$\displaystyle d_{0,k}$ $\displaystyle=\int_{-1}^{1}\text{d}\lambda
p(\lambda)U_{k}(\lambda)\lambda,$ (S57) $\displaystyle d_{1,k}$
$\displaystyle=\int_{-1}^{1}\text{d}\lambda
p(\lambda)U_{k}(\lambda)\sqrt{1-\lambda^{2}},$ (S58)

where $U_{k}(\lambda)$ are the Chebyshev polynomials of the second kind Boyd
(1989); Abramowitz  _et al._ (1998). Computing these integrals yields

$\displaystyle d_{0,k}$ $\displaystyle=\frac{1}{2}\delta_{1,k},$ (S59)
$\displaystyle d_{1,k}$ $\displaystyle=-\frac{4}{\pi}\frac{1+\cos(\pi
k)}{(k+3)(k^{2}-1)},$ (S60)

such that our initial ensemble state is discretized as

$\ket{\Psi_{0}}=\frac{1}{2}\ket{0,1}-\frac{8}{\pi}\sum_{k=0}^{\infty}\frac{1}{(2k+3)(4k^{2}-1)}\ket{1,2k}.$
(S61)

Figure S1: Lattice coupling strength behaviour for an ensemble with Gaussian
disorder, with and without an energy cutoff at $5\sigma$. We are able to
simulate dynamics for longer times by introducing an energy cutoff.

##  III Introducing an energy cutoff

The mapping breaks down for cases of disorder which follow a distribution
whose moments are undefined, e.g. Cauchy and Lévy distributions, which then in
turn result in undefined recurrence coefficients. The divergence implicit in
the integral form of $\hat{H}_{\text{Ens}}$ is made explicit in its lattice
form. We can remedy the situation by introducing an energy cutoff which then
leads to well defined recurrence coefficients and recover approximate averaged
dynamics whose accuracy depends on the choice of cutoff.

It is also useful to apply an energy cutoff even in cases where all moments of
the distribution are well defined, as doing so will ensure that the recurrence
coefficients will reach asymptotic values. Take for example a Gaussian
distribution with mean zero and standard deviation $\sigma$ such that the
recurrence coefficients will be $\alpha_{k}=0$ for all $k$ and
$\beta_{k+1}=\sigma^{2}(k+1)$. Due to the strictly increasing value of the
lattice couplings we find that there are diminishing returns in increasing the
lattice truncation, making it more and more computationally intensive to
simulate dynamics to longer times. A $5\sigma$ cutoff will lead to the
couplings having an asymptotic value of
$\lim_{k\to\infty}\sqrt{\beta_{k+1}}=\frac{5}{2}\sigma$, thus curbing the
monotonic increase of the lattice couplings whilst simultaneously preserving
the fidelity of the distribution up to $\pm 5\sigma$. Figure S1 demonstrates
this behaviour of the lattice couplings corresponding to a Gaussian disorder
distribution with and without an energy cutoff.

Figure S2: Disordered qubit ensemble dephasing dynamics error.

##  IV Table of quantities

Name | $p(\lambda)$ | $\mathcal{D}$ | $\alpha_{k}$ | $\beta_{k+1}$ | $\varphi_{\Lambda}(t)$
---|---|---|---|---|---
Gaussian | $\frac{1}{\sqrt{2\pi}\sigma}e^{-\lambda^{2}/2\sigma^{2}}$ | $\mathbb{R}$ | 0 | $\sigma^{2}(k+1)$ | $e^{-\sigma^{2}t^{2}/2}$
Cauchy | $\frac{1}{\pi\theta}\frac{\theta^{2}}{\lambda^{2}+\theta^{2}}$ | $\mathbb{R}$ | N/A | N/A | $e^{-\theta|t|}$
Semicircle | $\frac{2}{\pi w^{2}}\sqrt{w^{2}-\lambda^{2}}$ | $[-w,w]$ | 0 | $\frac{w}{2}$ | $2\frac{J_{1}(wt)}{wt}$
Uniform | $\frac{1}{2v}$ | $[-v,v]$ | 0 | $\frac{v(k+1)}{\sqrt{4(k+1)^{2}-1}}$ | $\frac{\sin(vt)}{vt}$

For the qubit dephasing considered in the main text we have set
$\sigma=\theta=w=v=E_{1}-E_{0}$. An energy cutoff of $E_{\pm}=\pm 30\theta$ is
introduced for the Cauchy distribution in order to achieve well defined
recurrence coefficients.

## References

  * Anderson (1958) P. W. Anderson, Phys. Rev. 109, 1492 (1958).
  * Gornyi  _et al._ (2005) I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Phys. Rev. Lett. 95, 206603 (2005).
  * Paredes  _et al._ (2005) B. Paredes, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 95, 140501 (2005).
  * Basko  _et al._ (2006) D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Ann. Phys. 321, 1126 (2006).
  * Bermudez  _et al._ (2010) A. Bermudez, M. A. Martin-Delgado, and D. Porras, New J. Phys. 12, 123016 (2010).
  * Pal and Huse (2010) A. Pal and D. A. Huse, Phys. Rev. B 82, 174411 (2010).
  * Bardarson  _et al._ (2012) J. H. Bardarson, F. Pollmann, and J. E. Moore, Phys. Rev. Lett. 109, 017202 (2012).
  * Breitweiser  _et al._ (2020) S. A. Breitweiser, A. L. Exarhos, R. N. Patel, J. Saouaf, B. Porat, D. A. Hopper, and L. C. Bassett, ACS Photonics 7, 288 (2020).
  * Maldacena and Stanford (2016) J. Maldacena and D. Stanford, Phys. Rev. D 94, 106002 (2016).
  * Pollack  _et al._ (2020) J. Pollack, M. Rozali, J. Sully, and D. Wakeham, Phys. Rev. Lett. 125, 021601 (2020).
  * Peng (2021) C. Peng, Phys. Rev. D 103, L061901 (2021).
  * Schlawin  _et al._ (2012) F. Schlawin, N. Cherroret, and A. Buchleitner, Europhys. Lett. 99, 14001 (2012).
  * Gilead  _et al._ (2015) Y. Gilead, M. Verbin, and Y. Silberberg, Phys. Rev. Lett. 115, 133602 (2015).
  * Bai  _et al._ (2016) Y. F. Bai, P. Xu, L. L. Lu, M. L. Zhong, and S. N. Zhu, J. Opt. 18, 055201 (2016).
  * Lib and Bromberg (2022) O. Lib and Y. Bromberg, Nat. Phys. 18, 986 (2022).
  * Popescu  _et al._ (2006) S. Popescu, A. J. Short, and A. Winter, Nat. Phys. 2, 754 (2006).
  * Cotler  _et al._ (2017) J. Cotler, N. Hunter-Jones, J. Liu, and B. Yoshida, J. High Energy Phys. 2017, 1 (2017).
  * Nation and Porras (2018) C. Nation and D. Porras, New J. Phys. 20, 103003 (2018).
  * Foini and Kurchan (2019) L. Foini and J. Kurchan, Phys. Rev. E 99, 042139 (2019).
  * Nation and Porras (2020) C. Nation and D. Porras, Phys. Rev. E 102, 042115 (2020).
  * Dabelow and Reimann (2020) L. Dabelow and P. Reimann, Phys. Rev. Lett. 124, 120602 (2020).
  * Frost  _et al._ (2010) J. M. Frost, M. A. Faist, and J. Nelson, Adv. Mater. 22, 4881 (2010).
  * Gelin  _et al._ (2021) M. F. Gelin, A. Velardo, and R. Borrelli, J. Chem. Phys. 155 (2021), 10.1063/5.0065896.
  * Olaya-Castro  _et al._ (2008) A. Olaya-Castro, C. F. Lee, F. F. Olsen, and N. F. Johnson, Phys. Rev. B 78, 085115 (2008).
  * Kolli  _et al._ (2012) A. Kolli, E. J. O’Reilly, G. D. Scholes, and A. Olaya-Castro, J. Chem. Phys. 137, 174109 (2012).
  * Stross  _et al._ (2016) C. Stross, M. W. Van der Kamp, T. A. A. Oliver, J. N. Harvey, N. Linden, and F. R. Manby, J. Phys. Chem. B 120, 11449 (2016).
  * Cohn  _et al._ (2022) B. Cohn, S. Sufrin, A. Basu, and L. Chuntonov, J. Phys. Chem. Lett. 13, 8369 (2022).
  * Gneiting  _et al._ (2016) C. Gneiting, F. R. Anger, and A. Buchleitner, Phys. Rev. A 93, 032139 (2016).
  * Kropf  _et al._ (2016) C. M. Kropf, C. Gneiting, and A. Buchleitner, Phys. Rev. X 6, 031023 (2016).
  * Chen  _et al._ (2018) H.-B. Chen, C. Gneiting, P.-Y. Lo, Y.-N. Chen, and F. Nori, Phys. Rev. Lett. 120, 030403 (2018).
  * Gneiting (2020) C. Gneiting, Phys. Rev. B 101, 214203 (2020).
  * Breuer and Petruccione (2007) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Bulla  _et al._ (2005) R. Bulla, H.-J. Lee, N.-H. Tong, and M. Vojta, Phys. Rev. B 71, 045122 (2005).
  * Bulla  _et al._ (2008) R. Bulla, T. A. Costi, and T. Pruschke, Rev. Mod. Phys. 80, 395 (2008).
  * Prior  _et al._ (2010) J. Prior, A. W. Chin, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 105, 050404 (2010).
  * Chin  _et al._ (2010) A. W. Chin, Á. Rivas, S. F. Huelga, and M. B. Plenio, J. Math. Phys. 51, 092109 (2010).
  * Ismail (2005) M. E. H. Ismail,  _Classical and Quantum Orthogonal Polynomials in One Variable_ (Cambridge University Press, Cambridge, England, UK, 2005).
  * Teschl (2000) G. Teschl,  _Jacobi Operators and Completely Integrable Nonlinear Lattices (Mathematical Surveys and Monographs)_, Vol. 72 (American Mathematical Society, 2000).
  * de Vega  _et al._ (2015) I. de Vega, U. Schollwöck, and F. A. Wolf, Phys. Rev. B 92, 155126 (2015).
  * Parker  _et al._ (2019) D. E. Parker, X. Cao, A. Avdoshkin, T. Scaffidi, and E. Altman, Phys. Rev. X 9, 041017 (2019).
  * Mück and Yang (2022) W. Mück and Y. Yang, Nucl. Phys. B 984, 115948 (2022).
  * Kropf (2020) C. M. Kropf, Phys. Rev. Res. 2, 033311 (2020).
  * Lukacs (1972) E. Lukacs, Adv. Appl. Probab. 4, 1 (1972).
  * Cogdell  _et al._ (2006) R. J. Cogdell, A. Gall, and J. Köhler, Q. Rev. Biophys. 39, 227 (2006), 17038210 .
  * Novoderezhkin  _et al._ (2003) V. Novoderezhkin, M. Wendling, and R. van Grondelle, J. Phys. Chem. B 107, 11534 (2003).
  * Boyd (1989) J. P. Boyd,  _Chebyshev & Fourier Spectral Methods_ (Springer, Berlin, Germany, 1989).
  * Abramowitz  _et al._ (1998) M. Abramowitz, I. A. Stegun, and R. H. Romer,  _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_, Vol. 56 (American Association of Physics Teachers, 1998) Chap. 22.