# A Universal End-to-End Approach to
Portfolio Optimization via Deep Learning

Chao Zhang, Zihao Zhang††footnotemark: , Mihai Cucuringu, Stefan Zohren
University of Oxford Equal contribution. Correspondence to: Zihao Zhang
<EMAIL_ADDRESS>
###### Abstract

We propose a universal end-to-end framework for portfolio optimization where
asset distributions are directly obtained. The designed framework circumvents
the traditional forecasting step and avoids the estimation of the covariance
matrix, lifting the bottleneck for generalizing to a large amount of
instruments. Our framework has the flexibility of optimizing various objective
functions including Sharpe ratio, mean-variance trade-off etc. Further, we
allow for short selling and study several constraints attached to objective
functions. In particular, we consider cardinality, maximum position for
individual instrument and leverage. These constraints are formulated into
objective functions by utilizing several neural layers and gradient ascent can
be adopted for optimization. To ensure the robustness of our framework, we
test our methods on two datasets. Firstly, we look at a synthetic dataset
where we demonstrate that weights obtained from our end-to-end approach are
better than classical predictive methods. Secondly, we apply our framework on
a real-life dataset with historical observations of hundreds of instruments
with a testing period of more than 20 years.

##  1 Introduction

Portfolio optimization is the process of selecting the best asset distribution
to invest in out of all the considered assets, with the aim of optimizing a
suitably chosen objective function, such as maximizing returns while
minimizing the level of risk. The theory, also known as modern portfolio
theory (MPT), was introduced by Markowitz [23] and has become one of the
corner stones of quantitative finance. The main benefit of MPT stems from
diversification, as the variance of a portfolio decreases with the addition of
less correlated instruments, leading to improvements in the portfolio’s Sharpe
Ratio [30], one of the popular performance metrics employed by practitioners.
Higher Sharpe Ratios can be combined with leverage, leading to higher returns
for a given level of risk preference.

Despite the wide popularity and practical usage of MPT, the classical mean-
variance (MV) optimization approach employed in MPT severely suffers from a
number of limitations. In particular, the traditional method of optimizing a
portfolio is a two-step optimization problem, where the first step aims to
predict future returns, and a constrained optimization problem is then solved
to derive the optimal portfolio weights. However, financial returns are
notoriously stochastic with extremely low signal-to-noise ratio. The
construction of predictive signals remains extremely difficult, and the
generalization and robustness of such signals is often questioned.
Furthermore, the MV method involves the estimation of the covariance matrix
for all the considered instruments, and computing its inverse [26]. However,
such estimations are highly unstable, due to the large number of pairwise
coefficients to estimate, which hinders the scalability of the method and can
lead to concentration risk111Portfolio weights are only assigned to a few
instruments, which decreases the diversification power..

In this work, we propose an end-to-end framework that directly optimizes a
portfolio by utilizing deep learning models [11]. Compared to the classical
two-step method, we directly optimize portfolio weights circumventing the
requirement of predicting future returns. Furthermore, our framework enjoys
the flexibility of being able to optimize _various objective functions_ ,
including global minimum variance, maximum Sharpe Ratio, mean-variance trade-
off, etc. This flexibility allows us to study portfolios with distinct
characteristics, and to adjust risk preference accordingly. Given the chosen
objective function, our setup implicitly takes into account the estimation of
the covariance matrix, thus ensuring the scalability of our approach.

Our framework allows for short-selling, and we study _different constraints_
attached to objective functions. Specifically, we consider cardinality222We
limit the number of traded assets in a portfolio., maximum position333We limit
the maximum position for individual instrument preventing concentration risk.,
short-selling and leverage. The general MV approach is a standard quadratic
programming problem, which can be solved exactly optimally. However, the
addition of constraints can increase the degree of complexity; for example,
the cardinality constraint renders the problem non-convex, and no longer
solvable by exact methods [14]. Our approach allows us to incorporate these
constraints into objective functions by utilizing several neural layers, and
gradient ascent is adopted for optimization, thus naturally dealing with
issues due to the non-convexity of the solution space.

We test our method by constructing a portfolio with hundreds of instruments,
and considering a testing period of more than 20 years. The experiments show
superior performance for our approach, ensuring the scalability and
generalization of the end-to-end framework. We also investigate the model
performance under realistic transaction costs, and carry out comparisons to
discuss the usage of different deep neural networks. In addition, we compare
the classical two-step MV problem with our approach on a simulated data set,
where the optimal portfolio weights are known. The results indicate that our
method delivers solutions that are closer to optimal weights, showcasing the
efficiency of our setup.

#### Outline

The remainder of this paper is structured as follows: A literature review is
included to introduce related works in Section 2. Section 3 presents our
method and details various constraints. We then describe our experiments in
Section 4, and conclude the findings in Section 5, along with future research
directions.

##  2 Literature Review

The portfolio optimization theory, formally introduced by Markowitz [23],
plays a significant role in both research and practice. The basic MV model
considers only long positions and requires budget constrains; for example, the
sum of weights equals to 1. It can be formulated as a quadratic programming
problem [9], and the solution only depends on the expected mean and covariance
matrix of asset returns. The simplest way of estimating the expected mean and
covariance matrix is to use the sample mean and covariance matrix. However,
the gain from plugging the sample estimates into the optimization problem is
offset by the error between population parameters and their sample estimates,
also known as “Markowitz optimization enigma", revealed by [24, 18, 16].

To resolve the “Markowitz optimization enigma", a major direction is to find
better estimators of expected means, such as Bayes and shrinkage estimators
[15, 18, 27]. However, the dynamics of financial returns is rather stochastic
and the prediction of returns is a highly non-trivial task. There also exists
a large volume of literature applying machine learning models to predict asset
returns, including random forest [1], vanilla neural networks [29], long
short-term memory (LSTM) networks [8, 7], and others. Furthermore, the
estimation of the covariance matrix is unstable when we have a large
collection of variables [20, 21]. To remedy this, our framework directly
outputs portfolio weights and optimizes objective functions bypassing the
estimation of expected means and the covariance matrix of returns.

Another drawback of the basic MV model is that it neglects financial
restrictions that are essential in practice. For example, the cardinality
constraint limits the number of assets to be included in a portfolio and the
requirement of maximum position for individual instrument is necessary for
diversification. In the literature, heuristic mechanisms [13, 5] are proposed
to study the cardinality constraint, such as genetic algorithms, simulated
annealing, and branch-and-bound algorithms [10, 19]. With respect to the
maximum position constraint, it can be solved efficiently by some specialized
methods such as the simplex method [35]. However, the computational complexity
of the portfolio optimization problem becomes much greater than the basic
model when considering these restrictions. In particular, the cardinality
constraint makes the problem fall into the class of NP-hard problems [31, 5].
In our case, we formulate these constraints by utilizing corresponding
specialized neural layers, including a differentiable sorting algorithm [12,
6, 2], paving the way for gradient ascent to be used for training, and thus
avoiding the computation difficulties.

The idea of combining the prediction and optimization tasks in an end-to-end
model has been previously explored in the context of portfolio optimization
[37, 4, 34], as well as in other contexts, such as momentum strategies [22].
Zhang et al. [37] adopt neural networks to directly optimize the portfolio’s
Sharpe Ratio by circumventing the prediction of future returns. However, the
portfolio in [37] only considers the basic MV problem, i.e. requiring long-
only and budget-constrained. Recent work of [4] uses neural networks with
differentiable optimization layers to find the solutions of MV problems with
certain constraints. The work of [34] learns the optimal risk contribution on
each asset and allocates the weights according to the risk-budgeting strategy
[3]. It is worth to emphasize that all aforementioned portfolios can be
optimized and unified in our proposed framework, and we provide a
comprehensive study of end-to-end frameworks by optimizing various objective
functions with several practical constraints.

##  3 Methods

In this section, we introduce our end-to-end framework and start with a
general mean-variance problem (MVP) that allows for _short-selling_

$\begin{split}\max_{\mathbf{w}_{t}}\,\,\mathbb{E}(r_{p,t+1})-\frac{\lambda}{2}\operatorname{Var}(r_{p,t+1}),&\text{\quad
s.t.\quad}\|\mathbf{w}_{t}\|_{1}=\sum_{i=1}^{N}|w_{i,t}|=1,\\\
r_{p,t+1}&={\mathbf{w}_{t}}^{\prime}\mathbf{r}_{t+1},\end{split}$ (1)

where $r_{p,t+1}$ is the portfolio return and $\lambda$ is the risk aversion
rate. $\mathbf{w}_{t}=(w_{1,t},\dots,w_{N,t})^{\prime}$ represents the
fraction of the portfolio value invested in the each asset at time $t$.
$\mathbf{r}_{t}=(r_{1,t},\dots,r_{N,t})^{\prime}$ is the collection of returns
of $N$ assets at time $t$ and $r_{i,t}$ denotes the percentage return of asset
$i\,\,(i=1,\dots,N)$. Note that $t$ can denote any arbitrary period, such as
one-minute intervals, days, months, etc.

#### Objective functions.

Besides the MVP defined in Equation (1), we study the following objective
functions

  * •

Global Minimum Variance Portfolio (GMVP)

$\min_{\mathbf{w}_{t}}\,\operatorname{Var}(r_{p,t+1}),$ (2)

  * •

Maximum Sharpe Ratio Portfolio (MSRP)

$\max_{\mathbf{w}_{t}}\,\frac{\mathbb{E}(r_{p,t+1})}{\operatorname{Var}(r_{p,t+1})},$
(3)

Note that we let all objective functions have the same constraints as in MVP
(Equation (1)), and allow for short-selling.

#### Constraints.

In terms of the constraints, we consider the following cases

  * •

Long-only (LONG): $w_{i,t}\geq 0$, as short-selling is forbidden in certain
markets, for example, in Chinese stock markets,

  * •

Maximum position (MAX): $w_{i,t}\leq u$ specifies the upper bound $u$ (with
$u>0$) of weights allocated to each asset in a portfolio which helps with
concentration risk, ensuring diversification,

  * •

Cardinality (CAR): $\|\mathbf{w}_{t}\|_{0}=K$ restricts the number of assets
in a portfolio,

  * •

Leverage (LEV): $\|\mathbf{w}_{t}\|_{1}=L$ with $L\geq 1$ allows investors to
increase the exposure to markets.

###  3.1 The Proposed End-to-End Framework

We denote by $\mathcal{R}_{t}$ the set of current information, for example,
previous returns $\\{\mathbf{r}_{1},\dots,\mathbf{r}_{t}\\}$. The classical
approach can be essentially divided into two parts: prediction and
optimization. The predictive step is to estimate future returns by
$\hat{r}_{i,t+1}=f\left(\mathcal{R}_{t}\right)$ and the predicted returns are
then inserted into the objective function in Equation (1) to derive the
portfolio weights $\mathbf{w}_{t}$. In this work, we directly model the
portfolio weights $\mathbf{w}_{t}$ with deep neural networks ($f$) as

$\mathbf{w}_{t}=f\left(\mathcal{R}_{t}\right),$ (4)

where $\mathcal{R}_{t}$ denotes the set of current information as inputs to
the networks. Figure 1 depicts our proposed end-to-end framework, which is
composed of two ingredients: a score block and a portfolio block.

Figure 1: Architecture of our proposed End-to-End approach. $h_{1}$ represents
a neural network, which transforms the input into fitness scores for each
asset. $h_{2}$ represents the specialized differential layers to convert the
scores to the portfolio weights satisfying the required constraints.

#### Score block.

The input $\mathcal{R}_{t}$ consists of current market information, for
example, returns $\mathbf{r}_{t-p},\dots,\mathbf{r}_{t}$ up to lag $p$. We use
a neural network to transform the input into fitness scores for each asset. A
higher score implies that the corresponding asset is more likely to be
assigned higher weights (without cardinality constraints) or selected (with
cardinality constraints), and vice versa. We denote this network function as
$h_{1}$ and the fitness scores as

$\mathbf{s}_{t}=(s_{1,t},\dots,s_{N,t})^{\prime}=h_{1}(\mathcal{R}_{t}).$ (5)

#### Portfolio block.

In the portfolio block, we convert the previous obtained scores
$\mathbf{s}_{t}$ into the portfolio weights satisfying the required
constraints through corresponding differential functions $h_{2}(\cdot)$. We
then compute the realized portfolio returns $r_{p,t+1}$ based on future
returns $\mathbf{r}_{t+1}$ and calculate the loss value according to the
selected objective function, for example, Equation (1). To allow for short
selling and to meet the requirement of the absolute values of the weights to
sum to one, we transform the fitness scores as follows,

$\begin{split}w_{i,t}&=h_{2}(s_{i,t})\\\
&=\text{sign}(s_{i,t})\times\frac{e^{s_{i,t}}}{\sum_{j=1}^{N}e^{s_{j,t}}}.\end{split}$
(6)

###  3.2 The Design of the Differential Function $h_{2}$ for Constraints

Similar to Equation (6), we now demonstrate how the differential function
$h_{2}$ can be constructed in order to meet the previous discussed
constraints, by considering the fitness scores $\bm{s}_{t}$ from the score
block. We first showcase how a single constraint can be incorporated by re-
parameterizing the scores, and then demonstrate how multiple constraints can
be combined.

#### Re-parametrisation of $\mathbf{s}_{t}$ for a single constraint

Constraint (1) Long-only and $\|\mathbf{w}_{t}\|_{1}=1$: We apply the
$\operatorname{softmax}$ activation function to the scores, in order to obtain
the non-negative weights. For $i=1,\ldots,N$, we consider

$w_{i,t}=\frac{e^{s_{i,t}}}{\sum_{j=1}^{N}e^{s_{j,t}}}.$ (7)

Constraint (2) Maximum and $\|\mathbf{w}_{t}\|_{1}=1$: We use a generalized
sigmoid function, $\phi_{a}(x)=a+\frac{1}{1+e^{-x}}$ (with $a\geq 0$) to
transform the scores, so that weights can automatically meet the upper bound
constraint $u$. Upon setting $a=\frac{1-u}{Nu-1}$, we obtain $w_{i,t}\leq u$
as follows

$w_{i,t}=\operatorname{sign}(s_{i,t})\times\frac{\phi_{a}(|s_{i,t}|)}{\sum_{j=1}^{N}\phi_{a}(|s_{j,t}|)}.$
(8)

If we let the maximum position $u$ equal to 1, $\phi_{a}(x)$ degenerates to
the sigmoid function.

Constraint (3) Cardinality and $\|\mathbf{w}_{t}\|_{1}=1$: To achieve
cardinality, we first employ a sorting operator $\Pi(\cdot)$ that maps
$\mathbf{s}_{t}\in\mathbb{R}^{N}$ into a permutation matrix
$\Pi(\mathbf{s}_{t})\in\mathbb{R}^{N\times N}$, such that
$\tilde{\mathbf{s}}_{t}=\Pi(\mathbf{s}_{t})$ is the sorted vector of
$\mathbf{s}_{t}$ in descending order. We then long the top $n$ instruments and
short the bottom $n$ instruments to construct our portfolio

$\displaystyle\begin{split}w_{i,t}&=\frac{1}{2}\times\frac{{1}_{\\{{s}_{i,t}>d_{u}\\}}e^{|s_{i,t}|}}{\sum_{j=1}^{N}{1}_{\\{{s}_{j,t}>d_{u}\\}}e^{|s_{j,t}|}}-\frac{1}{2}\times\frac{{1}_{\\{{s}_{j,t}<d_{l}\\}}e^{|s_{i,t}|}}{\sum_{j=1}^{N}{1}_{\\{{s}_{i,t}<d_{l}\\}}e^{|s_{j,t}|}},\\\
d_{u}&=\tilde{\mathbf{s}}_{t}[n],\quad d_{l}=\tilde{\mathbf{s}}_{t}[N-n],\quad
n=\lfloor K/2\rfloor+1,\\\ \end{split}$ (9)

where $\tilde{\mathbf{s}}_{t}[n]$ means the $n$-th element in vector
$\tilde{\mathbf{s}}_{t}$, i.e. the $n$-th largest element of $\mathbf{s}_{t}$.
To compute the sorting operator, we first define a square matrix
$\Lambda_{i,j}^{t}$ on the fitness score $\mathbf{s}_{t}$ as

$\Lambda_{i,j}^{t}=(N+1-2i)s_{j,t}-\sum_{m}|s_{j,t}-s_{m,t}|.$ (10)

According to previous works [25, 12, 2, 6], the permutation matrix
$\Pi(\mathbf{s}_{t})$ can be constructed as

$\Pi(\mathbf{s}_{t})_{i,j}=\begin{cases}1,&\text{if}\
j=\operatorname{argmax}(\Lambda_{i,:}^{t})\\\ 0,&\text{otherwise}.\end{cases}$
(11)

Since the $\operatorname{argmax}$ function is not differentiable, the authors
of [12] proposes a NeuralSort layer that replaces the $\operatorname{argmax}$
operator by the $\operatorname{softmax}$ operator, to arrive at a relaxed
version of $\Pi(\mathbf{s}_{t})$

$\widehat{\Pi(\mathbf{s}_{t})}_{i,:}=\operatorname{softmax}(\Lambda_{i,:}^{t}).$
(12)

Then Equation (9) is differentiable and the standard gradient descent method
can be applied.

Constraint (4) Leverage, i.e. $\|\mathbf{w}_{t}\|_{1}=L$: Similar to Equation
(6), we increase the total exposure of the positions by a factor of $L$

$w_{i,t}=L\times\operatorname{sign}(s_{i,t})\times\frac{e^{|s_{i,t}|}}{\sum_{j=1}^{N}e^{|s_{j,t}|}}.$
(13)

#### Re-parametrisation of $\mathbf{s}_{t}$ for multiple constraints.

We handle multiple constraints by combining the corresponding techniques in
the above cases. Consider the constraints Maximum & Cardinality & Leverage as
an example, using the techniques in the Cases (2), (3) and (4), we have, for
$i=1,\ldots,N$,

$\displaystyle\begin{split}w_{i,t}&=\frac{L}{2}\times\frac{{1}_{\\{{s}_{i,t}>d_{u}\\}}\phi_{a}(|s_{i,t}|)}{\sum_{j=1}^{N}{1}_{\\{{s}_{j,t}>d_{u}\\}}\phi_{a}(|s_{j,t}|)}-\frac{L}{2}\times\frac{{1}_{\\{{s}_{i,t}<d_{l}\\}}\phi_{a}(|s_{i,t}|)}{\sum_{j=1}^{N}{1}_{\\{{s}_{j,t}<d_{l}\\}}\phi_{a}(|s_{j,t}|)},\\\
d_{u}&=\tilde{\mathbf{s}}_{t}[n],\quad d_{l}=\tilde{\mathbf{s}}_{t}[N-n],\quad
n=\lfloor K/2\rfloor+1,\\\ \end{split}$ (14)

where $a=\frac{1-u/2L}{nu/2L-1}$, and one may analogously perform the re-
parametrisation for other combinations of constraints.

As an example, we can consider a long-short equity portfolio. This is a
popular strategy where we go long a certain percentile of top performing
stocks, and short the same percentile of bottom stocks. Commonly, in the
literature, one first predicts the assets’ returns and sorts assets in order.
Interestingly, this strategy can be achieved through our framework by
incorporating the third Cardinality constraint. As fitness scores are sorted
from the _Score block_ in Figure 1, we can simply long the top best assets and
short the bottom worst ones. This also provides an alternative to recently
proposed learning-to-rank algorithms to construct such long-short strategies
[28].

##  4 Experiments

###  4.1 Descriptions of Data sets

In this section, we demonstrate the efficacy of our methods on two data sets.
The first one is a synthetic data set generated from multivariate normal
distributions, where we compare the classical two-step method with our end-to-
end framework, in terms of the similarity of the optimal weights. The second
data set contains daily observations for 735 stocks from the Wharton Research
Data Services (WRDS) [36], and we test the generalization of our model with a
testing period of more than 20 years.

#### Synthetic data set.

We simulate daily returns from multivariate normal distributions with
parameters calibrated from real data. Specifically, in terms of the mean and
covariance matrix, we use the sample mean and (shrunk) sample covariance
matrix corresponding to each year from 1984 to 2020 of the WRDS data.

#### The WRDS data set

We extract daily observations for 735 stocks from the Russell 3000 Index,
ranging from 1984/01/03 to 2021/07/06. A rolling window approach is adopted
for training; in our case, we consider data between 1984 and 1999 as the first
training set, 2000 as the validation set, and 2001 as the testing set. We then
roll forward one year and repeat the process until the end date, gradually
increasing the training size. Overall, our testing period ranges from 2001 to
2021.

###  4.2 Training and Testing Procedure

We denote our end-to-end framework as E2E and compare it with the classical
setup (CS). Several models are tested in the experiments, including a linear
model (LM), a single fully connected neural network layer with 64 units, also
referred to as multilayer perceptron (MLP), a single LSTM layer with 64 units
and a convolutional neural network (CNN) model. The CNN model consists of 4
layers where the first 3 layers are one-dimensional convolutional layers with
filters of size 32, 64, 128, and each filter has the same kernel size (3,1).
The last layer of the CNN model is a single LSTM with 64 units. We set the
learning rate to $10^{-4}$, batch size to $64$ and number of epochs to $1000$.

When reporting performance during the testing period, we include transaction
costs and calculate the portfolio return as

$\begin{split}R_{p,t}&=\sum_{i=1}^{N}R_{i,t}\\\
R_{i,t}&=w_{i,t-1}r_{i,t}-C|w_{i,t}(1+w_{i,t-1}r_{i,t})-w_{i,t-1}|,\end{split}$
(15)

where $R_{i,t}$ is the trade return for the instrument $i$ and $r_{i,t}$ its
price return at time $t$. We denote the cost rate, $C$, in terms of basis
points ($1bp=10^{-4}$) and report the following metrics to evaluate model
performance

  * •

E(R): annualized expected return,

  * •

Std(R): annualized standard deviation,

  * •

Sharpe: annualized Sharpe Ratio,

  * •

DD(R): annualized downside deviation, where only the standard deviation of
negative returns is calculated,

  * •

Sortino: annualized Sortino Ratio,

  * •

MDD: maximum drawdown,

  * •

% of +Ret: the ratio between positive and negative returns,

  * •

Frobenius norm: the similarity between the estimated portfolio weights and the
optimal weights. The smaller the value, the better the estimation.

  * •

Turnover: an indicator that represents the turnover rate,

  * •

Beta: it measures a portfolio’s correlation relative to S&P 500 Index.

###  4.3 Experimental Results for the Synthetic Data Set

In Table 2, we present the results for the synthetic data set when maximizing
the Sharpe Ratio, with the constraint requiring that the summation of the
absolute values of the weights equals to one. We denote the classical two-step
method as CS, and our end-to-end approach as E2E. Each framework is tested by
four models described in Section 4.2, and we include a naive estimation for CS
by calculating the predicted return with just past returns (CS-SAMPLE).

We focus on Sharpe Ratio and Frobenius norm as the main metrics. The Sharpe
Ratio essentially describes the return per unit of risk and our E2E approach
demonstrates superior results compared to the CS framework. Additionally, all
E2E models deliver better portfolio weights in terms of Frobenius norm,
meaning that the weights derived from the end-to-end setup are closer to the
optimal weights, when compared to the classical predictive approach.

| E(R) | Std(R) | Sharpe | DD(R) | Sortino | MDD | $\frac{\text{\% of}}{\text{+Ret}}$ | $\frac{\text{Frobenius}}{\text{norm}}$ | Turnover
---|---|---|---|---|---|---|---|---|---
MSRP
CS-SAMPLE | 0.011 | 0.018 | 0.639 | 0.012 | 0.932 | 0.072 | 0.517 | 4.695 | 2.006
CS-LM | 0.001 | 0.011 | 0.087 | 0.008 | 0.123 | 0.037 | 0.501 | 4.700 | 9.635
CS-MLP | 0.001 | 0.017 | 0.011 | 0.012 | 0.016 | 0.082 | 0.497 | 4.873 | 9.518
CS-LSTM | 0.001 | 0.011 | 0.053 | 0.008 | 0.074 | 0.086 | 0.504 | 4.684 | 6.913
CS-CNN | 0.002 | 0.011 | 0.160 | 0.008 | 0.227 | 0.053 | 0.500 | 4.685 | 7.214
E2E-LM | 0.004 | 0.031 | 0.146 | 0.022 | 0.206 | 0.119 | 0.513 | 4.253 | 7.203
E2E-MLP | 0.003 | 0.012 | 0.238 | 0.009 | 0.339 | 0.056 | 0.516 | 4.280 | 5.849
E2E-LSTM | 0.014 | 0.012 | 1.216 | 0.008 | 1.810 | 0.008 | 0.529 | 4.407 | 0.692
E2E-CNN | 0.014 | 0.014 | 0.977 | 0.009 | 1.441 | 0.027 | 0.534 | 4.426 | 3.604

Table 2: Experimental results for the synthetic data set.

###  4.4 Experimental Results for the WRDS Data Set

Table 3 shows the experimental results for the WRDS data set, with a testing
period of more than 20 years ranging from 2001 to 2021. The table is split
into four parts, where the first part (Baselines) includes five benchmark
models: the performance of the S&P 500 Index, a portfolio with equal weights
(EWP), the maximum diversification (MD) portfolio [33], the diversity-weighted
portfolio (DWP) from stochastic portfolio theory [17] and the global minimum
variance portfolio (GMVP) [33]. The second block of Table 3 compares our
approach with the classical predictive methods. The objective function is the
same as in the previous Section 4.3, where we maximize the Sharpe Ratio. The
third block studies various other objective functions, and the last block
presents the effects of different constraints. In particular, we demonstrate
the following cases

  * •

MSRP-LONG: a long-only portfolio that maximizes the Sharpe Ratio,

  * •

LEV: a portfolio with leverage, where the leverage is set to 5,

  * •

CAR: a portfolio with a cardinality constraint set to 148, which amounts to
20% of our total number of instruments; in other words, we long the top decile
and short the bottom decile,

  * •

MAX: a portfolio with a maximum position for each individual instrument, which
we set to 5%,

  * •

LEV-CAR-MAX: a combination of above constraints.

We draw the following conclusions. From the second block of Table 3, we
observe that our end-to-end methods deliver better results than the classical
approach and outperform all the baseline models. The third block demonstrates
how different objective functions affect the model performances. Specifically,
GMVP attains the lowest variance, and we can achieve the desired risk level by
adjusting the risk aversion parameter ($\lambda$) in MVP, where a higher
$\lambda$ penalizes the risk more, leading to a portfolio with a smaller
variance. The last block shows the results for different constraints, and we
include them as examples to showcase the flexibility of our method. The exact
parameters of these constraints can be customized by the end-user, in order to
achieve their preferred conditions and adjust to the specificity of their
trading environment. Table 4 has the same layout as Table 3, but we study
model performances under transaction costs ($C=2$ bps).

| E(R) | Std(R) | Sharpe | DD(R) | Sortino | MDD | $\frac{\text{\% of}}{\text{+Ret}}$ | Beta | Turnover
---|---|---|---|---|---|---|---|---|---
Baselines
S&P 500 | 0.061 | 0.196 | 0.402 | 0.140 | 0.563 | 0.568 | 0.541 | 1.000 | 0.000
EWP | 0.130 | 0.212 | 0.682 | 0.148 | 0.973 | 0.548 | 0.546 | 1.000 | 0.000
MD | 0.439 | 0.239 | 1.641 | 0.141 | 2.785 | 0.519 | 0.548 | 0.599 | 2.347
DWP | 0.117 | 0.209 | 0.633 | 0.147 | 0.899 | 0.554 | 0.545 | 1.000 | 0.009
GMVP | 0.080 | 0.081 | 0.992 | 0.059 | 1.360 | 0.408 | 0.564 | 0.257 | 1.577
MSRP
CS-SAMPLE | -0.031 | 0.038 | -0.812 | 0.028 | -1.093 | 0.476 | 0.484 | 0.037 | 1.708
CS-LM | 0.004 | 0.015 | 0.290 | 0.011 | 0.414 | 0.062 | 0.504 | 0.009 | 9.157
CS-MLP | 0.008 | 0.027 | 0.299 | 0.019 | 0.424 | 0.140 | 0.515 | 0.036 | 9.319
CS-LSTM | 0.014 | 0.017 | 0.858 | 0.011 | 1.259 | 0.043 | 0.520 | 0.014 | 7.624
CS-CNN | 0.007 | 0.017 | 0.426 | 0.012 | 0.609 | 0.093 | 0.513 | 0.014 | 4.425
E2E-LM | 0.049 | 0.044 | 1.116 | 0.030 | 1.649 | 0.168 | 0.546 | 0.011 | 7.388
E2E-MLP | 0.044 | 0.026 | 1.688 | 0.016 | 2.657 | 0.073 | 0.552 | 0.008 | 7.124
E2E-LSTM | 0.060 | 0.023 | 2.604 | 0.013 | 4.448 | 0.017 | 0.578 | 0.017 | 2.840
E2E-CNN | 0.023 | 0.024 | 0.931 | 0.017 | 1.365 | 0.084 | 0.531 | 0.046 | 3.881
Other Objective Functions
LSTM-GMVP | 0.001 | 0.011 | 0.047 | 0.008 | 0.067 | 0.060 | 0.504 | -0.004 | 0.170
LSTM-MVPλ=10 | 0.064 | 0.317 | 0.349 | 0.207 | 0.534 | 0.878 | 0.518 | 0.342 | 4.209
LSTM-MVPλ=20 | 0.179 | 0.169 | 1.055 | 0.115 | 1.555 | 0.380 | 0.548 | 0.195 | 3.910
LSTM-MVPλ=30 | 0.168 | 0.116 | 1.394 | 0.076 | 2.149 | 0.187 | 0.552 | 0.060 | 3.483
Constraints
LSTM-MSRP-LONG | 0.368 | 0.197 | 1.691 | 0.125 | 2.666 | 0.253 | 0.566 | 0.767 | 5.309
LSTM-MSRP-LEV | 0.321 | 0.112 | 2.540 | 0.068 | 4.203 | 0.151 | 0.576 | 0.132 | 13.29
LSTM-MSRP-CAR | 0.032 | 0.056 | 0.588 | 0.039 | 0.844 | 0.167 | 0.520 | -0.011 | 8.638
LSTM-MSRP-MAX | 0.057 | 0.021 | 2.683 | 0.012 | 4.459 | 0.026 | 0.578 | 0.021 | 2.671
LSTM-MSRP-LEV-CAR-MAX | 0.121 | 0.140 | 0.885 | 0.093 | 1.340 | 0.296 | 0.518 | 0.029 | 18.49

Table 3: Experimental results for the WRDS dataset with C = 0 bp.

| E(R) | Std(R) | Sharpe | DD(R) | Sortino | MDD | $\frac{\text{\% of}}{\text{+Ret}}$ | Beta | Turnover
---|---|---|---|---|---|---|---|---|---
Baselines
S&P 500 | 0.061 | 0.196 | 0.402 | 0.140 | 0.563 | 0.568 | 0.541 | 1.000 | 0.000
EWP | 0.129 | 0.212 | 0.678 | 0.148 | 0.968 | 0.549 | 0.546 | 1.000 | 0.000
MD | 0.414 | 0.239 | 1.569 | 0.141 | 2.649 | 0.531 | 0.546 | 0.599 | 2.347
DWP | 0.116 | 0.209 | 0.630 | 0.147 | 0.895 | 0.555 | 0.545 | 1.000 | 0.009
GMVP | 0.068 | 0.081 | 0.849 | 0.060 | 1.156 | 0.421 | 0.562 | 0.257 | 1.577
MSRP
CS-SAMPLE | -0.043 | 0.038 | -1.136 | 0.029 | -1.505 | 0.593 | 0.473 | 0.037 | 1.708
CS-LM | -0.060 | 0.015 | -4.082 | 0.013 | -4.777 | 0.719 | 0.382 | 0.008 | 9.160
CS-MLP | -0.058 | 0.027 | -2.175 | 0.021 | -2.760 | 0.710 | 0.440 | 0.036 | 9.319
CS-LSTM | -0.040 | 0.017 | -2.439 | 0.013 | -3.065 | 0.573 | 0.424 | 0.014 | 7.624
CS-CNN | -0.025 | 0.017 | -1.486 | 0.013 | -1.935 | 0.407 | 0.461 | 0.014 | 4.425
E2E-LM | -0.005 | 0.044 | -0.090 | 0.031 | -0.127 | 0.257 | 0.487 | 0.011 | 7.388
E2E-MLP | -0.008 | 0.026 | -0.295 | 0.018 | -0.424 | 0.223 | 0.477 | 0.007 | 7.124
E2E-LSTM | 0.039 | 0.023 | 1.693 | 0.014 | 2.743 | 0.019 | 0.548 | 0.017 | 2.840
E2E-CNN | -0.006 | 0.024 | -0.217 | 0.017 | -0.302 | 0.261 | 0.487 | 0.046 | 3.881
Other Objective Functions
LSTM-GMVP | -0.001 | 0.011 | -0.105 | 0.008 | -0.147 | 0.070 | 0.499 | -0.004 | 0.170
LSTM-MVPλ=10 | 0.032 | 0.317 | 0.253 | 0.208 | 0.385 | 0.917 | 0.512 | 0.342 | 4.209
LSTM-MVPλ=20 | 0.146 | 0.170 | 0.888 | 0.116 | 1.296 | 0.422 | 0.542 | 0.195 | 3.910
LSTM-MVPλ=30 | 0.139 | 0.117 | 1.177 | 0.077 | 1.792 | 0.238 | 0.542 | 0.060 | 3.483
Constraints
LSTM-MSRP-LONG | 0.317 | 0.199 | 1.498 | 0.126 | 2.337 | 0.254 | 0.558 | 0.767 | 5.309
LSTM-MSRP-LEV | 0.200 | 0.113 | 1.679 | 0.071 | 2.647 | 0.169 | 0.548 | 0.131 | 13.29
LSTM-MSRP-CAR | -0.030 | 0.056 | -0.529 | 0.041 | -0.722 | 0.564 | 0.480 | -0.011 | 8.638
LSTM-MSRP-MAX | 0.037 | 0.021 | 1.738 | 0.013 | 2.741 | 0.030 | 0.551 | 0.021 | 2.671
LSTM-MSRP-LEV-CAR-MAX | -0.019 | 0.140 | -0.063 | 0.097 | -0.091 | 0.618 | 0.490 | 0.029 | 18.49

Table 4: Experimental results for the WRDS dataset with C = 2 bps.

To assess the robustness of our results, we look at the performance of
E2E-LSTM across the entire testing period. Figure 5 shows the cumulative
returns, rolling Sharpe Ratio, rolling beta against S&P 500 and a drawdown
plot for E2E-LSTM ($C=0$ bps). In general, we obtain decent monthly returns
without consistent losses, and the rolling Sharpe Ratio is mostly above zero.
Since our method allows for short selling, the rolling beta is always close to
zero indicating a market-neural strategy; therefore, our model is agnostic to
market movements. This observation can be also verified in the drawdown plot –
we observe that drawdown is well controlled during the 2008 crisis and 2020
market crash due to the pandemic, demonstrating robustness of the approach
even throughout turbulent periods.

Figure 5: Trading performance for E2E-LSTM. Top: Cumulative returns; Top
Middle: Rolling Sharpe ratio over time; Middle Bottom: Rolling beta over time;
Bottom: Drawdown plot.

##  5 Conclusion

In this work we propose a universal end-to-end framework for portfolio
optimization where portfolio weights are directly computed via a deep learning
architecture.

Our proposed pipeline bypasses the traditional two-stage procedure of first
forecasting returns and then computing portfolio weights based on the inverse
of the estimated covariance matrix. It thus entirely avoids the challenging
problem of estimating the asset covariance matrix. Within our deep learning
framework we have the flexibility of optimizing various objective functions
including the Sharpe Ratio, mean-variance trade-off and others. Furthermore,
several neural network layers are designed to impose constraints on the
portfolio weights, such as cardinality, maximum position and leverage. These
constraints are popular in practice but suffer from computational complexity
when using classical approaches as commonly discussed in the literature. We
incorporate these constraints into the objective functions through special
neural network layers and optimize the entire framework via gradient ascent.

We test the efficacy of our methods on two data sets. The first one is a
synthetic data set where we compare the estimated weights from the classical
predictive approach and our end-to-end framework against the optimal weights
which are known from the data generating process. The results showcase the
accuracy of our method, and demonstrates that the weights derived by our
approach are “closer” to the optimal weights. The second data set consists of
hundreds of stocks from WRDS, with a testing period of more than 20 years. Our
model delivers superior results in terms of Sharpe Ratio, and we demonstrate
how to use constraints to achieve the desired portfolio, for example, using
leverage to increase returns or imposing maximum position for individual
instrument to prevent concentration of risk.

In subsequent continuations of this work, we aim to integrate transaction
costs into the objective functions in order for the turnover of a portfolio to
be controlled during training process; we expect the resulting method to be
able to tolerate much higher cost rates. This idea is closely related to
Reinforcement Learning (RL) [32], and future extensions of our framework will
incorporate insights from the RL literature.

## Acknowledgements

The authors would like to thank Oxford-Man Institute of Quantitative Finance
for computation support and data access.

## References

  * [1] Michel Ballings, Dirk Van den Poel, Nathalie Hespeels, and Ruben Gryp.  Evaluating multiple classifiers for stock price direction prediction.  Expert systems with Applications, 42(20):7046–7056, 2015.
  * [2] Mathieu Blondel, Olivier Teboul, Quentin Berthet, and Josip Djolonga.  Fast differentiable sorting and ranking.  In International Conference on Machine Learning, pages 950–959. PMLR, 2020.
  * [3] Benjamin Bruder and Thierry Roncalli.  Managing risk exposures using the risk budgeting approach.  Available at SSRN 2009778, 2012.
  * [4] Andrew Butler and Roy Kwon.  Integrating prediction in mean-variance portfolio optimization.  Available at SSRN 3788875, 2021.
  * [5] T-J Chang, Nigel Meade, John E Beasley, and Yazid M Sharaiha.  Heuristics for cardinality constrained portfolio optimisation.  Computers & Operations Research, 27(13):1271–1302, 2000.
  * [6] Marco Cuturi, Olivier Teboul, and Jean-Philippe Vert.  Differentiable ranking and sorting using optimal transport.  Advances in Neural Information Processing Systems, 32, 2019.
  * [7] Mirco Fabbri and Gianluca Moro.  Dow jones trading with deep learning: The unreasonable effectiveness of recurrent neural networks.  In Data, pages 142–153, 2018.
  * [8] Thomas Fischer and Christopher Krauss.  Deep learning with long short-term memory networks for financial market predictions.  European Journal of Operational Research, 270(2):654–669, 2018\.
  * [9] Marguerite Frank, Philip Wolfe, et al.  An algorithm for quadratic programming.  Naval research logistics quarterly, 3(1-2):95–110, 1956.
  * [10] Jianjun Gao and Duan Li.  Optimal cardinality constrained portfolio selection.  Operations research, 61(3):745–761, 2013.
  * [11] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep learning.  MIT press, 2016.
  * [12] Aditya Grover, Eric Wang, Aaron Zweig, and Stefano Ermon.  Stochastic optimization of sorting networks via continuous relaxations.  In International Conference on Learning Representations, 2018.
  * [13] John H Holland.  Genetic algorithms.  Scientific american, 267(1):66–73, 1992.
  * [14] Yan Jin, Rong Qu, and Jason AD Atkin.  Constrained portfolio optimisation: The state-of-the-art markowitz models.  ICORES, 2:388–395, 2016.
  * [15] Philippe Jorion.  Bayes-stein estimation for portfolio analysis.  Journal of Financial and Quantitative analysis, 21(3):279–292, 1986\.
  * [16] Raymond Kan and Guofu Zhou.  Optimal portfolio choice with parameter uncertainty.  Journal of Financial and Quantitative Analysis, 42(3):621–656, 2007\.
  * [17] Yves-Laurent Kom Samo and Alexander Vervuurt.  Stochastic Portfolio Theory: A machine learning perspective.  In Proceedings of the 32nd Conference on Uncertainty in Artificial Intelligence, 2016.
  * [18] Tze Leung Lai, Haipeng Xing, and Zehao Chen.  Mean–variance portfolio optimization when means and covariances are unknown.  The Annals of Applied Statistics, 5(2A):798–823, 2011.
  * [19] Eugene L Lawler and David E Wood.  Branch-and-bound methods: A survey.  Operations research, 14(4):699–719, 1966.
  * [20] Olivier Ledoit and Michael Wolf.  Improved estimation of the covariance matrix of stock returns with an application to portfolio selection.  Journal of empirical finance, 10(5):603–621, 2003.
  * [21] Olivier Ledoit and Michael Wolf.  Honey, I shrunk the sample covariance matrix.  The Journal of Portfolio Management, 30(4):110–119, 2004.
  * [22] Bryan Lim, Stefan Zohren, and Stephen Roberts.  Enhancing time-series momentum strategies using deep neural networks.  The Journal of Financial Data Science, 1(4):19–38, 2019.
  * [23] Harry Markowitz.  Portfolio selection.  The journal of finance, 7(1):77–91, 1952.
  * [24] Richard O Michaud.  The markowitz optimization enigma: Is ‘optimized’optimal?  Financial analysts journal, 45(1):31–42, 1989.
  * [25] Wlodzimierz Ogryczak and Arie Tamir.  Minimizing the sum of the k largest functions in linear time.  Information Processing Letters, 85(3):117–122, 2003.
  * [26] Ester Pantaleo, Michele Tumminello, Fabrizio Lillo, and Rosario N Mantegna.  When do improved covariance matrix estimators enhance portfolio optimization? An empirical comparative study of nine estimators.  Quantitative Finance, 11(7):1067–1080, 2011.
  * [27] L’uboš Pástor and Robert F Stambaugh.  Comparing asset pricing models: An investment perspective.  Journal of Financial Economics, 56(3):335–381, 2000.
  * [28] Daniel Poh, Bryan Lim, Stefan Zohren, and Stephen Roberts.  Building cross-sectional systematic strategies by learning to rank.  The Journal of Financial Data Science, 3(2):70–86, 2021.
  * [29] Eberhard Schöneburg.  Stock price prediction using neural networks: A project report.  Neurocomputing, 2(1):17–27, 1990.
  * [30] William F Sharpe.  The sharpe ratio.  Journal of portfolio management, 21(1):49–58, 1994.
  * [31] Dong X Shaw, Shucheng Liu, and Leonid Kopman.  Lagrangian relaxation procedure for cardinality-constrained portfolio optimization.  Optimisation Methods & Software, 23(3):411–420, 2008.
  * [32] Richard S Sutton and Andrew G Barto.  Reinforcement learning: An introduction.  MIT press, 2018.
  * [33] Ludan Theron and Gary Van Vuuren.  The maximum diversification investment strategy: A portfolio performance comparison.  Cogent Economics & Finance, 6(1):1427533, 2018.
  * [34] Sinem Uysal, Xiaoyue Li, and John M Mulvey.  End-to-end risk budgeting portfolio optimization with neural networks.  Available at SSRN 3883614, 2021.
  * [35] Philip Wolfe.  The simplex method for quadratic programming.  Econometrica: Journal of the Econometric Society, pages 382–398, 1959.
  * [36] Wharton Research Data Services.  https://wrds-www.wharton.upenn.edu/, Accessed: 2021-09-09.
  * [37] Zihao Zhang, Stefan Zohren, and Stephen Roberts.  Deep learning for portfolio optimization.  The Journal of Financial Data Science, 2(4):8–20, 2020.


# An invariant measure of chiral quantum transport

K. Ziegler
Institut für Physik, Universität Augsburg
D-86135 Augsburg, Germany

Abstract:

We study the invariant measure of the transport correlator for a chiral
Hamiltonian and analyze its properties. The Jacobian of the invariant measure
is a function of random phases. Then we distinguish the invariant measure
before and after the phase integration. In the former case we found quantum
diffusion of fermions and a uniform zero mode that is associated with particle
conservation. After the phase integration the transport correlator reveals two
types of evolution processes, namely classical diffusion and back-folded
random walks. Which one dominates the other depends on the details of the
underlying chiral Hamiltonian and may lead either to classical diffusion or to
the suppression of diffusion.

##  1 Introduction

We consider the quantum evolution of a particle from the site ${\bf
r}^{\prime}$ to the site ${\bf r}$ during the time period $t$ in a particle
conserving system. Then the building block for quantum transport is provided
by the transition probability $|\langle{\bf r}|e^{-iHt}|{\bf
r}^{\prime}\rangle|^{2}$, which is obtained from the transition amplitude
$\langle{\bf r}|e^{-iHt}|{\bf r}^{\prime}\rangle$ over time $t$ for a system
defined by the Hamiltonian $H$. We can Laplace transform the amplitude with
respect to time and calculate the resulting transition probability, the
transport correlator, as

$\Big{|}\int_{0}^{\infty}\langle{\bf r}|e^{-iHt}|{\bf r}^{\prime}\rangle
e^{-\epsilon t}dt\Big{|}^{2}=|G_{{\bf r}{\bf r}^{\prime}}(\epsilon)|^{2},$ (1)

where $G(\epsilon)=(H-i\epsilon)^{-1}$ is the one-particle Green’s function.
Besides the spatial coordinates ${\bf r}$, the Hamiltonian can also depend on
additional quantum numbers, such as a spinor or band index. Quantum transport
in multiband systems has been a popular subject in recent years due to the
discovery of new materials, beginning with graphene [1], topological
insulators [2, 3, 4] and including more specialized models such as two-
dimensional models with random gauge field [5].

To describe realistic transport in a disordered material we must include
random scattering. This can be attributed to a random Hamiltonian $H$. A very
successful approach to random scattering is based on the random matrix theory
[6, 7], where $H$ is a symmetric or Hermitian matrix with independently and
identically distributed matrix elements $H_{ij}$ for $i\leq j$. Many
interesting properties (e.g., level repulsion statistics) can be calculated
from the invariant measure. This is based on the idea that the fluctuations of
the $N$ eigenvalues of, for instance, a symmetric $N\times N$ random matrix
$H$, which are invariant under the orthogonal transformation of the matrix,
are separated from the remaining $N(N-1)/2$ random degrees of freedom of the
matrix elements. Then the Jacobian of the latter provides the invariant
measure. This concept can also been applied to the more general class of
random matrices, which are defined on a lattice and describe hopping of
quantum particles on this lattice [8]. Instead of fixing the eigenvalues of
these matrices, one can also use a uniform approximation of the matrix that is
fixed by a variational approach (saddle-point approximation). Then the
fluctuations induced by the chiral symmetry are represented as the invariant
measure.

The goal of the following analysis is to start from the disorder averaged
transport correlator $\langle|G_{{\bf r}{\bf
r}^{\prime}}(\epsilon)|^{2}\rangle_{d}$ and reduce the general disorder
average $\langle\ldots\rangle_{d}$ by the integration with respect to its
invariant measure. The justification of this reduction is that the long-range
behavior is associated with the invariant measure of the underlying chiral
property of a particle-hole symmetric Hamiltonian $H$. It will be shown that
the invariant measure provides the long-range properties even before
performing the phase integration.

The paper is organized as follows. After a brief discussion of the chiral
Hamiltonian $H$, the related Green’s functions, and the corresponding chiral
invariance in Sect. 2, we give a summary of the invariant measure, which was
originally derived in a previous paper, in Sect. 3. This is followed by an
analysis of the corresponding Jacobian and a discussion of the effects caused
by phase averaging. In Sect. 4 the different mappings, related to the
derivation of the invariant measure, are summarized and an example for a
simple two-band Hamiltonian is studied. Finally, we provide some conclusions
with open problems and a brief outlook in Sect. 5.

##  2 Model

The transport correlator in Eq. (1) is a fundamental quantity for the
description of quantum transport [9], consisting of products of matrix
elements of the advanced and the retarded one-particle Green’s function
$G(\pm\epsilon)=(H\pm i\epsilon)^{-1}$ in the form of $G_{{\bf r}{\bf
r}^{\prime}}(\epsilon)G_{{\bf r}^{\prime}{\bf
r}}(-\epsilon)=(H-i\epsilon)^{-1}_{{\bf r}{\bf
r}^{\prime}}(H+i\epsilon)^{-1}_{{\bf r}^{\prime}{\bf r}}$. In the following we
consider a multiband Hamiltonian $H$ on a $d$-dimensional lattice $\Lambda$
with $|\Lambda|<\infty$ sites and assume that it is related to its transposed
$H^{T}$ through the chiral relation $UH^{T}U^{\dagger}=-H$, where $U$ is a
spatially uniform unitary matrix. Then the chiral relation implies that the
trace and the determinant for the advanced/retarded Green’s functions obey
${\rm Tr}[G(-\epsilon)]=-{\rm Tr}[G(\epsilon)]$ and
$\det[G(-\epsilon)]=\det[-G(\epsilon)]$. Using the Pauli matrices
$\\{\sigma_{j}\\}_{j=1,2,3}$ and the $2\times 2$ unit matrix $\sigma_{0}$, a
two-band Hamiltonian can be expressed by an expansion in terms of these
matrices. An example is the Hamiltonian
$H=h_{1}\sigma_{1}+h_{2}\sigma_{2}+m\sigma_{3}$ with a random $m$ that is
independently and identically distributed on the lattice, while the matrices
$h_{1,2}$ act on the lattice. For symmetric matrices $h_{1,2}$ this
Hamiltonian satisfies the chiral relation for $U=\sigma_{2}$.

The pair of the advanced and retarded Green’s functions can also be written as
blocks of an advanced Green’s function when we replace the Hamiltonian $H$ by
the block-diagonal Hamiltonian ${\hat{H}}=diag(H,H^{T})$. Then the
corresponding advanced Green’s function reads
${\hat{G}}(\epsilon)=({\hat{H}}-i\epsilon)^{-1}$ and we get
$G(\epsilon)={\hat{G}}_{11}(\epsilon)$ and
$G(-\epsilon)=-U{\hat{G}}_{22}(\epsilon)U^{\dagger}$. For

${\hat{S}}=\pmatrix{0&sU\cr s^{\prime}U^{\dagger}&0\cr}$ (2)

with two independent parameters $s$ and $s^{\prime}$ we obtain for ${\hat{H}}$
the symmetry relation

$e^{\hat{S}}{\hat{H}}e^{\hat{S}}={\hat{H}},$ (3)

since ${\hat{H}}$ and ${\hat{S}}$ anticommute. It was shown that this symmetry
is associated with an invariant measure [8].

##  3 Invariant measure

In general, for a random chiral Hamiltonian there exists an invariant measure
of the disorder average. To demonstrate this we consider the example of the
previous section $H=h_{1}\sigma_{1}+h_{2}\sigma_{2}+m\sigma_{3}$ and assume a
random $m_{\bf r}$ at each lattice site ${\bf r}$ that obeys a Gaussian
distribution $\exp(-m_{\bf r}^{2}/g)dm_{\bf r}$. To control the strong
fluctuations of $|G_{{\bf R}{\bf R}^{\prime}}(\epsilon)|^{2}$ near the poles
of the Green’s functions, we transform the random field $m_{\bf r}$ to another
random field $Q_{\bf r}$ such that this mapping provides the relation

$\langle|G_{{\bf R}{\bf R}^{\prime}}(\epsilon)|^{2}\rangle_{d}=\langle Q_{\bf
R}Q_{{\bf R}^{\prime}}\rangle_{Q}.$ (4)

Thus, the average transport correlator is a correlation function of the random
field $Q_{\bf r}$ [10]. The mapping $m_{\bf r}\to Q_{\bf r}$ yields a
Jacobian, and the Gaussian weight $\exp(-m_{\bf r}^{2}/g)$ transforms to the
corresponding weight $P(\\{Q_{\bf r}\\})$ for the field $Q_{\bf r}$, where the
latter is correlated on the lattice. Separating the integration over $Q_{\bf
r}$ in two parts, namely one that leaves $P(\\{Q_{\bf r}\\})$ invariant and
one that does not. The former, together with the corresponding Jacobian $J$,
defines an invariant measure. More specific, it was shown that for a chiral
Hamiltonian the invariant measure is given by an integration over random
phases, where the Jacobian is $J=\det C$ with the $|\Lambda|\times|\Lambda|$
random-phase matrix $C$. The elements of this matrix are [8]

$C_{{\bf r}{\bf r}^{\prime}}=2\delta_{{\bf r}{\bf
r}^{\prime}}-\sum_{\mu,\mu^{\prime}=1,2}z_{{\bf r}\mu}h_{{\bf r}\mu,{\bf
r}^{\prime}\mu^{\prime}}\sum_{{\bf
r}^{\prime\prime}}\sum_{\mu^{\prime\prime}=1,2}h^{\dagger}_{{\bf
r}^{\prime}\mu^{\prime},{\bf r}^{\prime\prime}\mu^{\prime\prime}}z^{*}_{{\bf
r}^{\prime\prime}\mu^{\prime\prime}}\ {\rm with}\ z_{{\bf
r}\mu}=e^{i\alpha_{{\bf r}\mu}}\ (0\leq\alpha_{{\bf r}\mu}<2\pi)$ (5)

with the non-random $2|\Lambda|\times 2|\Lambda|$ matrix

$h={\bf 1}+2i\eta(\bar{H}-i\eta-i\epsilon)^{-1},$ (6)

which represents an effective Green’s function. $\eta$ ($0<\eta<\infty$) is an
effective scattering rate, usually obtained either from experimental
observations, from numerical simulations or from a self-consistent Born
approximation of the average one-particle Green’s function. It is proportional
to the average density of states, indicating that on average there are state
with non-zero density for $\eta>0$. For our approach it must be positive but
its actual value is not important for the subsequent discussion. Finally,
${\bar{H}}$ is the average Hamiltonian: ${\bar{H}}=\langle H\rangle_{d}$. It
is crucial to note that in the limit $\epsilon\to 0$ the matrix $h$ is unitary
(i.e., $hh^{\dagger}={\bf 1}$). $\epsilon>0$ is only necessary to avoid the
uniform zero mode with $\sum_{{\bf r}^{\prime}}C_{{\bf r}{\bf r}^{\prime}}=0$.
After removing it from the spectrum (cf. App. A), we can take the limit
$\epsilon\to 0$. The phase integration provides the relation

$\langle|G_{{\bf R}{\bf
R}^{\prime}}(\epsilon)|^{2}\rangle_{d}\sim\frac{\langle adj_{{\bf R}{\bf
R}^{\prime}}C\rangle_{\alpha}}{\langle\det C\rangle_{\alpha}},$ (7)

where the adjugate matrix is $adj_{{\bf R}{\bf R}^{\prime}}C\equiv
C^{-1}_{{\bf R}{\bf R}^{\prime}}\det C$, which is obtained from the
determinant by differentation as

$\partial_{W_{{\bf R}^{\prime}{\bf R}}}\det(C+W)|_{W=0}=\sum_{\pi}\delta_{{\bf
R}\pi({\bf R}^{\prime})}{\rm sgn}\pi\prod_{{\bf r}\neq{\bf R}^{\prime}}C_{{\bf
r}\pi({\bf r})}.$ (8)

The relation $\sim$ preserves the asymptotic properties on large distances
$|{\bf R}-{\bf R}^{\prime}|$. The phase average
$\langle\ldots\rangle_{\alpha}$ is taken with respect to a statistically
independent and uniform distribution of the phases $\\{\alpha_{{\bf r}\mu}\\}$
on the interval $[0,2\pi)$. We will discuss the phase average only in a second
step but treat the invariant measure in the following for a general
realization of the random phases. Averaging can be performed explicitly and
leads to a loop expansion of the determinant and the adjugate matrix [8].

The random-phase matrix $C$ has some remarkable properties, which will be
analyzed subsequently, revealing eventually a spatial scaling behavior of
$\langle|G_{{\bf R}{\bf R}^{\prime}}(\epsilon)|^{2}\rangle_{d}$ and quantum
diffusion of fermions. To this end we consider the random unitary matrix
$u=zhz^{*}$ and rewrite $C$ as

$C_{{\bf r}{\bf r}^{\prime}}={\rm Tr}_{2}[{\hat{C}}_{{\bf r}{\bf
r}^{\prime}}E_{2}]\ \ \ {\rm with}\ \ {\hat{C}}_{{\bf r}\mu,{\bf
r}^{\prime}\mu^{\prime}}=\delta_{{\bf r}{\bf
r}^{\prime}}\delta_{\mu\mu^{\prime}}-u_{{\bf r}\mu,{\bf
r}^{\prime}\mu^{\prime}}\sum_{{\bf
r}^{\prime\prime},\mu^{\prime\prime}}u^{\dagger}_{{\bf
r}^{\prime}\mu^{\prime},{\bf r}^{\prime\prime}\mu^{\prime\prime}},$ (9)

where $E_{2}=\sigma_{0}+\sigma_{1}$ and ${\rm Tr}_{2}$ is the trace with
respect to $2\times 2$ matrices for a Hamiltonian $H$ with two bands. A
generalization to a multiband Hamiltonian is straightforward but not
considered explicitly. Moreover, with the $2|\Lambda|\times 2|\Lambda|$ matrix
$E$, whose elements are 1, we can write $[u^{\dagger}E]_{{\bf
r}\mu}=\sum_{{\bf r}^{\prime\prime},\mu^{\prime\prime}}u^{\dagger}_{{\bf
r}\mu,{\bf r}^{\prime\prime}\mu^{\prime\prime}}$ and get

${\hat{C}}={\bf 1}-u[u^{\dagger}E]_{d}={\bf
1}-uu^{\dagger}+u(u^{\dagger}-[u^{\dagger}E]_{d}),$ (10)

where $[A]_{d}$ is the diagonal part of the general matrix $A$ and ${\bf
1}-uu^{\dagger}=4\epsilon\eta[{\bar{H}}^{2}+(\eta+\epsilon)^{2}]^{-1}$. Now we
split $u^{\dagger}$ into its diagonal part $u^{\prime}:=[u^{\dagger}]_{d}$ and
its off-diagonal part $v:=u^{\dagger}-u^{\prime}$ and insert this into
$u^{\dagger}-[u^{\dagger}E]_{d}$. This gives

$u^{\prime}+v-[u^{\prime}E]_{d}-[vE]_{d}=v-[vE]_{d},$ (11)

since $[u^{\prime}E]_{d}=u^{\prime}$. Thus, the diagonal part $u^{\prime}$
drops out and reveals that
$z^{*}(v-[vE]_{d})z=h^{\dagger}-[h^{\dagger}]_{d}-[vE]_{d}$ can be understood
as a hopping matrix $z^{*}vz$ with an additional random diagonal part
$-[vE]_{d}$. In contrast to the conventional disordered tight-binding model
[11, 12, 13, 14], the diagonal term also scales with the hopping matrix $v$,
since it consists of a sum of hopping terms. In other words, when $v$ changes
as $v\to L^{p}v$ under a change of the lattice constant $a\to La$, the entire
matrix $v-[vE]_{d}$ scales with $L^{p}$.

We return to $C$ in Eq. (9) and ${\hat{C}}$ in Eq. (10) and write

$C=\epsilon C^{\prime}+C^{\prime\prime}\ \ {\rm with}\ \ C^{\prime}=4\eta{\rm
Tr}_{2}[(z{\bar{H}}^{2}z^{*}+(\eta+\epsilon)^{2})^{-1}E_{2}],\ \
C^{\prime\prime}={\rm Tr}_{2}[u(v-[vE]_{d})E_{2}].$ (12)

The properties of the determinant and the adjugate matrix of $C$ are linked to
the spectral properties of ${\bf 1}-{u^{\dagger}}^{-1}[u^{\dagger}E]_{d}$ in
Eq. (10). To analyze this matrix we rewrite it as a quadratic form in $z$:

$z_{{\bf r}\mu}(\delta_{{\bf r}{\bf r}^{\prime}}\delta_{{\bf r}{\bf
r}^{\prime\prime}}\delta_{\mu\mu^{\prime\prime}}-\sum_{\mu^{\prime}}h^{\dagger-1}_{{\bf
r}\mu,{\bf r}^{\prime}\mu^{\prime}}h^{\dagger}_{{\bf
r}^{\prime}\mu^{\prime},{\bf r}^{\prime\prime}\mu^{\prime\prime}})z^{*}_{{\bf
r}^{\prime\prime}\mu^{\prime\prime}}$ $=z_{{\bf
r}\mu}\sum_{\mu^{\prime}}(\delta_{{\bf r}{\bf r}^{\prime}}\sum_{\bar{{\bf
r}}}h^{\dagger-1}_{{\bf r}\mu,\bar{{\bf r}}\mu^{\prime}}h^{\dagger}_{\bar{{\bf
r}}\mu^{\prime},{\bf r}^{\prime\prime}\mu^{\prime\prime}}-h^{\dagger-1}_{{\bf
r}\mu,{\bf r}^{\prime}\mu^{\prime}}h^{\dagger}_{{\bf
r}^{\prime}\mu^{\prime},{\bf r}^{\prime\prime}\mu^{\prime\prime}})z^{*}_{{\bf
r}^{\prime\prime}\mu^{\prime\prime}},$ (13)

which yields

$C^{\prime\prime}_{{\bf r}{\bf r}^{\prime}}=2\sum_{\bar{{\bf r}}}(\delta_{{\bf
r}{\bf r}^{\prime}}-\delta_{\bar{{\bf r}}{\bf r}^{\prime}})K_{{\bf r}\bar{{\bf
r}}}=2(\delta_{{\bf r}{\bf r}^{\prime}}\sum_{\bar{{\bf r}}}K_{{\bf r}\bar{{\bf
r}}}-K_{{\bf r}{\bf r}^{\prime}}),$ (14)

where

$K_{{\bf r}{\bf
r}^{\prime}}=\frac{1}{2}\sum_{\mu,\mu^{\prime},\mu^{\prime\prime}}\sum_{{\bf
r}^{\prime\prime}}u^{\dagger-1}_{{\bf r}\mu,{\bf
r}^{\prime}\mu^{\prime}}u^{\dagger}_{{\bf r}^{\prime}\mu^{\prime},{\bf
r}^{\prime\prime}\mu^{\prime\prime}}\ \ {\rm with}\ \ \sum_{\bar{{\bf
r}}}K_{{\bf r}\bar{{\bf r}}}=1.$ (15)

This indicates that $K_{{\bf r}{\bf r}^{\prime}}$ is formally a random
transition amplitude for ${\bf r}^{\prime}\to{\bf r}$ with particle
conservation due to $\sum_{\bar{{\bf r}}}K_{{\bf r}\bar{{\bf r}}}=1$. Then
$C^{\prime\prime}$, defined in Eq. (14), can be understood as a model for
quantum diffusion. In contrast to classical diffusion its transition rates are
complex numbers rather than probabilities. On the other hand, the average
amplitude is a positive number, given by the isotropic mean

$\langle K_{{\bf r}{\bf
r}^{\prime}}\rangle_{\alpha}=\frac{1}{2}\sum_{\mu,\mu^{\prime}}h^{\dagger-1}_{{\bf
r}\mu,{\bf r}^{\prime}\mu^{\prime}}h^{\dagger}_{{\bf
r}^{\prime}\mu^{\prime},{\bf r}\mu}=\frac{1}{2}\sum_{\mu,\mu^{\prime}}|h_{{\bf
r}\mu,{\bf r}^{\prime}\mu^{\prime}}|^{2}$ (16)

for $h^{\dagger-1}=h$ in the limit $\epsilon\to 0$. This implies

$\langle C^{\prime\prime}_{{\bf r}{\bf
r}^{\prime}}\rangle_{\alpha}=2(\delta_{{\bf r}{\bf r}^{\prime}}-\langle
K_{{\bf r}{\bf r}^{\prime}}\rangle_{\alpha})=2(\delta_{{\bf r}{\bf
r}^{\prime}}-\sum_{\mu,\mu^{\prime}}|h_{{\bf r}\mu,{\bf
r}^{\prime}\mu^{\prime}}|^{2}).$ (17)

These results can be interpreted as follows. Before phase average the
expansion of the determinant $\det C$ or the adjugate matrix $adj_{{\bf R}{\bf
R}^{\prime}}C$ provide quantum walks on the lattice in the form of loops,
where a site is only visited once. This reflects the fact that the determinant
represents fermions, which obey Pauli’s exclusion principle. Thus, the
invariant measure describes quantum diffusion of fermions.

After phase average the physics is quite different though, since the average
connects lattice sites ${\bf r}$ and ${\bf r}^{\prime}$, visited by the
quantum walk, with the transition amplitude $h^{\dagger}_{{\bf r}{\bf
r}^{\prime}}$. The latter decays on large distances exponentially, while on
short distances it depends strongly on the specific structure of the
Hamiltonian ${\bar{H}}$. This is as if we connect the quantum walk between the
two lattice sites by an elastic rubber band with an exponential strength for
large stretches. In other words, phase averaging is subject to quantum
interference effects, which is sensitive to spatial directions and can cause
an anisotropic evolution. The average transition amplitude in Eq. (16) is
isotropic, since $\langle K_{{\bf r}{\bf r}^{\prime}}\rangle_{\alpha}$ is
equal in all lattice directions. Therefore, the length of the rubber band is
the same in all directions. The situation is different though for two
consecutive steps ${\bf r}\to{\bf r}^{\prime}\to{\bf r}^{\prime\prime}$:

$\langle K_{{\bf r}{\bf r}^{\prime}}K_{{\bf r}^{\prime}{\bf
r}^{\prime\prime}}\rangle_{\alpha}=\langle K_{{\bf r}{\bf
r}^{\prime}}\rangle_{\alpha}\langle K_{{\bf r}^{\prime}{\bf
r}^{\prime\prime}}\rangle_{\alpha}+\frac{1}{4}\sum_{\mu,\mu^{\prime},\mu^{\prime\prime},\mu^{\prime\prime\prime}}h^{\dagger-1}_{{\bf
r}\mu,{\bf r}^{\prime}\mu^{\prime}}h^{\dagger}_{{\bf
r}^{\prime}\mu^{\prime},{\bf
r}^{\prime}\mu^{\prime\prime\prime}}h^{\dagger-1}_{{\bf
r}^{\prime}\mu^{\prime\prime\prime},{\bf
r}^{\prime\prime}\mu^{\prime\prime}}h^{\dagger}_{{\bf
r}^{\prime\prime}\mu^{\prime\prime},{\bf r}\mu}.$ (18)

While the first term on the right-hand side is isotropic, the second term
depends on the site ${\bf r}^{\prime\prime}$ relative to ${\bf r}$ and ${\bf
r}^{\prime}$. This is visualized in Fig. 1c) and calculated for a special
example in Sect. 4. While for the first term the length of the rubber band
does not change when going from ${\bf r}^{\prime}$ to any of the available
${\bf r}^{\prime\prime}$ with fixed length $|{\bf r}^{\prime\prime}-{\bf
r}^{\prime}|$, the rubber band has different lengths for the second term. The
latter favors a return of the quantum walk closer to ${\bf r}$ rather than to
move away from ${\bf r}$. This picture can be generalized to more steps,
connected by rubber bands between pairs of sites along the quantum walk. In
general, there is one specific pairing of sites, where the lengths of the
rubber bands are not affected for all possible conformations of the random
walks, namely when the sites of the quantum walk are connected in each step.
This corresponds to the factorization of the phase average

$\langle K_{{\bf r}{\bf r}_{1}}\rangle_{\alpha}\cdots\langle K_{{\bf
r}_{n-1}{\bf r}_{n}}\rangle_{\alpha}\ ,$

which is visualized in the lower graph of Fig. 1b). It represents a classical
random walk with transition probabilities $\langle K_{{\bf r}{\bf
r}^{\prime}}\rangle_{\alpha}\geq 0$ and $\sum_{\bar{{\bf r}}}\langle K_{{\bf
r}\bar{{\bf r}}}\rangle_{\alpha}=1$. The co-existence of the classical random
walk and the folded walks with rubber bands indicates a competition between
those walks. To understand this point better we rewrite the right-hand side of
Eq. (7) as

$\frac{\langle adj_{{\bf R}{\bf R}^{\prime}}C\rangle_{\alpha}}{\langle\det
C\rangle_{\alpha}}=\frac{adj_{{\bf R}{\bf R}^{\prime}}\langle
C\rangle_{\alpha}+\sum_{j}a_{j;{\bf R}{\bf R}^{\prime}}}{\det\langle
C\rangle_{\alpha}+\sum_{j}d_{j}},$ (19)

where the two sums are finite on a finite lattice and are created by a
systematic expansion in terms of truncated correlations (linked cluster
expansion) as

$\langle\prod_{\bf r}C_{{\bf r}\pi({\bf r})}\rangle=\prod_{\bf r}\langle
C_{{\bf r}\pi({\bf r})}\rangle+\sum_{{\bf r},{\bf r}^{\prime}}\left[\langle
C_{{\bf r}\pi({\bf r})}C_{{\bf r}^{\prime}\pi({\bf
r}^{\prime})}\rangle_{\alpha}-\langle C_{{\bf r}\pi({\bf
r})}\rangle_{\alpha}\langle C_{{\bf r}^{\prime}\pi({\bf
r}^{\prime})}\rangle_{\alpha}\right]\prod_{\bar{{\bf r}}\neq{\bf r},{\bf
r}^{\prime}}\langle C_{\bar{{\bf r}}\pi(\bar{{\bf
r}})}\rangle_{\alpha}+\cdots$ (20)

for the determinant and the adjugate matrix. Then it depends on the details of
the Hamiltonian ${\bar{H}}$ and on the scattering rate $\eta$, which term
dominates the others. On an elementary level, we can compare the first and the
second term on the right-hand side of Eq. (18), as discussed in Sect. 4. For
large distances the classical random walk would be important. On the other
hand, for a given number of steps the sum over all realizations of the
classical random walk may have a much smaller weight than the sum over all
compactly folded walks, which are created by the rubber bands. Thus, we must
compare the weights of all walks with a given number of steps and extract
those with the highest weight as dominant.

Figure 1: a) A realization of a random walk from site ${\bf r}_{1}$ to site
${\bf r}_{5}$ before phase average, consisting of four steps: $K_{{\bf
r}_{1}{\bf r}_{2}}K_{{\bf r}_{2}{\bf r}_{3}}K_{{\bf r}_{3}{\bf r}_{4}}K_{{\bf
r}_{4}{\bf r}_{5}}$. According to Eq. (15), the black links represent
$u^{\dagger-1}$, while the red arrows represent the factors $\sum_{{\bf
r}^{\prime\prime},\mu^{\prime\prime}}u^{\dagger}_{{\bf r}\mu,{\bf
r}^{\prime\prime}\mu^{\prime\prime}}$. b) The same random walk after phase
average: Two typical contributions are shown for $\langle K_{{\bf r}_{1}{\bf
r}_{2}}K_{{\bf r}_{2}{\bf r}_{3}}K_{{\bf r}_{3}{\bf r}_{4}}K_{{\bf r}_{4}{\bf
r}_{5}}\rangle_{\alpha}$, where the red connection is $h^{\dagger}$ and
corresponds to a “rubber band”. The lower graph represents a classical random
walk $\langle K_{{\bf r}_{1}{\bf r}_{2}}\rangle_{\alpha}\langle K_{{\bf
r}_{2}{\bf r}_{3}}\rangle_{\alpha}\langle K_{{\bf r}_{3}{\bf
r}_{4}}\rangle_{\alpha}\langle K_{{\bf r}_{4}{\bf r}_{5}}\rangle_{\alpha}$. c)
This is the two step walk $\langle K_{{\bf r}{\bf r}^{\prime}}K_{{\bf
r}^{\prime}{\bf r}^{\prime\prime}}\rangle_{\alpha}-\langle K_{{\bf r}{\bf
r}^{\prime}}\rangle_{\alpha}\langle K_{{\bf r}^{\prime}{\bf
r}^{\prime\prime}}\rangle_{\alpha}$ of Eq. (18).

##  4 Discussion of the Results

This work on the invariant measure of the transport properties in chiral
systems is based on the random phase matrix $C=\epsilon
C^{\prime}+C^{\prime\prime}$ of Eq. (12), whose properties were studied in the
previous section. $C$ is connected with the original random Hamiltonian $H$
through several mappings between different random matrices. First, there is
the mapping between Hamiltonians

$H\to z{\bar{H}}z^{*}-i\eta$ (21)

with the phenomenological scattering rate $\eta>0$. This is accompanied by the
mapping of the Green’s function $(H-i\epsilon)^{-1}\to u$ with the effective
Green’s function of the Hamiltonian $z{\bar{H}}z^{*}-i\eta$

$u=zhz^{*}={\bf 1}+2i\eta(z{\bar{H}}z^{*}-i\eta-i\epsilon)^{-1},$ (22)

from which we eventually get $C$ in Eq. (12) and the transition amplitude $K$
in Eq. (15). $u$ is unitary in the limit $\epsilon\to 0$, which implies the
conservative transfer $\sum_{{\bf r}^{\prime}}K_{{\bf r}{\bf r}^{\prime}}=1$
in Eq. (14). The eigenvalues of $u$ are identical with those of $h$ in Eq.
(6), which is not random. The corresponding eigenvectors $\Psi_{{\bf k},{\bf
r}\mu}$ of $u$ are random though, given by the plane wave $\psi_{{\bf k},{\bf
r}}=\exp(i{\bf k}\cdot{\bf r})/\sqrt{|\Lambda|}$ with a random phase factor as
$\Psi_{{\bf k},{\bf r}\mu}=z_{{\bf r}\mu}\psi_{{\bf k},{\bf r}}$.

For strong scattering ($\eta\sim\infty$) we get $u\sim-{\bf 1}$, whereas for
weak scattering ($\eta\sim 0$) we obtain $u\sim{\bf 1}$. In both limits the
transition matrix becomes diagonal with $K={\bf 1}$, implying that the quantum
diffusion disappears. All these results are valid before phase averaging. The
transport correlator, defined through the relation in Eq. (7), requires phase
averaging of $\det C$ and the related adjugate matrix. Although the average
can be performed easily, the result leads to complex expressions, which can be
associated with classical random walks and folded random walks with rubber
bands. On a more elementary level we get for $u$ the diagonal mean $\langle
u_{{\bf r}\mu,{\bf r}^{\prime}\mu^{\prime}}\rangle_{\alpha}=\delta_{{\bf
r}{\bf r}^{\prime}}\delta_{\mu\mu^{\prime}}h_{{\bf r}\mu,{\bf r}\mu}$ and a
vanishing variance, while the mean $\langle C^{\prime\prime}\rangle_{\alpha}$
in Eq. (17) reveals classical diffusion. The fluctuations of
$C^{\prime\prime}$ are characterized by the variance

$\langle{C^{\prime\prime}_{{\bf r}{\bf
r}^{\prime}}}^{2}\rangle_{\alpha}-\langle C^{\prime\prime}_{{\bf r}{\bf
r}^{\prime}}\rangle_{\alpha}^{2}=\sum_{\mu}\left(\sum_{\mu^{\prime}}|h_{{\bf
r}\mu,{\bf r}^{\prime}\mu^{\prime}}|^{2}\right)^{2}.$ (23)

Then the ratio of the standard deviation and the mean value is

$R:=\frac{\sqrt{\langle C_{{\bf r}{\bf
r}^{\prime}}^{2}\rangle_{\alpha}-\langle C_{{\bf r}{\bf
r}^{\prime}}\rangle_{\alpha}^{2}}}{\langle C_{{\bf r}{\bf
r}^{\prime}}\rangle_{\alpha}}=\sqrt{1-2p_{1}p_{2}/(p_{1}+p_{2})^{2}}\ \ {\rm
with}\ \ p_{\mu}=\sum_{\mu^{\prime}}|h_{{\bf r}\mu,{\bf
r}^{\prime}\mu^{\prime}}|^{2},$ (24)

which is equal or less than 1. In other words, the fluctuations of $C$ around
the mean value $\langle C\rangle$ are restricted. Moreover, for $N$ bands with
equal $p_{\mu}\equiv{\bar{p}}$ we get $R=N^{-1/2}$, indicating a non-random
limit for a large number of bands $N\sim\infty$. In contrast, the fluctuations
of $|G_{{\bf r}{\bf r}^{\prime}}(\epsilon)|^{2}$ are unbounded for
$\epsilon\sim 0$ due to the proximity to the poles of the Green’s functions.
This demonstrates that the invariant measure is better controlled than the
original transport correlator.

In order to give an example for the anisotropic effect of the phase averaging,
we expand $h$ in powers of $1/\eta$:

$h=-{\bf
1}+i\frac{1}{\eta}{\bar{H}}+\frac{1}{\eta^{2}}{\bar{H}}^{2}+O(1/\eta^{3}).$
(25)

Then we choose a band-diagonal matrix ${\bar{H}}=\Delta\sigma_{3}$, where
$\Delta$ is nearest-neighbor hopping: $\Delta_{{\bf r}-{\bf r}^{\prime}}=1$
for ${\bf r}$, ${\bf r}^{\prime}$ nearest-neighbor sites and zero otherwise.
For the nearest-neighbor vectors ${\bf r}^{\prime}={\bf r}\pm{\bf e}_{1,2}$ on
a square lattice with lattice unit vectors ${\bf e}_{1,2}$ we obtain, after
neglecting terms of order $1/\eta^{3}$,

$h_{{\bf r},\mu}=\cases{-1+4/\eta^{2}&${\bf r}=0$\cr(-1)^{\mu}i/\eta&${\bf
r}=\pm{\bf e}_{1,2}$\cr 1/\eta^{2}&${\bf r}=\pm({\bf e}_{1}+{\bf
e}_{2}),\pm({\bf e}_{1}-{\bf e}_{2}),\pm 2{\bf e}_{1,2}$\cr 0&otherwise\cr}$
(26)

with $h_{{\bf r}\mu,{\bf r}^{\prime}\mu}\equiv h_{{\bf r}-{\bf
r}^{\prime},\mu}$. The weight $|h_{{\bf r}\mu}|$ decreases in powers of
$1/\eta$ as we increase the distance $|{\bf r}|$, with the maximal weight
$1-4/\eta^{2}$ for ${\bf r}=0$. This yields for the classical diffusion in two
steps with the special choice ${\bf r}^{\prime}={\bf r}+{\bf e}_{1}$ and ${\bf
r}^{\prime\prime}={\bf r}^{\prime}+{\bf R}={\bf r}+{\bf e}_{1}+{\bf R}$

$\langle K_{{\bf r}{\bf r}^{\prime}}\rangle_{\alpha}\langle K_{{\bf
r}^{\prime}{\bf
r}^{\prime\prime}}\rangle_{\alpha}=\frac{1}{4}\sum_{\mu,\mu^{\prime}}|h_{{\bf
r}-{\bf r}^{\prime},\mu}|^{2}\sum_{\mu}|h_{{\bf r}^{\prime}-{\bf
r}^{\prime\prime},\mu^{\prime}}|^{2}=\frac{1}{4}\sum_{\mu,\mu^{\prime}}|h_{-{\bf
e}_{1},\mu}|^{2}\sum_{\mu}|h_{-{\bf R},\mu^{\prime}}|^{2}$
$=\frac{1}{\eta^{4}}\cases{1&${\bf R}=\pm{\bf e}_{1,2}$\cr 1/\eta^{2}&${\bf
r}=\pm({\bf e}_{1}+{\bf e}_{2}),\pm({\bf e}_{1}-{\bf e}_{2}),\pm 2{\bf
e}_{1,2}$\cr}$ (27)

and for the corresponding two step term $\langle K_{{\bf r}{\bf
r}^{\prime}}K_{{\bf r}^{\prime}{\bf r}^{\prime\prime}}\rangle_{\alpha}-\langle
K_{{\bf r}{\bf r}^{\prime}}\rangle_{\alpha}\langle K_{{\bf r}^{\prime}{\bf
r}^{\prime\prime}}\rangle_{\alpha}$ of Eq. (18)

$\frac{1}{4}h_{-{\bf e}_{1},\mu}h^{\dagger}_{0,\mu}h_{-{\bf
R},\mu}h^{\dagger}_{{\bf R}+{\bf
e}_{1},\mu}=\frac{1-4/\eta^{2}}{4\eta^{2}}\cases{(-1+4/\eta^{2})&${\bf
R}=-{\bf e}_{1}$\cr 1/\eta^{2}&${\bf R}={\bf e}_{1},\pm{\bf e}_{2},-{\bf
e}_{1}\pm{\bf e}_{2},-2{\bf e}_{1}$\cr 0&otherwise\cr},$ (28)

which is real and does not depend on the band index $\mu$. It should be noted
that (i) ${\bf R}=-{\bf e}_{1}$ implies ${\bf r}^{\prime\prime}={\bf r}$,
which does not contribute to the determinant or to the adjugate matrix, and
(ii) the two step term vanishes for ${\bf R}={\bf e}_{1}\pm{\bf e}_{2}$,
$2{\bf e}_{1,2}$ and $-2{\bf e}_{2}$, indicating a complete suppression of the
transfer to these sites. This, in contrast to the non-vanishing value at ${\bf
R}=-{\bf e}_{1}\pm{\bf e}_{2}$ and ${\bf R}=-2{\bf e}_{1}$, represents the
anisotropy due to the rubber-band effect of the phase averaging.

The two examples in Fig. 1b) can be directly extended to larger random walks
by adding repeatedly either a link with a single-site loop (upper graph) or a
two-site loop $\langle K_{{\bf r}{\bf r}^{\prime}}\rangle_{\alpha}$ (lower
graph). With Eq. (26) the weight of the upper graph then changes by a factor
of order $1/\eta$, while the lower graph changes by a factor of order
$1/\eta^{2}$ for each additional link. This indicates that for this example
the weight of the exponentially decaying upper graph dominates the classical
random walk graph as the number of steps increases.

##  5 Conclusions

The analysis of the invariant measure has revealed that transport in systems
with a chiral Hamiltonian is linked to quantum diffusion of fermions, where
the evolution is described by the random transition amplitude $K_{{\bf r}{\bf
r}^{\prime}}$ with particle conservation $\sum_{{\bf r}^{\prime}}K_{{\bf
r}{\bf r}^{\prime}}=1$. This conservation law is reflected by a uniform zero
mode. Central for the derivation of the transition amplitude is the effective
Green’s function $u$ defined in Eq. (22). This can be seen as an approximation
of the original product of Green’s functions $G(\pm\epsilon)$. The
fluctuations of $u$ are bounded, in contrast to the unbounded fluctuations of
$G(\pm\epsilon)$. Through phase averaging, i.e., integration with respect to
the invariant measure, two types of competing walks are created, namely
classical random walks and back-folded random walks. The latter are formed as
if the walker is pulled back by rubber bands to its previous positions. The
competition of these contributions can lead to a suppression of classical
diffusion due to interference effects in phase averaging. Details of this
phenomenon for specific models are still open and should be addressed in the
future.

The limit $\eta\to 0$ is connected with a transition from the regime of
quantum diffusion to a non-diffusive regime. The latter is characterized by a
vanishing average density of states. Assuming a gap in the spectrum of
${\bar{H}}$, disorder in the form of, e.g., a random gap term $m$ in
$H=h_{1}\sigma_{1}+h_{2}\sigma_{2}+m\sigma_{3}$, can create localized (bound)
states at zero energy (midgap states), which are not visible in the average
density of states. The properties of the transition and the transport
properties in the regime with $\eta=0$, where the invariant measure approach
is not valid any more, should be studied in a separate project.

##  Appendix A Separation of the uniform zero mode

The role of the uniform zero mode can be discussed in term of the spectral
representation of the (normalized) correlator before averaging

$C^{-1}_{{\bf R}{\bf R}^{\prime}}=\sum_{k\geq 0}\frac{\psi^{*}_{k}({\bf
R})\psi_{k}({\bf R}^{\prime})}{\lambda_{k}}\ \ (\lambda_{0}=O(\epsilon))$

with eigenmodes $\psi_{k}$ and eigenvalues $\lambda_{k}$ of $C$. Since the
uniform zero mode $\psi_{0}=1/\sqrt{|\Lambda|}$ is not random and
$\lambda_{k}\neq 0$ for $k>0$, we separate this mode and split the sum into
$k=0$ and $k>0$:

$C^{-1}_{{\bf R}{\bf
R}^{\prime}}=\frac{1}{|\Lambda|}{\bar{K}}_{rec}+\sum_{k>0}\frac{\psi^{*}_{k}({\bf
R})\psi_{k}({\bf R}^{\prime})}{\lambda_{k}}\ ,\ \
\frac{1}{|\Lambda|}{\bar{K}}_{rec}=\frac{\psi^{*}_{0}({\bf R})\psi_{0}({\bf
R}^{\prime})}{\lambda_{0}}=\frac{1}{|\Lambda|\lambda_{0}},$ (29)

where the recurrent part with ${\bar{K}}_{rec}$ does not decay but vanishes
uniformly for a large system size $|\Lambda|$. The recurrent part provides the
$1/\epsilon$ behavior for the normalized $\sum_{{\bf
R}^{\prime}}{\hat{K}}_{{\bf R}{\bf R}^{\prime}}$, while the remaining sum is
independent of $\lambda_{0}$. After phase averaging we get for the right-hand
side of Eq. (7)

$\frac{\langle C^{-1}_{{\bf R}{\bf R}^{\prime}}\prod_{k^{\prime}\geq
0}\lambda_{k^{\prime}}\rangle}{\langle\prod_{k^{\prime}\geq
0}\lambda_{k^{\prime}}\rangle}=\frac{\lambda_{0}\langle C^{-1}_{{\bf R}{\bf
R}^{\prime}}\prod_{k^{\prime}>0}\lambda_{k^{\prime}}\rangle}{\lambda_{0}\langle\prod_{k^{\prime}>0}\lambda_{k^{\prime}}\rangle}=\frac{1}{|\Lambda|}{\bar{K}}_{rec}+\frac{1}{\langle\prod_{k^{\prime}>0}\lambda_{k^{\prime}}\rangle}\Big{\langle}\sum_{k>0}\frac{\psi^{*}_{k}({\bf
R})\psi_{k}({\bf
R}^{\prime})\prod_{k^{\prime}>0}\lambda_{k^{\prime}}}{\lambda_{k}}\Big{\rangle},$
(30)

since $\lambda_{0}$ does not depend on the random variables. This means that
the uniform zero mode gives an additive term to the correlator, which leads to
a divergence of $\sum_{{\bf R}^{\prime}}{\hat{K}}_{{\bf R}{\bf
R}^{\prime}}\sim 1/\lambda_{0}$ even for a finite system.

## References

  * [1] Novoselov K S, Geim A K, Morozov S V, Jiang D, Katsnelson M I, Grigorieva I V, Dubonos S V and Firsov A A 2005 Nature 438 197–200 ISSN 1476-4687 URL https://doi.org/10.1038/nature04233
  * [2] Qi X L, Hughes T L and Zhang S C 2008 Phys. Rev. B 78(19) 195424 URL https://link.aps.org/doi/10.1103/PhysRevB.78.195424
  * [3] Bernevig B and Hughes T 2013 Topological Insulators and Topological Superconductors (Princeton University Press) ISBN 9781400846733
  * [4] Haim A, Ilan R and Alicea J 2019 Phys. Rev. Lett. 123(4) 046801 URL https://link.aps.org/doi/10.1103/PhysRevLett.123.046801
  * [5] Li C A, Zhang S B, Budich J C and Trauzettel B 2022 Phys. Rev. B 106(8) L081410 URL https://link.aps.org/doi/10.1103/PhysRevB.106.L081410
  * [6] Dyson F J 1962 Journal of Mathematical Physics 3 140–156 (Preprint https://doi.org/10.1063/1.1703773) URL https://doi.org/10.1063/1.1703773
  * [7] Mehta M 2004 Random Matrices (Pure and applied mathematics no v. 142) (Elsevier/Academic Press) ISBN 9780120884094 URL https://books.google.de/books?id=1vKgAQAACAAJ
  * [8] Ziegler K 2015 Journal of Physics A: Mathematical and Theoretical 48 055102 ISSN 1751-8121 URL http://dx.doi.org/10.1088/1751-8113/48/5/055102
  * [9] Das Sarma S, Adam S, Hwang E H and Rossi E 2011 Rev. Mod. Phys. 83(2) 407–470 URL https://link.aps.org/doi/10.1103/RevModPhys.83.407
  * [10] Ziegler K 2009 Phys. Rev. B 79(19) 195424
URL https://link.aps.org/doi/10.1103/PhysRevB.79.195424

  * [11] Anderson P W 1958 Phys. Rev. 109(5) 1492–1505
URL https://link.aps.org/doi/10.1103/PhysRev.109.1492

  * [12] Abrahams E, Anderson P W, Licciardello D C and Ramakrishnan T V 1979 Phys. Rev. Lett. 42(10) 673–676 URL https://link.aps.org/doi/10.1103/PhysRevLett.42.673
  * [13] Wegner F 1979 Zeitschrift für Physik B Condensed Matter 35 207–210 ISSN 1431-584X URL https://doi.org/10.1007/BF01319839
  * [14] Vollhardt D and Wölfle P 1980 Phys. Rev. B 22(10) 4666–4679 URL https://link.aps.org/doi/10.1103/PhysRevB.22.4666