11institutetext: University College London, Department of Computer Science
London, United Kingdom
11email<EMAIL_ADDRESS>
# Portfolio Transformer for Attention-Based
Asset Allocation

Damian Kisiel  Denise Gorse

###### Abstract

Traditional approaches to financial asset allocation start with returns
forecasting followed by an optimization stage that decides the optimal asset
weights. Any errors made during the forecasting step reduce the accuracy of
the asset weightings, and hence the profitability of the overall portfolio.
The _Portfolio Transformer_ (PT) network, introduced here, circumvents the
need to predict asset returns and instead directly optimizes the Sharpe ratio,
a risk-adjusted performance metric widely used in practice. The PT is a novel
end-to-end portfolio optimization framework, inspired by the numerous
successes of attention mechanisms in natural language processing. With its
full encoder-decoder architecture, specialized time encoding layers, and
gating components, the PT has a high capacity to learn long-term dependencies
among portfolio assets and hence can adapt more quickly to changing market
conditions such as the COVID-19 pandemic. To demonstrate its robustness, the
PT is compared against other algorithms, including the current LSTM-based
state of the art, on three different datasets, with results showing that it
offers the best risk-adjusted performance.

###### Keywords:

Transformers Deep Learning Portfolio Optimization.

##  1 Introduction

Portfolio optimization algorithms aim to select the optimal weighting of
financial assets in a given portfolio as a means to maximize or minimize some
specific metric of interest. It is arguably the most important phase in the
entire investment lifecycle, without which investors would be exposed to
unacceptable levels of risk. Markowitz formally formulated this problem in
what is now known as _Modern Portfolio Theory_ (MPT) [14]. The risk-return
trade-off pioneered by Markowitz was very influential at the time and became
in effect a go-to tool for the vast majority of industry practitioners.
However, despite its rigorous theoretical foundations and wide popularity, the
MPT has significant shortcomings when applied in practice. One of these
limitations is the assumption that future investment returns of individual
assets can be predicted with a reasonable level of precision. This task has,
however, been shown to be extremely difficult due to the highly stochastic
nature of financial markets [15].

Recent advances in computational power and the wider availability of market
data have allowed machine learning architectures to be used for portfolio
optimization. For example, by using ensembles of gradient-boosted trees, one
can reduce estimation errors in the returns predictions of standard Markowitz-
style optimization [3], and the XGBoost model has been successfully used
within a meta-allocation framework to switch between different risk-based
strategies in order to achieve better risk-adjusted performance [12]. Most
recently, deep learning architectures have also started to play a major role
[7]. However, a substantial drawback of the above portfolio selection
methodologies is that they follow the classical two-step procedure in which
errors in the parameter estimations of the first step are translated into
inaccurate asset weightings in the second step.

Moody et al. [16] pioneered the contrasting idea of combining prediction and
performance optimization in a single step. This work was later extended by
that of Zhang et al. [24], who introduced an LSTM-based architecture that
showed significant performance improvements over classical asset allocation
techniques. In this paper, we introduce the _Portfolio Transformer_ (PT),
which combines prediction and optimization in a novel end-to-end deep learning
architecture based on an attention mechanism, directly outputting portfolio
weights that optimize the Sharpe ratio, a measure of risk-adjusted return
widely used in practice, under the specified transaction cost penalties.
Additionally, the PT makes use of specialized gating mechanisms to determine
the ideal level of non-linearity when optimizing each portfolio. Results
demonstrate that the Portfolio Transformer is able to outperform a number of
other methodologies, ranging from a classical optimization method to the
current LSTM-based state of the art [24], on three different datasets
encompassing ETFs, commodities, and stocks.

##  2 Background & Related Work

###  2.1 Long Short-Term Memory (LSTM)

There exists a large volume of literature applying recurrent neural networks,
such as simple RNNs [18] or Gated Recurrent Unit (GRU) networks [17], to
financial time series prediction problems. However, the main drawback of
standard recurrent neural networks, observed in multiple domains, including
finance, is the so-called ’vanishing gradient problem’ [8], whereby gradients
corresponding to long-term dependencies become very small, effectively
preventing the model from further training. LSTM networks [9] tackle this
problem by introducing gate mechanisms that allow gradients to flow unchanged;
through a process of filtering and summarizing they can ignore irrelevant past
information. The vanilla LSTM architecture proposed by Zhang et al. [24] for
portfolio optimization is used in this work as one of the benchmark
algorithms.

###  2.2 The Transformer Model

Despite their efficacy at learning time-localized patterns, LSTMs struggle to
capture meaningful dependencies when the length of a sequence is relatively
large. The Transformer architecture [19] was developed to address this issue
and has now been established as state of the art in most work in natural
language processing [20]. Several studies also show the successful application
of this architecture in the financial domain. For example, Wood et al. [21]
show that their Transformer model outperforms an LSTM network in time-series
momentum strategies, and Xu et al. [22] apply the model to portfolio policy
learning in a reinforcement learning setting (though to maximize cumulative
return rather than the industry-preferred Sharpe ratio).

At the heart of every Transformer architecture lies a mechanism called ’self-
attention’, which replaces recurrence and allows for simultaneous processing
of all sequence elements. The Portfolio Transformer implements ’scaled dot-
product attention’ [19], as given by

$Attention\left(Q,K,V\right)=softmax\left(\frac{QK^{T}}{\sqrt{d_{model}}}+M\right)V.$
(1)

The _value matrix_ $V\in\mathbb{R}^{\tau\times d_{v}}$ of equation (1) is
weighted by a set of ’scores’ obtained from the softmax operation, which
determines how much emphasis each time step from the _key matrix_
$K\in\mathbb{R}^{\tau\times d_{k}}$ should receive when encoding sequence
positions from the _query matrix_ $Q\in\mathbb{R}^{\tau\times d_{k}}$. The dot
product of $Q$ and $K$ is divided by the square root of the encoding dimension
($d_{model}$) to counteract problems associated with small gradients.
Additionally, the Portfolio Transformer implements masking via matrix
$M\in\mathbb{R}^{\tau\times\tau}$ in the first attention block of each decoder
layer to ensure it can only attend to preceding time steps and hence maintain
its autoregressive property. The operation defined by equation (1) is repeated
$h$ times in what is known as _multi-headed attention_ (MHA),

$MHA\left(Q,K,V\right)=Concatenate\left(head_{1},\cdots,head_{h}\right)W^{O},$
(2)

$head_{i}=Attention\left(QW_{i}^{Q},KW_{i}^{K},VW_{i}^{V}\right),$ (3)

which allows the model to extract information from multiple representation
subspaces, where each ${head}_{i}$ is implemented using its own set of learned
linear projection matrices $W_{i}^{Q}\in\mathbb{R}^{d_{model}\times d_{k}}$,
$W_{i}^{K}\in\mathbb{R}^{d_{model}\times d_{k}}$ and
$W_{i}^{V}\in\mathbb{R}^{d_{model}\times d_{v}}$. Outputs from all attention
heads are then concatenated and again linearly projected using a learned
parameter matrix $W^{O}\in\mathbb{R}^{{hd_{v}}\times d_{model}}$ to obtain
final values.

###  2.3 Gated Residual Network (GRN)

In the original Transformer implementation of [19] each attention layer is
followed by a simple feed-forward network. The Portfolio Transformer adopts a
more flexible approach and instead makes use of a _Gated Residual Network_
(GRN) [13], acting as a gating mechanism that determines the extent of non-
linear processing required for a particular portfolio, defined by

$GRN\left(z\right)=LayerNorm\left(z+GLU\left(g_{1}\right)\right),$ (4)

$g_{1}=W_{1}g_{2}+b_{1},$ (5)

$g_{2}=ELU\left(W_{2}z+b_{2}\right),$ (6)

in which the GRN’s input is given by vector $z$, $g_{1}$ and $g_{2}$ are
intermediate layers, and $ELU$ is the Exponential Linear Unit [4] activation
function, and in which the process of filtering non-linear contributions is
carried out via a Gated Linear Unit (GLU) [5], which provides the Portfolio
Transformer with the ability to scale down the amount of non-linear processing
and default to a simpler model when, for example, the dataset is small or
highly noisy.

###  2.4 Time2Vec Embedding

Since the attention layers of a Transformer do not make use of recurrence,
they cannot inherently capture any information about the relative position of
each element in a sequence. In the original model [19] this information is
injected via positional encoding. The Portfolio Transformer, however,
implements the time encoding proposed by Kazemi et al. [10], that takes the
following form:

$Time2Vec\left(t\right)[i]=\begin{cases}\omega_{i}t+\varphi_{i}&\text{if
$i=0$}\\\ sin\left(\omega_{i}t+\varphi_{i}\right)&\text{if 1 $\leq$ i $\leq$
k.}\end{cases}$ (7)

The temporal signal represented by $t$ in equation (7) is decomposed into a
set of frequencies $\omega$ and phase shifts $\varphi$. This time
decomposition technique is closely related to Fourier transforms, but instead
of using a fixed set of values, all frequencies and phase shifts are learnable
parameters. It should be noted that the use of $sine$ as the activation
function enables the Portfolio Transformer to capture periodic behaviors in
data.

##  3 Methodology

###  3.1 Portfolio Transformer Architecture

The network architecture of the Portfolio Transformer is shown in Fig. 1. It
consists of four main building blocks: input layer (Time2Vec embedding,
section 2.4), encoder and decoder (sections 2.2 and 2.3), and output layer.
Each of these blocks will now be discussed in turn.

Figure 1: The Portfolio Transformer: model architecture.

####  3.1.1 Input Layer.

Each sequence position, denoted by vector $x_{t}$ in Fig. 1, contains
concatenated returns of all $N$ assets in a given portfolio on day $t$. There
are $\tau$ such vectors (per encoder and decoder block) stacked together to
form the input matrix $X$ of dimension ($\tau$, $N$). Time2Vec embedding is
then used to extract time-encoded cross-sectional features.

####  3.1.2 Encoder.

The Portfolio Transformer uses a stack of four identical encoder layers.
Inside each of these layers, the time-encoded input is first processed by a
multi-headed attention mechanism, where the number of attention heads $h$ is
selected during hyperparameter optimization. A series of gating mechanisms are
then applied, using a GRN module that determines the ideal amount of non-
linearity. A residual connection [6], followed by layer normalization [1], are
additionally applied to these two sub-layers, as shown in Fig. 1.

####  3.1.3 Decoder.

The decoder block is also composed of four identical layers, each of which
contains two multi-headed attention modules. The first one uses masking to
ensure predictions made by the Portfolio Transformer depend only on data from
preceding time steps. The second one allows the decoder to attend to the
output of the encoder stack, which provides a much more nuanced representation
of the data and the ability to learn longer-term dependencies among portfolio
assets. As in the case of the encoder, a GRN is used to remove any unnecessary
complexity and there is a residual connection around each sub-layer followed
by layer normalization.

####  3.1.4 Output Layer.

The Portfolio Transformer allows for short-selling via a specialized output
layer that implements a compound function proposed by Zhang et al. [23].
First, the output of the final decoder layer is processed by a fully-connected
layer. The resulting vector ($s_{i,t}$) is then used to compute final
portfolio weights, given by

$w_{i,t}=sign\left(s_{i,t}\right)\times softmax\left(s_{i,t}\right)\triangleq
sign\left(s_{i,t}\right)\times\frac{e^{s_{i,t}}}{\sum_{j=1}^{N}e^{s_{j,t}}}.$
(8)

The use of the softmax operation in equation (8) ensures that, while portfolio
weights can be either positive or negative (the latter allowing for short-
selling), the sum of their absolute values always remains equal to one.

###  3.2 Loss Function

The objective of the PT model as currently implemented (other objective
functions being possible) is to learn the asset distribution that maximizes
risk-adjusted returns as measured by the Sharpe ratio, which is defined below
as expected portfolio return divided by its volatility:

$SR=\frac{E(R_{P})}{\sqrt{E(R_{P}^{2})-(E(R_{P}))^{2}}}.$ (9)

Since transaction costs can significantly diminish the performance of
allocation strategies with high turnover, the Portfolio Transformer uses cost-
adjusted portfolio returns

$R_{P,t}=\sum_{i}^{N}w_{i,t-1}\times
r_{i,t}-C\times\sum_{i}^{N}\left|w_{i,t-1}-w_{i,t-2}\right|$ (10)

in order to find solutions that account for trading costs, where $C$ is a
constant cost rate, set to a realistic value of two basis points (2 bps),
$w_{i,t-1}$ represents the weight of asset $i$ on day $t-1$, and $r_{i,t}$
denotes the realized arithmetic return of asset $i$ from day $t-1$ to day $t$,
computed using asset prices $P_{i,t}$ and $P_{i,t-1}$ as follows:

$r_{i,t}=\frac{P_{i,t}}{P_{i,t-1}}-1.$ (11)

The expected portfolio return, denoted by $E\left(R_{P}\right)$ in equation
(9), is obtained by taking an average of all portfolio returns over a trading
period of length $\tau$:

$E\left(R_{P}\right)=\frac{1}{\tau}\sum_{t=1}^{\tau}R_{P,t}.$ (12)

Finally, since the PT allows for short-selling but no leverage, the portfolio
positions are constrained by $w_{i,t}\in[-1,1]$ and
$\sum_{i}^{N}\left|w_{i,t}\right|=1$, which is achieved through the use of the
compound function in equation (8).

###  3.3 Training & Model Calibration

Data is split into training and test segments using an expanding window
approach, where initially all data points before the end of 2015 are used for
training and the out-of-sample test is carried out on observations recorded in
2016. The training window is then extended to include the year 2016 and the
model is tested on the subsequent year (2017), and so on. This way, the model
is retrained every year, with all available historical data being used to
update the network parameters. Portfolio positions are adjusted on a daily
basis, and a transaction cost rate of 2 bps is used during performance
evaluation.

The PT network is trained via mini-batch stochastic gradient descent (with
batch size being among the hyperparameters) using the Adam optimizer [11]. For
model calibration purposes, and to control for overfitting, 10% of any
training segment is set aside as a separate validation set. Hyperparameter
optimization, conducted using 100 iterations of random grid search, is
performed only on the validation set, ensuring that the model has access to
test data only during the performance evaluation stage. To further improve the
model’s ability to generalize, early stopping is implemented. The Portfolio
Transformer was developed using the TensorFlow framework and all experiments
were conducted on NVIDIA’s Tesla P100 16GB GPU with 55GB of RAM memory.

###  3.4 Datasets Used

The efficacy of the Portfolio Transformer is demonstrated on three datasets
containing daily price observations. The first of these datasets starts in
2006 and is composed of the same four Exchange Traded Funds (ETFs) used in the
LSTM-based experiments carried out by Zhang et al. [24]: AGG (aggregate bond
index), DBC (commodity index), VIX (volatility index), and VTI (US stocks
index). The second dataset, which starts in 2002, is composed of 24 continuous
commodity futures contracts, including metals, agricultural products, and
energy commodities such as oil and natural gas. Finally, the PT model is
tested on daily observations of 500 stocks based in the US and listed on
NASDAQ. This last dataset starts in 1996 and aims to demonstrate the model’s
performance on a large portfolio of hundreds of instruments.

###  3.5 Benchmark Models

Four algorithms are implemented as benchmarks: (1) mean-variance optimization
(MV) [14], a classical two-step portfolio selection procedure with a moving
window of 50 days used to estimate expected asset returns and covariances; (2)
XGBoost [2], included because gradient-boosted decision trees perform very
well in applied machine learning competitions; (3) a multilayer perceptron
(MLP), as a universal function approximator that can capture highly non-linear
dependencies; and (4) the LSTM architecture of [24], a high-performance
recurrent architecture that represents the previous state of the art.

##  4 Results

###  4.1 Performance Comparison: Full Investment Horizon

The performance of the Portfolio Transformer is compared to that of the
benchmark algorithms of section 3.5 using a number of metrics that aim to
capture portfolio risk level through annualized volatility (Vol.) and maximum
drawdown (MDD), profitability via annualized returns (Returns) and percentage
of positive returns (% of + Ret), and risk-adjusted performance using
annualized Sharpe, Sortino and Calmar ratios.

Table 1: Experimental results for different algorithms and datasets. | Returns | Vol. | Sharpe | Sortino | MDD | Calmar |  % of +Ret
---|---|---|---|---|---|---|---
Panel A: ETFs |  |  |  |  |  |  |
MV | 0.004 | 0.122 | 0.012 | 0.270 | 0.120 | 0.836 | 0.497
XGBoost | 0.100 | 0.140 | 0.657 | 1.240 | 0.116 | 2.591 | 0.496
MLP | 0.135 | 0.128 | 0.923 | 1.789 | 0.087 | 2.225 | 0.498
LSTM | 0.215 | 0.133 | 1.539 | 2.830 | 0.096 | 3.621 | 0.535
Portfolio Transformer | 0.138 | 0.067 | 2.252 | 4.093 | 0.036 | 4.773 | 0.548
Panel B: Commodities |  |  |  |  |  |  |
MV | 0.008 | 0.053 | 0.174 | 0.342 | 0.059 | 0.682 | 0.519
XGBoost | 0.015 | 0.059 | 0.277 | 0.420 | 0.063 | 0.267 | 0.505
MLP | 0.026 | 0.056 | 0.479 | 0.727 | 0.057 | 0.797 | 0.515
LSTM | 0.038 | 0.031 | 1.182 | 1.852 | 0.023 | 2.108 | 0.528
Portfolio Transformer | 0.174 | 0.108 | 1.506 | 2.304 | 0.077 | 2.272 | 0.543
Panel C: Stocks |  |  |  |  |  |  |
MV | 0.079 | 0.126 | 0.694 | 1.126 | 0.106 | 1.386 | 0.523
XGBoost | 0.101 | 0.118 | 0.923 | 1.352 | 0.080 | 1.491 | 0.533
MLP | 0.089 | 0.121 | 0.767 | 1.102 | 0.087 | 1.236 | 0.534
LSTM | 0.111 | 0.077 | 1.456 | 2.155 | 0.056 | 2.561 | 0.565
Portfolio Transformer | 0.334 | 0.147 | 2.001 | 3.440 | 0.091 | 4.824 | 0.566
Figure 2: Comparison of cumulative returns.

Panel A in Table 1 shows the results for the portfolio of ETFs, where the
Portfolio Transformer outperforms the benchmarks on all but one metric,
annualized returns. These are higher for the LSTM model, but this comes at a
cost of increased volatility; when this is accounted for, by using the Sharpe
ratio, the PT offers much higher _risk-adjusted_ returns. Results presented in
Panels B (commodities) and C (stocks) show that in both cases the PT achieves
the best risk-adjusted performance, while also delivering the highest
annualized returns. The PT is beaten by the LSTM in these cases on only two
out of seven metrics, volatility and maximum drawdown. In the former case,
this is compensated by the PT’s higher Sharpe ratio and in the latter by the
higher Calmar ratio, relative to the LSTM; the Sharpe ratio is a measure of
portfolio return adjusted by volatility, and Calmar ratio of portfolio return
adjusted by maximum drawdown, arguably of more relevance than volatility and
maximum drawdown per se.

The cumulative return plots of Fig. 2 demonstrate the superior performance of
the Portfolio Transformer over the whole investment horizon. The PT generates
the highest cumulative returns and offers a reasonable risk profile for all
three datasets. The second-best performing model is the LSTM, suggesting that
time dependencies learned through recurrence, in case of the LSTM model, or
through an attention mechanism, in case of the PT, are very useful in a
portfolio optimization setting. However, it can be seen that the LSTM model
struggles during the COVID-19 crisis (first quarter of 2020), while the
attention-driven Portfolio Transformer shows a much quicker response to this
sudden market regime change. The performances of the XGBoost and MLP models
are comparable, but lag considerably behind those of the LSTM and the PT,
while mean-variance optimization (MV) is by far the worst-performing
algorithm, suggesting that highly inaccurate asset weightings are generated by
this classical two-step procedure.

###  4.2 Performance Comparison: COVID-19 Crisis

Figure 3: Rolling Sharpe ratio (12-month) on the ETF dataset.

The above-mentioned difference between the performances of the PT and LSTM
models during this period of extreme market volatility is further illustrated
in Fig. 3, which shows the 12-month rolling Sharpe ratio of the two models on
the ETF dataset. The LSTM suffers a large fall in its risk-adjusted returns
during this period and there are times when it even drops down below zero. The
Portfolio Transformer, on the other hand, shows a more stable behavior, with
its rolling Sharpe ratio staying mostly well above one and delivering
outstanding risk-adjusted performance during the Bull market that followed the
crisis.

##  5 Conclusions

This work has introduced the Portfolio Transformer (PT), which directly
optimizes risk-adjusted returns using a novel end-to-end attention-based
architecture with specialized time-encoding layers and gating mechanisms. By
incorporating transaction costs directly into its loss function the PT model
can account for trading cost constraints faced by investors. The results
demonstrate that the PT model delivers exceptional risk-adjusted performance,
in this respect outperforming all benchmark algorithms on three different
datasets with varying portfolio sizes. Due to its full encoder-decoder
configuration and its attention mechanism the Portfolio Transformer is able to
learn long-term dependencies and as a result can react more quickly to
changing market regimes, as demonstrated by its response to the COVID-19
crisis. Turning to future work, one extension of the current model could study
the Portfolio Transformer’s performance under an attention mechanism different
from the scaled dot-product attention currently used. In addition, while it is
popular with industry practitioners, the Sharpe ratio is only one of many
possible objective functions that could be optimized, and subsequent work will
consider alternative metrics and their effect on the overall portfolio
performance.

## References

  * [1] Ba, J.L., Kiros, J.R., Hinton, G.E.: Layer normalization. arXiv preprint arXiv:1607.06450 (2016)
  * [2] Chen, T., Guestrin, C.: XGBoost: A scalable tree boosting system. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. pp. 785–794 (2016)
  * [3] Chen, W., Zhang, H., Mehlawat, M.K., Jia, L.: Mean-variance portfolio optimization using machine learning-based stock price prediction. Applied Soft Computing 100, 106943 (2021)
  * [4] Clevert, D.A., Unterthiner, T., Hochreiter, S.: Fast and accurate deep network learning by exponential linear units (elus). arXiv preprint arXiv:1511.07289 (2015)
  * [5] Dauphin, Y.N., Fan, A., Auli, M., Grangier, D.: Language modeling with gated convolutional networks. In: International Conference on Machine Learning. pp. 933–941 (2017)
  * [6] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 770–778 (2016)
  * [7] Heaton, J.B., Polson, N.G., Witte, J.H.: Deep learning for finance: deep portfolios. Applied Stochastic Models in Business and Industry 33(1), 3–12 (2017)
  * [8] Hochreiter, S., Bengio, Y., Frasconi, P., Schmidhuber, J.: Gradient flow in recurrent nets: the difficulty of learning long-term dependencies. In: A field guide to dynamical recurrent neural networks, chap. 14, pp. 237–374. IEEE Press (2001)
  * [9] Hochreiter, S., Schmidhuber, J.: Long short-term memory. Neural Computation 9(8), 1735–1780 (1997)
  * [10] Kazemi, S.M., Goel, R., Eghbali, S., Ramanan, J., Sahota, J., Thakur, S.: Time2Vec: Learning a vector representation of time. arXiv preprint arXiv:1907.05321 (2019)
  * [11] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014)
  * [12] Kisiel, D., Gorse, D.: A meta-method for portfolio management using machine learning for adaptive strategy selection. In: International Conference on Computational Intelligence and Intelligent Systems 2021. pp. 67–71. ACM (2021)
  * [13] Lim, B., Arık, S.Ö., Loeff, N., Pfister, T.: Temporal fusion transformers for interpretable multi-horizon time series forecasting. International Journal of Forecasting 37(4), 1748–1764 (2021)
  * [14] Markowitz, H.: Portfolio selection. The Journal of Finance 7(1), 77–91 (1952)
  * [15] Michaud, R.O., Michaud, R.O.: Efficient asset management: a practical guide to stock portfolio optimization and asset allocation. Oxford University Press (2008)
  * [16] Moody, J., Wu, L., Liao, Y., Saffell, M.: Performance functions and reinforcement learning for trading systems and portfolios. Journal of Forecasting 17(5-6), 441–470 (1998)
  * [17] Shen, G., Tan, Q., Zhang, H., Zeng, P., Xu, J.: Deep learning with gated recurrent unit networks for financial sequence predictions. Procedia Computer Science 131, 895–903 (2018)
  * [18] Tino, P., Schittenkopf, C., Dorffner, G.: Financial volatility trading using recurrent neural networks. IEEE Transactions on Neural Networks 12(4), 865–874 (2001)
  * [19] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N.: Attention is all you need. Advances in Neural Information Processing Systems 30 (2017)
  * [20] Wolf, T., Debut, L., Sanh, V., Chaumond, J., Delangue, C., Moi, A.: Huggingface’s transformers: State-of-the-art natural language processing. arXiv preprint arXiv:1910.03771 (2019)
  * [21] Wood, K., Giegerich, S., Roberts, S., Zohren, S.: Trading with the momentum transformer: An intelligent and interpretable architecture. arXiv preprint arXiv:2112.08534 (2021)
  * [22] Xu, K., Zhang, Y., Ye, D., Zhao, P., Tan, M.: Relation-aware transformer for portfolio policy learning. In: Proceedings of the Twenty-Ninth International Conference on International Joint Conferences on Artificial Intelligence. pp. 4647–4653 (2021)
  * [23] Zhang, C., Zhang, Z., Cucuringu, M., Zohren, S.: A universal end-to-end approach to portfolio optimization via deep learning. arXiv preprint arXiv:2111.09170 (2021)
  * [24] Zhang, Z., Zohren, S., Roberts, S.: Deep learning for portfolio optimization. The Journal of Financial Data Science 2(4), 8–20 (2020)


# Stochastic hydrodynamics and hydro-kinetics: Similarities and differences

Aritra De  School of Physics & Astronomy, University of Minnesota,
Minneapolis, MN 55455, USA  Chun Shen  Department of Physics and Astronomy,
Wayne State University, Detroit, MI 48201, USA  RIKEN BNL Research Center,
Brookhaven National Laboratory, Upton, NY 11973, USA  Joseph I. Kapusta
School of Physics & Astronomy, University of Minnesota, Minneapolis, MN 55455,
USA

###### Abstract

The hydro-kinetic formalism has been used as a complementary approach to
solving the Stochastic Differential Equations (SDE) corresponding to noisy
hydrodynamics. The hydro-kinetic formalism consists of a deterministic set of
relaxation type equations that tracks the evolution of 2-point correlation
functions of stochastic hydrodynamic quantities. Hence they are comparatively
easier to solve than the SDEs, which are computationally intensive and need to
deal with arbitrarily large gradients. This work compares the two approaches
for the propagation and diffusion of conserved charge fluctuations in the
Bjorken hydrodynamic model. For white noise, the two approaches agree. For
colored Catteneo noise, which is causal, the two approaches diverge. This is
because white noise only induces two-point correlations, while Catteneo noise
also induces higher-order correlations. This difference is quantified from the
effects of causal evolution and influence from higher-order correlations
induced by the Catteneo noise.

##  I Introduction

Experimental facilities such as the Relativistic Heavy-Ion Collider (RHIC) and
the Large Hadron Collider (LHC) collide nuclei at high energies to study and
quantify properties of the quark-gluon plasma. Relativistic hydrodynamics
provides a macroscopic and quantitative description of the fireball evolution
in heavy-ion collisions Heinz and Snellings (2013); Gale  _et al._ (2013);
Shen and Yan (2020). The applicability of hydrodynamics is justified if the
mean free paths of the particles are small compared to the distances over
which thermodynamic quantities vary. It has been successfully used to
constrain various transport properties of the quark-gluon plasma such as
specific shear and bulk viscosity and thermal conductivity Pratt  _et al._
(2015); Bernhard  _et al._ (2016, 2019); Everett  _et al._ (2021a, b); Nijs
_et al._ (2021a, b). Hybrid frameworks that combine hydrodynamics with a
hadronic transport model can reproduce various flow measurements with
successful predictive power Song  _et al._ (2011); Shen  _et al._ (2016);
Putschke  _et al._ (2019); Schenke  _et al._ (2020).

Relativistic hydrodynamics involves differential equations for the energy-
momentum tensor and the charge currents. Ripples caused by stochastic thermal
fluctuations will give corrections to macroscopic variables of hydrodynamics
Young  _et al._ (2015). The fluctuation-dissipation theorem Kapusta  _et al._
(2012) relates dissipative properties to the amplitude of hydrodynamic
fluctuations. Hence, study of fluctuations can be leveraged to study quark-
gluon plasma transport properties Pratt and Young (2017); Pratt and Plumberg
(2019). Specifically, fluctuations are known to become important near critical
points Plumberg and Kapusta (2017); Bzdak  _et al._ (2020); Bluhm  _et al._
(2020); An  _et al._ (2022). This has led to studies that investigate
fluctuations in conserved quantities, such as electric charge, baryon number,
and strangeness on an event-by-event basis Kapusta and Plumberg (2018); Pratt
_et al._ (2018); Pratt and Plumberg (2020); De  _et al._ (2020). Recently,
there have also been studies of dynamic critical phenomena near a QCD critical
point Nahrgang and Bluhm (2020); Rajagopal  _et al._ (2020); Du  _et al._
(2020); An  _et al._ (2022); Du  _et al._ (2021). Non-Gaussianities (higher-
order cumulants) also become important near critical points as they are more
sensitive to large correlation lengths Bzdak  _et al._ (2020).

It is important to systematically study fluctuations in heavy-ion collision
systems. Adding noise to hydrodynamic equations transforms the hydrodynamical
equations from partial differential equations to stochastic differential
equations. Historically, stochastic differential equations have been notorious
to deal with Young  _et al._ (2015); Murase (2015); Singh  _et al._ (2019);
Sakai  _et al._ (2020, 2021). Broadly, there have been two methods to study
the effect of fluctuations. The first is the direct method of studying
stochastic differential equations through analytical and numerical techniques
Kapusta and Plumberg (2018); Young  _et al._ (2015); Singh  _et al._ (2019);
Sakai  _et al._ (2020); Pratt  _et al._ (2018); Pratt and Plumberg (2020); De
_et al._ (2020). One can simulate stochastic differential equations like usual
partial differential equations over millions of different realizations of the
noise. Statistical properties of those solutions are studied to understand the
underlying physical process. Alternatively, some authors have suggested the
use of a set of deterministic kinetic equations to capture the behavior of
time evolution of $n$-point correlations of hydrodynamic variables Akamatsu
_et al._ (2017); Stephanov and Yin (2018); Akamatsu  _et al._ (2019); An  _et
al._ (2019); Martinez  _et al._ (2019); An  _et al._ (2020); Martinez and
Schäfer (2017). This approach has come to be known as the hydro-kinetic
approach. The method relies on the truncation of higher moments from noisy
hydrodynamic calculations and holds the advantage over direct numerical
studies of SDEs by being computationally less intensive.

It seems imperative that a direct comparison between the hydro-kinetic
approach and the stochastic hydro formalism needs to be made which is lacking
in the current literature. In this paper, we compare the two approaches for a
system with just one conserved charge in an expanding fluid background, namely
the Bjorken system Bjorken (1983). At the top RHIC and LHC energies, the
charge current decouples from the energy-momentum conservation because the net
charge is small and therefore can be studied separately. The system with a
single conserved charge is simple enough to allow for analytical calculations.
Additionally, the conserved baryon density is believed to be an order
parameter for a possible critical point in the QCD phase diagram Bzdak  _et
al._ (2020). The usual diffusion equation has an infinite signal propagation
speed. To make it finite, we also study the Catteneo equation and which
introduces so-called colored noise in the hydrodynamical equations Kapusta and
Young (2014a, b). In the presence of colored noise, the hydrodynamics become
causal. We show the comparison of hydro-kinetics to causal hydrodynamics as
well.

The paper is organized as follows. We give a concise description of the
dynamical evolution of one conserved charge in a 1+1D Bjorken system in Sec.
II. In Sec. III, we derive the hydro-kinetic equations and compare them with
the stochastic calculation in the white noise limit. We extend the comparison
for the colored noise case to hydro-kinetics in Sec. IV. Sec. V summarizes our
concluding remarks. Complementary calculations without transverse dynamics are
provided in the Appendix A.

##  II Stochastic diffusion in boost-invariant hydrodynamics

The assumption of longitudinal boost invariance implies that the initial
conditions for local variables are functions of the proper time
$\tau\equiv\sqrt{t^{2}-z^{2}}$ and the transverse coordinates $x$ and $y$.
They do not depend on spacetime rapidity
$\xi\equiv\tanh^{-1}\left(z/t\right)$. We define fluid velocity $u^{\mu}$ to
be the velocity of the energy-momentum current, known as the Landau-Lifshitz
frame,

$T^{\mu\nu}u_{\nu}=eu^{\mu}.$ (1)

We consider a single species of conserved charge in the fluid. This charge
current takes the form

$J^{\mu}=n_{Q}u^{\mu}+\Delta J^{\mu}+I^{\mu},$ (2)

where $n_{Q}$ is the charge density in the local rest frame of the fluid cell,
$\Delta J^{\mu}$ is the diffusion current, and $I^{\mu}$ is the noise current.
The fluctuation-dissipation theorem dictates that dissipation is always
accompanied by stochastic noises for any equilibrium system. In the
conventional, lowest order in derivatives diffusion equation, $\Delta J^{\mu}$
takes the form

$\Delta
J^{\mu}=D_{Q}\nabla^{\mu}n_{Q}=\sigma_{Q}\,T\,\nabla^{\mu}\left(\frac{\mu_{Q}}{T}\right),$
(3)

where $\mu_{Q}$ is the charge chemical potential, $\sigma_{Q}$ is the charge
conductivity, and $\nabla^{\mu}$ is the transverse derivative

$\nabla^{\mu}\equiv\Delta^{\mu\nu}\partial_{\nu}=\partial^{\mu}-u^{\mu}(u\cdot\partial)$
(4)

with the projection operator $\Delta^{\mu\nu}=g^{\mu\nu}-u^{\mu}u^{\nu}$. One
can verify that the charge current conservation equation
$\partial_{\mu}J^{\mu}=0$ in the local rest frame where
$u^{\mu}=(1,\mbox{\boldmath$0$})$ follows the usual diffusion equation

$\left(\frac{\partial}{\partial t}-D_{Q}\nabla^{2}\right)n_{Q}=0\,.$ (5)

The diffusion constant $D_{Q}$ and charge conductivity are related by the
Einstein relation $D_{Q}=\sigma_{Q}/\chi_{Q}$, where $\chi_{Q}$ is the
electric charge susceptibility defined by

$\chi_{Q}=\frac{\partial n_{Q}(T,\mu_{Q})}{\partial\mu_{Q}}\,.$ (6)

Following Refs. Kapusta and Plumberg (2018); De  _et al._ (2020), we assume
$\chi_{Q}\propto T^{2}$. The amplitudes for dissipation and stochastic noise
are related by the fluctuation-dissipation theorem which is justified because
the fluid description assumes the system is close to local equilibrium. In
other words, modes disturbed by $I^{\mu}$ are equilibrated by the viscous term
$\Delta J^{\mu}$. The 1-point function of the noise $\langle I^{\mu}\rangle$
vanishes while the 2-point functions are determined by the fluctuation-
dissipation theorem. Assuming white noise

$\langle I^{\mu}(x)\rangle=0,\quad\langle
I^{\mu}(x_{1})I^{\nu}(x_{2})\rangle=2\sigma_{Q}T\,\Delta^{\mu\nu}\delta(x_{1}-x_{2}).$
(7)

Such a noise function whose 2-point correlation is proportional to the Dirac
$\delta$-function in both space and time is referred to as white noise.

We adopt Milne coordinates which are proper time and space-time rapidity,

$\displaystyle\begin{aligned} t&=&\tau\cosh\xi\qquad z&=\tau\sinh\xi\\\
\tau&=&\sqrt{t^{2}-z^{2}}\qquad\xi&=\tanh^{-1}\left(\frac{z}{t}\right).\\\
\end{aligned}$ (8)

The transverse coordinates $x$ and $y$ are the same as in Cartesian
coordinates. The non-zero components of the background flow velocity are

$u^{0}=\cosh\xi\qquad u^{z}=\sinh\xi\,.$ (9)

The transverse derivatives in the $t$ and $z$ directions are

$\displaystyle\nabla^{0}=-\frac{\sinh\xi}{\tau}\frac{\partial}{\partial\xi}\qquad\nabla^{3}=-\frac{\cosh\xi}{\tau}\frac{\partial}{\partial\xi}$
(10)

with $\,u\cdot\partial=\partial/\partial\tau$. The entropy density $s$ can be
factored out from the noise current $I^{\mu}$ so as to capture the fluctuating
part in the random dimensionless scalars $f$, $g_{1}$ and $g_{2}$.

$\displaystyle I^{t}$ $\displaystyle=$ $\displaystyle
s(\tau)\;f(\xi,\tau,x,y)\sinh\xi$ $\displaystyle I^{z}$ $\displaystyle=$
$\displaystyle s(\tau)\;f(\xi,\tau,x,y)\cosh\xi$ $\displaystyle I^{x}$
$\displaystyle=$ $\displaystyle s(\tau)\;g_{1}(\xi,\tau,x,y)$ $\displaystyle
I^{y}$ $\displaystyle=$ $\displaystyle s(\tau)\;g_{2}(\xi,\tau,x,y)$ (11)

Here $f,g_{1},g_{2}$ are uncorrelated random scalar functions with the
following 2-point correlation in Cartesian coordinates

$\displaystyle\langle f(x_{1})f(x_{2})\rangle$ $\displaystyle=$
$\displaystyle\langle g_{1}(x_{1})g_{1}(x_{2})\rangle=\langle
g_{2}(x_{1})g_{2}(x_{2})\rangle$ (12) $\displaystyle=$
$\displaystyle\frac{2\sigma_{Q}T}{s^{2}}\delta^{4}(x_{1}-x_{2})\;\;\;$

and

$\langle f(x_{1})g_{1}(x_{2})\rangle=\langle
g_{2}(x_{1})g_{1}(x_{2})\rangle=\langle f(x_{1})g_{2}(x_{2})\rangle=0.$ (13)

The last equation denotes that the random functions $f$, $g_{1}$ and $g_{2}$
are uncorrelated with each other.

The background conservation equations for the charge density and entropy
density are

$\displaystyle\frac{ds}{d\tau}+\frac{s}{\tau}=0\;\;$
$\displaystyle\Rightarrow$
$\displaystyle\;\;s(\tau)=\frac{s_{i}\tau_{i}}{\tau}$ (14)
$\displaystyle\frac{dn_{Q}}{d\tau}+\frac{n_{Q}}{\tau}=0\;\;$
$\displaystyle\Rightarrow$
$\displaystyle\;\;n_{Q}(\tau)=\frac{n_{i}\tau_{i}}{\tau}\,.$ (15)

The $s_{i}$ and $n_{i}$ are the entropy and charge densities at some initial
time $\tau_{i}$. We take the initial charge density $n_{i}$ to be zero, hence
the average charge density for subsequent times is zero as well. Please note
that because the stochastic fluctuations are small perturbations compared to
the background quantities, we neglect their feedback corrections to the
background fluid velocity in this study.

We now have all the pieces of the charge current $J^{\mu}$ to write the charge
current conservation equation $\partial_{\mu}J^{\mu}=0$. It is convenient to
define the variable $X=\tau\delta n_{Q}$. The charge conservation equation
then becomes

$\displaystyle\frac{1}{\tau}\frac{\partial
X}{\partial\tau}-\frac{D_{Q}}{\tau}\left[\nabla_{\perp}^{2}+\frac{1}{\tau^{2}}\frac{\partial^{2}}{\partial\xi^{2}}\right]X=-\frac{s}{\tau}\frac{\partial
f}{\partial\xi}-s\frac{\partial g_{1}}{\partial x}-s\frac{\partial
g_{2}}{\partial y}$ (16)

where $\nabla^{2}_{\perp}=\partial^{2}/\partial x^{2}+\partial^{2}/\partial
y^{2}$. We do a Fourier transform

$X(\xi,x,y,\tau)=\int\frac{dk_{\xi}}{2\pi}\frac{d^{2}\mathbf{k_{\perp}}}{2\pi}e^{i(k_{\xi}\xi+\mathbf{k_{\perp}.x_{\perp}})}\;\tilde{X}(k_{\xi},\mathbf{k_{\perp}},\tau).$
(17)

In Fourier space the charge conservation equation becomes

$\frac{\partial\tilde{X}}{\partial\tau}+D_{Q}\mathbf{k}_{\perp}^{2}\tilde{X}+\frac{D_{Q}k_{\xi}^{2}}{\tau^{2}}\tilde{X}=-is\left[k_{\xi}\tilde{f}+\tau
k_{x}\tilde{g_{1}}+\tau k_{y}\tilde{g_{2}}\right].$ (18)

The Green’s function of the homogeneous part of the above equation is

$\tilde{G}(k_{\xi},\mathbf{k_{\perp}},\tau_{f},\tau_{i})=e^{-D_{Q}\mathbf{k}_{\perp}^{2}(\tau_{f}-\tau_{i})}e^{-D_{Q}k_{\xi}^{2}\left(\frac{1}{\tau_{i}}-\frac{1}{\tau_{f}}\right)}.$
(19)

Then

$\displaystyle\tilde{X}(k_{\xi},k_{x},k_{y},\tau_{f})$ $\displaystyle=$
$\displaystyle-i\int^{\tau_{f}}_{\tau_{i}}d\tau s\left[k_{\xi}\tilde{f}+\tau
k_{x}\tilde{g_{1}}+\tau k_{y}\tilde{g_{2}}\right]$ (20) $\displaystyle\times$
$\displaystyle
e^{-D_{Q}\mathbf{k}_{\perp}^{2}(\tau_{f}-\tau_{i})}e^{-D_{Q}k_{\xi}^{2}\left(\frac{1}{\tau_{i}}-\frac{1}{\tau_{f}}\right)}.$

The 2-point correlation function can now be calculated as

$\displaystyle\langle\delta\tilde{n}_{Q}(k_{\xi},\mathbf{k_{\perp}},\tau_{f})\,\delta\tilde{n}_{Q}(-k_{\xi}^{\prime},-\mathbf{k^{\prime}_{\perp}},\tau_{f})\rangle$
(21) $\displaystyle=$ $\displaystyle\frac{4\pi
D_{Q}\chi_{f}T_{f}}{\tau_{f}}\int
d\tau^{\prime}\left[\frac{k_{\xi}^{2}}{\tau^{\prime
2}}+k_{x}^{2}+k_{y}^{2}\right]\hskip 56.9055pt$ $\displaystyle\times$
$\displaystyle
e^{-2D_{Q}\mathbf{k}_{\perp}^{2}(\tau_{f}-\tau^{\prime})}e^{-2D_{Q}k_{\xi}^{2}\left(\frac{1}{\tau^{\prime}}-\frac{1}{\tau_{f}}\right)}$

The above integral can be solved analytically only if $k_{x}=k_{y}=0$, see
Appendix A. We will perform a numerical computation of Eq. (21) in the
following section.

##  III The hydro-kinetic equation for the 2-point correlation function

In this section, we derive and solve the equation of motion for the 2-point
function of stochastic fluctuations. This is referred to as the hydro-kinetic
approach. Let us begin with the derivation of the hydro-kinetic equation by
ignoring the transverse directions $x$ and $y$. The noisy charge conservation
equation in the absence of transverse noise in $k$-space is

$\displaystyle\frac{\partial\delta\tilde{n}_{Q}}{\partial\tau}=-\frac{\delta\tilde{n}_{Q}}{\tau}-\frac{D_{Q}k^{2}}{\tau^{2}}\delta\tilde{n}_{Q}-\frac{iks}{\tau}\tilde{f}.$
(22)

The noise correlation can be calculated from Eq. (12) as

$\langle\tilde{f}(k_{1},\tau_{1})\tilde{f}(k_{2},\tau_{2})\rangle=\frac{4\pi\sigma_{Q}T}{A\tau_{1}s^{2}}\delta(\tau_{1}-\tau_{2})\delta(k_{1}+k_{2}).$
(23)

where $A$ is the transverse area. The evolution of the 2-point function
satisfies

$\displaystyle\frac{d}{d\tau}\langle\delta\tilde{n}_{Q}$ $\displaystyle(k)$
$\displaystyle\delta\tilde{n}_{Q}(-k^{\prime})\rangle$ $\displaystyle=$
$\displaystyle-$
$\displaystyle\frac{2}{\tau}\langle\delta\tilde{n}_{Q}(k)\delta\tilde{n}_{Q}(-k^{\prime})\rangle$
(24) $\displaystyle-$ $\displaystyle\frac{D_{Q}}{\tau^{2}}(k^{2}+k^{\prime
2})\langle\delta\tilde{n}_{Q}(k)\delta\tilde{n}_{Q}(-k^{\prime})\rangle$
$\displaystyle-$
$\displaystyle\frac{iks}{\tau}\langle\tilde{f}(k)\delta\tilde{n}_{Q}(-k^{\prime})\rangle$
$\displaystyle+$
$\displaystyle\frac{ik^{\prime}s}{\tau}\langle\delta\tilde{n}_{Q}(k)\tilde{f}(-k^{\prime})\rangle.$

We can compute the $\langle\tilde{f}\tilde{n}_{Q}\rangle$ up to
$\mathcal{O}(\delta\tilde{n}_{Q}^{2})$ accuracy using the ideal part of the
equation of motion and ignoring the dissipative part.

$\displaystyle\delta\tilde{n}_{Q}=-\int\frac{d\tau}{\tau}(iks)\tilde{f}\quad\Rightarrow\quad\langle\delta\tilde{n}_{Q}\tilde{f}\rangle=-\int\frac{d\tau}{\tau}(iks)\langle\tilde{f}\tilde{f}\rangle$

Then

$\displaystyle\frac{ik^{\prime}s}{\tau}\langle\delta\tilde{n}_{Q}(k)\tilde{f}(-k^{\prime})\rangle$
$\displaystyle=$
$\displaystyle\frac{ik^{\prime}s}{\tau}\int\frac{d\tau}{\tau}(-iks)\langle\tilde{f}(k)\tilde{f}(-k^{\prime})\rangle$
(25) $\displaystyle=$
$\displaystyle\frac{kk^{\prime}4\pi\sigma_{Q}T}{A\tau^{3}}\delta(k-k^{\prime})$

We now define the charge correlation $N(\tau,k)$ as

$N(\tau_{1},k)\delta(k-k^{\prime})\equiv\frac{\langle\delta\tilde{n}_{Q}(k,\tau_{1})\delta\tilde{n}_{Q}(-k^{\prime},\tau_{1})\rangle}{2\pi}.$
(26)

Using this result in Eq. (24), we get the following evolution equation for the
2-point correlation function

$\displaystyle\frac{dN(\tau,k)}{d\tau}=-\frac{2D_{Q}k^{2}}{\tau^{2}}\left(N(\tau,k)-\frac{\chi_{Q}T}{A\tau}\right)-\frac{2}{\tau}N(\tau,k)\qquad$
(27)

Equation (27) is the hydro-kinetic equation for this situation. The term
$\chi_{Q}T/A\tau$ is the asymptotic high-$k$ limit. In coordinate space, this
term represents self-correlation.

We now generalize the hydro-kinetic equation to 3+1 D. The way we do it is by
recognizing that when we include the transverse directions they have
dimensions whereas the spacetime rapidity is dimensionless. Therefore the $k$
will be replaced by $\mathbf{K}=(\frac{k_{\xi}}{\tau},k_{x},k_{y})$. We obtain
the evolution equation for the 2-point correlation function in 3D as

$\displaystyle\frac{dN}{d\tau}=-2D_{Q}\mathbf{K}^{2}\left(N-\frac{\chi_{Q}T}{A\tau}\right)-\frac{2N}{\tau}.$
(28)

The form of this evolution equation matches that of Ref. Akamatsu  _et al._
(2017); Martinez and Schäfer (2017).

We solve the hydro-kinetic equation Eq. (28) numerically and compare it with
the integral Eq. (21). Because Eq. (28) shows the two-point correlation
function $N$ only depends on the magnitude of $K$, it is sufficient for us to
show the comparisons at two values of $k_{x}$ and $k_{y}$ in Fig. 1.

Figure 1: (Color online) Comparison of momentum-space 2-point correlation
functions for charge density fluctuations from the hydro-kinetics approach and
stochastic hydrodynamics. (a) is for $k_{x}=k_{y}=0.1$ fm-1 and (b) is for
$k_{x}=k_{y}=0.5$ fm-1. Both plot are at final time $\tau_{f}=6.3$ fm/$c$ with
$A=1$ fm2.

The hydro-kinetic result matches very well with the stochastic solutions. This
agreement between the two formalism is expected because the noise functions in
Eqs. (12) and (13) only introduce 2-point correlations in the white noise
limit. There is an emergent length scale at $k_{\xi}=k_{*}$ from the
competition between longitudinal expansion and local equilibration Akamatsu
_et al._ (2017). We can equate the viscous dissipation rate to the expansion
rate to get $k_{*}$

$\displaystyle
2D_{Q}\left(k_{x}^{2}+k_{y}^{2}+\frac{k_{*}^{2}}{\tau_{f}^{2}}\right)=\frac{2}{\tau_{f}}.$
(29)

We set the diffusion constant $D_{Q}$ to be $0.16$ fm. For $k_{x}=k_{y}=0.1$
fm-1 and for the parameters we chose, $k_{*}$ turns out to be 6.21. It is
dimensionless because space-time rapidity is dimensionless. For
$k_{x}=k_{y}=0.5$ fm-1, $k_{*}$ is 4.42. One can see in Fig. 1 that the
2-point correlator has attained the equilibrium value of $\chi_{Q}T/\tau_{f}$
for momentum scales $k_{\xi}\gtrsim 0.5k_{*}$. In the following section, we
will extend our analysis to causal forms of noise, which are also called
colored noise.

##  IV Stochastic noise in causal 3+1D relativistic hydrodynamics

The signal propagation speed in the conventional charge diffusion equation is
infinite, which violates causality in the relativistic regime. In order to
limit signal propagation speeds to less than the speed of light, we introduce
a relaxation time $\tau_{Q}$ in the diffusion equation,

$\left(\frac{\partial}{\partial
t}-D_{Q}\nabla^{2}+\tau_{Q}\frac{\partial^{2}}{\partial t^{2}}\right)n_{Q}=0$
(30)

This equation is called the Catteneo equation Cattaneo  _et al._ (1958) and
the noise associated with it is called the Catteneo noise. One can show that
high frequency waves travel at a speed of $v_{Q}=\sqrt{D_{Q}/\tau_{Q}}$
Kapusta and Young (2014b). The dissipative current is modified to

$\Delta
J^{\mu}=D_{Q}\Delta^{\mu}\left[\frac{1}{1+\tau_{Q}(u\cdot\partial)}\right]n_{Q}.$
(31)

The fluctuation-dissipation theorem dictates the amplitude of the stochastic
noise part of the charge current,

$\langle
I^{i}(x_{1})I^{j}(x_{2})\rangle=\frac{\sigma_{Q}T}{\tau_{Q}}\delta(\mbox{\boldmath$x$}_{1}-\mbox{\boldmath$x$}_{2})\,e^{-|t_{1}-t_{2}|/\tau_{Q}}\,\delta_{ij}.$
(32)

The Dirac $\delta$-function in time is replaced by an exponential decay
function. In the limit $\tau_{Q}\rightarrow 0$ this 2-point function becomes
the Dirac $\delta$-function for the white noise limit in Eq. (12). Such a
noise function with an exponential function replacing the Dirac
$\delta$-function in the 2-point correlation is an example of colored noise.
If we use the previous parametrization of $I^{\mu}$, the auxiliary function
$f$ will have the property

$\displaystyle\langle(\tau_{Q}\,\partial_{\tau_{1}}+1)\tilde{f}(k_{1},\tau_{1})(\tau_{Q}\,\partial_{\tau_{2}}+1)\tilde{f}(k_{2},\tau_{2})\rangle$
$\displaystyle=\frac{4\pi\sigma_{Q}(\tau_{1})T(\tau_{1})}{A\tau_{1}s^{2}(\tau_{1})}\,\delta(\tau_{1}-\tau_{2})\delta(k_{1}+k_{2}).$
(33)

The same is the case for $g_{1}$ and $g_{2}$. The source functions $f,g_{1},$
and $g_{2}$ are still uncorrelated with each other as Eq. (13) still holds.
Noise generated at the same rapidity but at different times is correlated.
Figure 2 illustrates what fluctuation propagation looks like in the presence
of the Catteneo noise De  _et al._ (2020).

Figure 2: (Color online) Illustration of colored noise propagation. The stars
denote a noise source and the bursts are the charge fluctuations resulting
from noise. Reproduced from Ref. De  _et al._ (2020).

Figure 2 also illustrates the correlation introduced in the system due to
colored noise. For a white noise profile, the two star noise sources will not
be correlated and hence the charge fluctuation that has traveled will not have
any correlation. But if we switch on colored noise, there will be a
correlation between the two bursts.

In the Catteneo equation of motion for the charge current, the viscous current
is given by Eq. (31). The evolution equation for charge fluctuations
represented by $X=\tau\delta n_{Q}$ in momentum space is

$\displaystyle\tau_{Q}\frac{\partial^{2}\tilde{X}}{\partial\tau^{2}}$
$\displaystyle+$
$\displaystyle\left(1+\frac{2\tau_{Q}}{\tau}\right)\frac{\partial\tilde{X}}{\partial\tau}+\frac{D_{Q}k_{\xi}^{2}}{\tau^{2}}\tilde{X}+(k_{x}^{2}+k_{y}^{2})D_{Q}\tilde{X}$
(34) $\displaystyle+$
$\displaystyle\frac{2D_{Q}\tau_{Q}(k_{x}^{2}+k_{y}^{2})}{\tau}(1+\tau_{Q}\partial_{\tau})^{-1}\tilde{X}$
$\displaystyle=$ $\displaystyle-
ik_{\xi}s\left[\tau_{Q}\frac{\partial\tilde{f}}{\partial\tau}+\left(1+\frac{\tau_{Q}}{\tau}\right)\tilde{f}\right]$
$\displaystyle-
is\bigg{[}(k_{x}\tilde{g}_{1}+k_{y}\tilde{g}_{2})(\tau+2\tau_{Q})$
$\displaystyle+\tau_{Q}\tau\left(k_{x}\frac{\partial\tilde{g}_{1}}{\partial\tau}+k_{y}\frac{\partial\tilde{g}_{2}}{\partial\tau}\right)\bigg{]}.$

To avoid a third-order differential equation we make the following
approximation

$(1+\tau_{Q}\partial_{\tau})^{-1}\simeq 1-\tau_{Q}\partial_{\tau}.$ (35)

This is a good approximation when $\tau_{Q}<\tau$. The approximated charge
conservation equation is

$\displaystyle\tau_{Q}\frac{\partial^{2}\tilde{X}}{\partial\tau^{2}}$
$\displaystyle+$
$\displaystyle\left(1+2\frac{\tau_{Q}}{\tau}\right)\frac{\partial\tilde{X}}{\partial\tau}+\frac{D_{Q}k_{\xi}^{2}}{\tau^{2}}\tilde{X}+\frac{2D_{Q}\tau_{Q}(k_{x}^{2}+k_{y}^{2})}{\tau}\tilde{X}$
(36) $\displaystyle-$
$\displaystyle\frac{2D_{Q}\tau_{Q}^{2}(k_{x}^{2}+k_{y}^{2})}{\tau}\frac{\partial\tilde{X}}{\partial\tau}+D_{Q}(k_{x}^{2}+k_{y}^{2})\tilde{X}\qquad\qquad$
$\displaystyle=$ $\displaystyle-
ik_{\xi}s\left[\tau_{Q}\frac{\partial\tilde{f}}{\partial\tau}+\left(1+\frac{\tau_{Q}}{\tau}\right)\tilde{f}\right]$
$\displaystyle-
is\bigg{[}(k_{x}\tilde{g}_{1}+k_{y}\tilde{g}_{2})(\tau+2\tau_{Q})$
$\displaystyle+\tau_{Q}\tau\left(k_{x}\frac{\partial\tilde{g}_{1}}{\partial\tau}+k_{y}\frac{\partial\tilde{g}_{2}}{\partial\tau}\right)\bigg{]}.$

Eq. (36) is a complex SDE and is challenging to solve analytically. If we
remove the dependence on the transverse coordinates it is possible to solve
the SDE analytically, as was done in Ref. Kapusta and Plumberg (2018). With
the inclusion of transverse coordinates, we will solve the SDE using direct
numerical simulations on a discrete lattice. To ensure numerical convergence,
we use $d\tau=0.003$ fm/$c$. Detailed numerical schemes are explained in
Appendix B.

To generate a 2-point correlation in the form of Eq. (32) for the noise term
$f(\tau,\xi)$ we need to solve an additional Langevin equation

$f+\tau_{Q}\frac{\partial f}{\partial\tau}=\zeta$ (37)

where $\zeta$ is a white noise such that

$\langle\zeta(\tau)\rangle=0,\qquad\langle\zeta(\tau_{1})\zeta(\tau_{2})\rangle=M(\tau_{1})\delta(\tau_{1}-\tau_{2}).$
(38)

Here $M(\tau)$ is the variance of the white noise $\zeta(\tau)$. The
analytical solution of the Langevin equation is given by

$\langle
f(\tau_{1})f(\tau_{2})\rangle=\frac{M}{2\tau_{Q}}\left[e^{|\tau_{1}-\tau_{2}|/\tau_{Q}}-e^{(2\tau_{i}-\tau_{1}-\tau_{2})/\tau_{Q}}\right].$
(39)

For illustration, Fig. 3 uses $M=1$ fm/$c$ and shows that the analytical
solution and the numerical solution match each other. The noise that we need
in this work is given by the 2-point correlation given in Eq. (33).

Figure 3: (Color online) The time evolution for the variance of fluctuations.
Calculations are averaged over one million events with $\tau_{f}=6.3$ fm/$c$,
$\tau_{i}=0.5$ fm/$c$, $\tau_{Q}=0.4$ fm/$c$ and $M=1$ fm/$c$.

Note that this is the same correlation we need for both $g_{1}$ and $g_{2}$.
Once we generate the noise terms we can simulate Eq. (36) for different values
of $\tau_{Q}$ as shown in Fig. 4.

Figure 4: (Color online) The 2-point correlation function for charge
fluctuations with causal evolution. Results with different $\tau_{Q}$ values
are compared with the white noise limit $\tau_{Q}=0$ fm/$c$. The result is for
$10^{5}$ events and is at final time $\tau_{f}=6.3$ fm/$c$ with $A=1$ fm2.
Figure 5: (Color online) The 4-point correlation functions for the Catteneo
noise $f$ and white noise $\zeta$ with $\tau_{Q}=0.4$ fm/$c$.

In the case of the Catteneo noise, the 2-point correlator shows noticeable
differences from the results in the white noise limit. It is because the
Catteneo noise introduces higher-order correlations that are not tracked by
the hydro-kinetic equation, as shown in Fig. 5111The actual values of the
four-point function depends on the lattice spacing used in the simulation.
With $\Delta x=0.013$ fm, we expect $\langle\zeta\zeta\zeta\zeta\rangle\propto
1/(\Delta x)^{2}\sim 6000$.. These higher-order correlations are generated
from the Langevin equation in Eq. (37). The 4-point correlator of white noise
shown in the figure should be strictly zero except when $\tau=\tau_{f}$. The
numerical simulations reproduce that mathematical statement. The same 4-point
correlator for Catteneo noise need not be zero because of higher-order
correlations not present in white noise. This particular correlator is
sufficient to exemplify that fact. Cumulants could also be used to illustrate
that, but calculating them takes about two orders of magnitude more
computational time to obtain the same accuracy. One can see that for Catteneo
noise there is no fixed equilibrium value for any finite value of $k_{\xi}$.
For a lower signal propagation speed, meaning larger $\tau_{Q}$, the magnitude
of the 2-point correlation function $N$ is smaller at larger values of
$k_{\xi}$. That arises from the nature of the self-correlation term as
described in Ref. De  _et al._ (2020).

The nature of the self-correlation is easiest to understand in the 1+1 D case.
The mathematical expression for self-correlation is given by (Refs. Kapusta
and Plumberg (2018); De  _et al._ (2020)).

$\Big{\langle}\delta n_{Q}(\xi_{1},\tau_{f})\delta
n_{Q}(\xi_{2},\tau_{f})\Big{\rangle}_{\text{self}}=\frac{\tau_{f}}{D_{Q}}\int\frac{dk}{2\pi}e^{ik(\xi_{1}-\xi_{2})}\\\
$
$\times\int\frac{d\tau^{\prime\prime}}{\tau^{\prime\prime}}\;\frac{\tilde{G}(-k,\tau_{f},\tau^{\prime\prime})}{ik}\;\langle
f(\tau_{f})f(\tau^{\prime\prime})\rangle$ (40)

The $\tilde{G}$ is a Green function, see Appendix A. Direct numerical
simulations in Ref. De  _et al._ (2020) showed how the self-correlation looks
for various values of $\tau_{Q}$, as depicted in the Fig. 6.

Figure 6: (Color online) Numerical results for self-correlations for colored
noise with different signal propagation speeds.

One can see that as $\tau_{Q}\to 0$, the self-correlation looks more and more
like a Dirac $\delta$\- function. The absence of the Dirac $\delta$-function
in the case of colored noise means that there is no fixed equilibrium value.
The definition of self-correlation is straightforwardly extendable to the 3+1
D.

Self-correlation can be interpreted as the correlation of a charge fluctuation
generated in $\xi_{1}$ at final time $\tau_{f}$ with another charge
fluctuation generated in $\xi_{1}$ at a previous time but which has traveled
to $\xi_{2}$ at $\tau_{f}$. It is nontrivial for colored noise because colored
noise fluctuations in the same $\xi$ are correlated in time. Figure 7 gives an
illustration of self-correlation.

Figure 7: (Color online) Description of self-correlation. The stars denote a
noise source and the bursts are the charge fluctuations resulting from the
noise. Reproduced from Ref. De  _et al._ (2020).

Higher point correlations are also introduced because of causal diffusion.
Specifically, the non-linear treatment of the white noise introduces the
higher-order correlations.

Given the behavior of the 2-point correlator in the presence of colored noise,
with a peak around $k_{*}$ and a decaying tail, it is interesting to
disentangle the effects of causal diffusion and causal noise. Figure 8 shows
the 2-point correlator as a function of $k_{\xi}$ for the following three
cases:

  1. 1.

Causal diffusion and white noise.

  2. 2.

Usual diffusion and causal noise.

  3. 3.

Usual diffusion and white noise

Figure 8: (Color Online) The 2-point correlation functions of charge
fluctuations in momentum space for three combinations of different types of
diffusion equations and noises to disentangle the effects of causal diffusion
and causal noise. The plot is computed at $\tau_{f}=6.3$ fm/$c$.

First, consider case 1. Causal diffusion is an inefficient process compared to
usual diffusion when it comes to distributing the charge buildup at a
particular space point to neighboring points. In addition, white noise
fluctuations at a particular point at a specific time are uncorrelated with
whatever noise fluctuations that have originated at that point in previous
times. This just means that there is continuous pumping of noise fluctuations
at a space point through its temporal journey irrespective of the amount of
charge buildup there. Causal diffusion does a limited redistribution of the
charge buildup. Hence there is a peak in the 2-point correlation around the
length scale $k_{*}$. This is further illustrated in Fig. 9 which is a plot of
the 2-point correlator in coordinate space.

Figure 9: (Color online) The effects of causal diffusion evolution with white
noise on 2-point charge density correlation in coordinate space with lattice
spacing $d\xi=0.1$.

Now consider case 2. When we incorporate causal noise, the noise fluctuations
that occur at a spatial point are correlated to the noise fluctuations that
have already originated there in the previous moments. Hence a spatial point
that has had a lot of charge buildup, because of recent large charge
fluctuations, will not have large fluctuations at subsequent times because of
the exponentially decaying correlation. In this case, the diffusion, either
usual or causal, has enough time to distribute the charge buildup effectively.
This is further illustrated in Fig. 10.

Figure 10: (Color online) The effects of causal diffusion evolution with
causal noise on 2-point charge density correlation in coordinate space with
lattice spacing $d\xi=0.1$.

Both curves in Fig. 10 use the same causal noise and hence the self-
correlation part for both the curves are the same. The presence of causal
diffusion narrows the inverted Gaussian (which signifies the diffusion). Case
3 is the standard white noise result from Fig. 10 of Ref. De  _et al._ (2020),
specifically, the Dirac $\delta$-function minus a Gaussian. Here, the noise
fluctuations pumped in by white noise is effectively dispersed away with
infinite signal propagation speed by the usual diffusion process.

This brings us to our final point. The fluctuation-dissipation theorem tells
us that causal noise and causal diffusion are inextricably present. Hence we
expect both a peak in the two point correlator and a decaying tail at large
$k_{\xi}$. Both may have phenomenological consequences.

##  V Conclusion

Hydrodynamics is an effective theory that describes the behavior of quark-
gluon plasma in a heavy-ion collision. With viscous hydrodynamics comes
thermal noise that affect the macroscopic fluid variables. In this paper, we
compared the two prevalent ways to study the evolution and effects of thermal
fluctuations in relativistic hydrodynamics. The direct stochastic calculation
is compared with the deterministic way of tracking the 2-point functions of
hydrodynamic variables which are called hydro-kinetics in the literature. Here
we considered the relatively simple case of a Bjorken evolution of a single
conserved charge, including transverse fluctuations, and focus on the
comparison of the 2-point charge correlation functions. We first perform
comparisons using white noise and find exact agreements between the two
approaches. This agreement is expected because white noise only introduces
2-point correlations. We then extend the stochastic hydrodynamics approach to
study the evolution of colored (Catteneo) noise with the causal diffusion
equation. The colored noise together with causal evolution results in
significant deviations from the results at the white noise limit.

The main result of our work is that the colored noise evolved with the causal
diffusion equation generates a peak structure around $k_{\xi}\sim k_{*}$ and
decaying tail at large $k_{\xi}$ in the 2-point correlation function in the
Fourier space as shown in Fig. 4. These features can be better understood in
coordinate space. The peak around $k_{*}$ is a consequence of the correlation
function being confined inside the light cone. The 2-point charge correlation
function at large $k_{\xi}$ is no longer a constant in the causal case because
the self-correlation has a finite extension in the Catteneo noise case. Even
if the system was non-expanding, the Catteneo noise would have induced this
feature. These features in the causal charge correlation functions could lead
to phenomenological consequences in final-state particle correlations.
Specifically, the enhancement and narrowing of the charge balance functions in
the presence of colored noise has been discussed in Ref. De  _et al._ (2020).

The quantity $k_{*}$ is an important emergent length scale in the dynamically
evolving system. Modes with $k_{\xi}\gtrsim 0.5k_{*}$ remain in equilibrium
whereas those with $k_{\xi}\lesssim 0.5k_{*}$ are pulled out of equilibrium
because of the background expansion of the system. The value of $k_{*}$ should
not change significantly in the Catteneo noise with causal diffusion. This is
because the competition between the quenching rate for the fluctuations and
the medium expansion rate is not affected a lot by the non-linearity
introduced by the colored noise and causal equation of motion. The $k_{*}$ for
the causal case will be a little less than that of the usual diffusion case.
Hence, a possible future direction would be to study observables associated
with $k_{*}$. Another interesting avenue would be to explore the emergent
scales when a critical point is introduced into the equation of state. This
question was investigated in Ref. Akamatsu  _et al._ (2019).

It has been shown in this work that the hydro-kinetic approach is equivalent
to solving the white noise SDE. It remains an open question whether the causal
evolution and correlated color noise can be implemented in the hydro-kinetic
approach. As shown in Fig. 5, the correlated Catteneo noise introduces higher-
order correlations. We also find that the causal diffusion evolution
introduces higher-order correlations because of its non-linear treatment of
noise. It would be interesting to compare the higher-point correlation
functions with the results from the hydro-kinetic approach in Ref. An  _et
al._ (2021).

Finally, it is clear from the discussion at the end of Sec. IV and Fig. 9 that
considering causal diffusion without colored noise leads to a buildup of
charge correlation in small $k_{\xi}$. The location of this buildup depends on
grid size, which is unphysical. Causal diffusion and colored noise need to be
considered as a package whenever we are considering the effects of stochastic
noise on hydrodynamic variables. In the full-fledged causal Israel-Stewart
hydrodynamics Young  _et al._ (2015); Murase (2015); Singh  _et al._ (2019);
Sakai  _et al._ (2020), stochastic fluctuations are sourced as white noises in
the relaxation-type of equation of motion for shear stress tensor to generate
the finite time correlations. Such treatment is consistent with the Catteneo
noise generation in our case. One can go further and consider the Gurtin-
Pipkin Gurtin and Pipkin (1968) noise which introduces a noise correlation in
the spatial direction in addition to the correlation in time. Gurtin-Pipkin
noise is in some sense more ‘physical’ and has been dealt with analytically in
Ref. Kapusta and Young (2014a). Exploration with a noise profile with a
spatial correlation is left to future work.

##  VI Acknowledgments

The authors thank Derek Teaney and Mayank Singh for their helpful discussions.
This work was supported by the U.S. Department of Energy, Office of Science,
Office of Nuclear Physics, under DOE Contract Nos. DE-FG02-87ER40328 and DE-
SC0021969. C.S. acknowledges a DOE Office of Science Early Career Award. We
acknowledge the Minnesota Supercomputing Institute (MSI) at the University of
Minnesota for providing resources that contributed to the research results
reported within this paper.

##  Appendix A Calculations in 1+1 D

It is instructive to understand noise propagation in 1+1 D dynamics. In this
simplified case, start with the white noise diffusion equation given by

$\frac{\partial
X}{\partial\tau}-\frac{D_{Q}}{\tau^{2}}\frac{\partial^{2}X}{\partial\xi^{2}}+s\frac{\partial
f}{\partial\xi}=0\,.$ (41)

The Fourier transform of $X$ is

$X(\xi,\tau)=\int_{-\infty}^{\infty}\frac{dk}{2\pi}e^{ik\xi}\tilde{X}(k,\tau)\,,$
(42)

and similarly for $f$. Then the SDE for white noise is

$\frac{\partial}{\partial\tau}\tilde{X}+\frac{D_{Q}k^{2}}{\tau^{2}}\tilde{X}=-iks\tilde{f}\,.$
(43)

Given a SDE of the form of Eq. (43) one can write

$\tau\delta\tilde{n}_{Q}(k,\tau)=-\int^{\tau_{f}}_{\tau_{i}}d\tau^{\prime}\,s(\tau^{\prime})\,\tilde{G}(k;\tau,\tau^{\prime})\,\tilde{f}(k,\tau^{\prime}).$
(44)

Then the 2-point correlation function can be written as

$\displaystyle\langle\tau_{1}\delta\tilde{n}_{Q}(k_{1},\tau_{1})\tau_{2}\delta\tilde{n}_{Q}(k_{2},\tau_{2})\rangle=\int
d\tau^{\prime}_{1}s(\tau^{\prime}_{1})\int\tau^{\prime}_{2}s(\tau^{\prime}_{2})\qquad$
$\displaystyle\times\tilde{G}(k_{1};\tau_{1},\tau^{\prime}_{1})\tilde{G}(k_{2};\tau_{2},\tau^{\prime}_{2})\langle\tilde{f}(k_{1},\tau^{\prime}_{1})\tilde{f}(k_{2},\tau^{\prime}_{2})\rangle.\qquad$
(45)

The Green’s function for the homogeneous part of Eq. (43) is

$\tilde{G}(k,\tau,\tau^{\prime})=ike^{D_{Q}k^{2}(\frac{1}{\tau}-\frac{1}{\tau^{\prime}})}.$
(46)

This Green’s function is for acausal white noise. After evaluating all the
integrals one gets

$N(\tau_{f},k)=\frac{\chi_{Q}T_{f}}{A\tau_{f}}\left[1-e^{-2D_{Q}k^{2}(\frac{1}{\tau_{i}}-\frac{1}{\tau_{f}})}\right]\,.$
(47)

The large $k$ limit is $\chi_{Q}T_{f}/A\tau_{f}$. We compare this to the
results of the SDE in Fig. 11.

Figure 11: (Color online) Comparison of the 2-point charge correlation
function in momentum space from the hydro-kinetics and stochastic
hydrodynamics with white noise in the 1+1D case and with $A=1$ fm2.

One can see from the figure that the agreement is excellent. This is expected
because white noise only contains 2-point functions and the hydro-kinetic
equation just tracks the effect of it.

Balancing the rate of dissipation to the rate of Bjorken expansion leads to

$D_{Q}k_{*}^{2}=\tau_{f}.$ (48)

For the parameter values used in the body of the paper, $k_{*}\approx 6.26$.
We can see from the plot that beyond $k_{*}$ the $k$-modes asymptote to the
value of $\chi_{Q}T/A\tau$.

The charge conservation equation to be solved for Catteneo noise is Kapusta
and Plumberg (2018)

$\frac{\tau_{Q}\tau}{D_{Q}\chi_{Q}T}\frac{\partial^{2}\tilde{X}}{\partial\tau^{2}}+\left[\frac{\tau}{D_{Q}\chi_{Q}T}+\tau_{Q}\frac{\partial}{\partial\tau}\left(\frac{\tau}{D_{Q}\chi_{Q}T}\right)\right]\frac{\partial\tilde{X}}{\partial\tau}.$
$+\frac{k^{2}}{\tau\chi_{Q}T}\tilde{X}=-i\frac{k\tau_{Q}\tau
s}{D_{Q}\chi_{Q}T}\frac{\partial\tilde{f}}{\partial\tau}$ $-ik\left[\frac{\tau
s}{D_{Q}\chi_{Q}T}+\tau_{Q}\frac{\partial}{\partial\tau}\left(\frac{\tau
s}{D_{Q}\chi_{Q}T}\right)\right]\tilde{f}.$ (49)

In Ref. Kapusta and Plumberg (2018) it was assumed that $\tau_{Q}$ and $D_{Q}$
were constant while $\chi_{Q}$ was proportional to $T^{2}$ and $s$ was
proportional to $T^{3}$. That allowed for the Green’s function to be
calculated and expressed in terms of special functions. One can use Eq. (45)
to compute the 2-point correlator in k-space. Alternatively, one can simulate
the SDE Eq. (49) numerically and obtain the 2-point correlator that way. The
procedure for simulating such stochastic differential equations was given in
Ref. De  _et al._ (2020). A brief summary of the numerical stochastic
machinery is provided in Appendix B.

Figure 12: (Color online) The 2-point charge correlation function in momentum
space with several values of $\tau_{Q}$ in the (1+1)D system. The plot is at
final time $\tau_{f}=6.3$ fm/$c$ with $A=1$ fm2.

A plot of $N$ for Catteneo noise for various values of $\tau_{Q}$ is shown in
Fig. 12. Additionally, we see that in the presence of transverse fluctuations,
even though the 2-point correlator for the small $k$ values differs, the large
$k$ values for the 2-point correlation for charge fluctuation are the same
with or without transverse diffusion. This is illustrated in Fig. 13.

Figure 13: (Color online) Effects of transverse fluctuations in the presence
of causal Catteneo noise on the 2-point charge correlation functions. (a) is
for $\tau_{Q}=0.1$ fm/$c$ and (b) is for $\tau_{Q}=0.4$ fm/$c$. Both sets are
computed with $10^{5}$ events.

##  Appendix B Stochastic calculus on a discrete lattice

The 2-point correlation of a function is given by

$\langle
f(x)f(x^{\prime})\rangle=\displaystyle{\frac{\displaystyle{\sum_{\mathrm{all\>random\>events}}}f(x)f(x^{\prime})\quad\quad}{\mathrm{number\>of\>random\>events}}}.$
(50)

In the continuum limit a white noise random function is defined as

$\langle f(x)f(x^{\prime})\rangle=M(x)\delta(x-x^{\prime})\text{ and }\langle
f(x)\rangle=0$ (51)

where $M(x)$ is a normalization factor and all higher-order cummulants are
required to vanish. The event-by-event distribution of $f$’s fluctuations is
therefore a normal distribution with finite variance and zero mean. One can
infer that, after discretization, we will have

$\langle f(x_{i})f(x_{i^{\prime}})\rangle=\frac{\delta_{ii^{\prime}}}{\Delta
x}.$ (52)

The normalization factor in the above equation has been set to unity. The
$\delta_{ii^{\prime}}/\Delta x$ becomes a Dirac delta function in the
continuum limit $\Delta x\rightarrow 0$. Therefore, we sample the white noise
function $f$ from a normal distribution of mean $0$ and standard deviation
$1/\sqrt{\Delta x}$. We use a random number generator for a large number of
instances (e.g. $10^{7}$) to minimize statistical noise. Subsequently, one
finds the following relations for correlation functions on a discrete lattice

$\displaystyle\langle
f(x_{i})f^{\prime}(x_{i^{\prime}})\rangle=\frac{\delta_{i+1,i^{\prime}}-\delta_{i,i^{\prime}}}{\Delta
x^{2}}\hskip 48.36958pt$ (53)

$\displaystyle\langle
f^{\prime}(x_{i})f^{\prime}(x_{i^{\prime}})\rangle=-\frac{\delta_{i,i^{\prime}+1}+\delta_{i,i^{\prime}-1}-2\delta_{i,i^{\prime}}}{\Delta
x^{3}}$ (54)

$\left\langle\int_{x_{i}}^{x_{1}}f(x^{\prime})dx^{\prime}\int_{x_{i}}^{x_{2}}f(x^{\prime})dx^{\prime}\right\rangle=\min(x_{1},x_{2})-x_{i}$
(55)

## References

  * Heinz and Snellings (2013) Ulrich Heinz and Raimond Snellings, “Collective flow and viscosity in relativistic heavy-ion collisions,” Ann. Rev. Nucl. Part. Sci. 63, 123–151 (2013), arXiv:1301.2826 [nucl-th] .
  * Gale  _et al._ (2013) Charles Gale, Sangyong Jeon, and Bjoern Schenke, “Hydrodynamic Modeling of Heavy-Ion Collisions,” Int. J. Mod. Phys. A 28, 1340011 (2013), arXiv:1301.5893 [nucl-th] .
  * Shen and Yan (2020) Chun Shen and Li Yan, “Recent development of hydrodynamic modeling in heavy-ion collisions,” Nucl. Sci. Tech. 31, 122 (2020), arXiv:2010.12377 [nucl-th] .
  * Pratt  _et al._ (2015) Scott Pratt, Evan Sangaline, Paul Sorensen, and Hui Wang, “Constraining the Eq. of State of Super-Hadronic Matter from Heavy-Ion Collisions,” Phys. Rev. Lett. 114, 202301 (2015), arXiv:1501.04042 [nucl-th] .
  * Bernhard  _et al._ (2016) Jonah E. Bernhard, J. Scott Moreland, Steffen A. Bass, Jia Liu, and Ulrich Heinz, “Applying Bayesian parameter estimation to relativistic heavy-ion collisions: simultaneous characterization of the initial state and quark-gluon plasma medium,” Phys. Rev. C 94, 024907 (2016), arXiv:1605.03954 [nucl-th] .
  * Bernhard  _et al._ (2019) Jonah E. Bernhard, J. Scott Moreland, and Steffen A. Bass, “Bayesian estimation of the specific shear and bulk viscosity of quark–gluon plasma,” Nature Phys. 15, 1113–1117 (2019).
  * Everett  _et al._ (2021a) D. Everett _et al._ (JETSCAPE), “Phenomenological constraints on the transport properties of QCD matter with data-driven model averaging,” Phys. Rev. Lett. 126, 242301 (2021a), arXiv:2010.03928 [hep-ph] .
  * Everett  _et al._ (2021b) D. Everett _et al._ (JETSCAPE), “Multisystem Bayesian constraints on the transport coefficients of QCD matter,” Phys. Rev. C 103, 054904 (2021b), arXiv:2011.01430 [hep-ph] .
  * Nijs  _et al._ (2021a) Govert Nijs, Wilke van der Schee, Umut Gürsoy, and Raimond Snellings, “Bayesian analysis of heavy ion collisions with the heavy ion computational framework Trajectum,” Phys. Rev. C 103, 054909 (2021a), arXiv:2010.15134 [nucl-th] .
  * Nijs  _et al._ (2021b) Govert Nijs, Wilke van der Schee, Umut Gürsoy, and Raimond Snellings, “Transverse Momentum Differential Global Analysis of Heavy-Ion Collisions,” Phys. Rev. Lett. 126, 202301 (2021b), arXiv:2010.15130 [nucl-th] .
  * Song  _et al._ (2011) Huichao Song, Steffen A. Bass, Ulrich Heinz, Tetsufumi Hirano, and Chun Shen, “200 A GeV Au+Au collisions serve a nearly perfect quark-gluon liquid,” Phys. Rev. Lett. 106, 192301 (2011), [Erratum: Phys.Rev.Lett. 109, 139904 (2012)], arXiv:1011.2783 [nucl-th] .
  * Shen  _et al._ (2016) Chun Shen, Zhi Qiu, Huichao Song, Jonah Bernhard, Steffen Bass, and Ulrich Heinz, “The iEBE-VISHNU code package for relativistic heavy-ion collisions,” Comput. Phys. Commun. 199, 61–85 (2016), arXiv:1409.8164 [nucl-th] .
  * Putschke  _et al._ (2019) J. H. Putschke _et al._ , “The JETSCAPE framework,” (2019), arXiv:1903.07706 [nucl-th] .
  * Schenke  _et al._ (2020) Bjoern Schenke, Chun Shen, and Prithwish Tribedy, “Running the gamut of high energy nuclear collisions,” Phys. Rev. C 102, 044905 (2020), arXiv:2005.14682 [nucl-th] .
  * Young  _et al._ (2015) C. Young, J. I. Kapusta, C. Gale, S. Jeon, and B. Schenke, “Thermally Fluctuating Second-Order Viscous Hydrodynamics and Heavy-Ion Collisions,” Phys. Rev. C 91, 044901 (2015), arXiv:1407.1077 [nucl-th] .
  * Kapusta  _et al._ (2012) J. I. Kapusta, B. Muller, and M. Stephanov, “Relativistic Theory of Hydrodynamic Fluctuations with Applications to Heavy Ion Collisions,” Phys. Rev. C 85, 054906 (2012), arXiv:1112.6405 [nucl-th] .
  * Pratt and Young (2017) Scott Pratt and Clint Young, “Relating Measurable Correlations in Heavy Ion Collisions to Bulk Properties of Equilibrated QCD Matter,” Phys. Rev. C 95, 054901 (2017), arXiv:1612.08983 [nucl-th] .
  * Pratt and Plumberg (2019) Scott Pratt and Christopher Plumberg, “Evolving Charge Correlations in a Hybrid Model with both Hydrodynamics and Hadronic Boltzmann Descriptions,” Phys. Rev. C 99, 044916 (2019), arXiv:1812.05649 [nucl-th] .
  * Plumberg and Kapusta (2017) Christopher Plumberg and Joseph I. Kapusta, “Hydrodynamic fluctuations near a critical endpoint and Hanbury-Brown–Twiss interferometry,” Phys. Rev. C 95, 044910 (2017), arXiv:1702.01368 [nucl-th] .
  * Bzdak  _et al._ (2020) Adam Bzdak, Shinichi Esumi, Volker Koch, Jinfeng Liao, Mikhail Stephanov, and Nu Xu, “Mapping the Phases of Quantum Chromodynamics with Beam Energy Scan,” Phys. Rept. 853, 1–87 (2020), arXiv:1906.00936 [nucl-th] .
  * Bluhm  _et al._ (2020) Marcus Bluhm _et al._ , “Dynamics of critical fluctuations: Theory – phenomenology – heavy-ion collisions,” Nucl. Phys. A 1003, 122016 (2020), arXiv:2001.08831 [nucl-th] .
  * An  _et al._ (2022) Xin An _et al._ , “The BEST framework for the search for the QCD critical point and the chiral magnetic effect,” Nucl. Phys. A 1017, 122343 (2022), arXiv:2108.13867 [nucl-th] .
  * Kapusta and Plumberg (2018) Joseph I. Kapusta and Christopher Plumberg, “Causal Electric Charge Diffusion and Balance Functions in Relativistic Heavy Ion Collisions,” Phys. Rev. C 97, 014906 (2018), [Erratum: Phys.Rev.C 102, 019901 (2020)], arXiv:1710.03329 [nucl-th] .
  * Pratt  _et al._ (2018) Scott Pratt, Jane Kim, and Christopher Plumberg, “Evolution of Charge Fluctuations and Correlations in the Hydrodynamic Stage of Heavy Ion Collisions,” Phys. Rev. C 98, 014904 (2018), arXiv:1712.09298 [nucl-th] .
  * Pratt and Plumberg (2020) Scott Pratt and Christopher Plumberg, “Determining the Diffusivity for Light Quarks from Experiment,” Phys. Rev. C 102, 044909 (2020), arXiv:1904.11459 [nucl-th] .
  * De  _et al._ (2020) Aritra De, Christopher Plumberg, and Joseph I. Kapusta, “Calculating Fluctuations and Self-Correlations Numerically for Causal Charge Diffusion in Relativistic Heavy-Ion Collisions,” Phys. Rev. C 102, 024905 (2020), arXiv:2003.04878 [nucl-th] .
  * Nahrgang and Bluhm (2020) Marlene Nahrgang and Marcus Bluhm, “Modeling the diffusive dynamics of critical fluctuations near the QCD critical point,” Phys. Rev. D 102, 094017 (2020), arXiv:2007.10371 [nucl-th] .
  * Rajagopal  _et al._ (2020) Krishna Rajagopal, Gregory Ridgway, Ryan Weller, and Yi Yin, “Understanding the out-of-equilibrium dynamics near a critical point in the QCD phase diagram,” Phys. Rev. D 102, 094025 (2020), arXiv:1908.08539 [hep-ph] .
  * Du  _et al._ (2020) Lipei Du, Ulrich Heinz, Krishna Rajagopal, and Yi Yin, “Fluctuation dynamics near the QCD critical point,” Phys. Rev. C 102, 054911 (2020), arXiv:2004.02719 [nucl-th] .
  * Du  _et al._ (2021) Lipei Du, Xin An, and Ulrich Heinz, “Baryon transport and the QCD critical point,” Phys. Rev. C 104, 064904 (2021), arXiv:2107.02302 [hep-ph] .
  * Murase (2015) Koichi Murase,  _Causal hydrodynamic fluctuations and their effects on high-energy nuclear collisions_ , Ph.D. thesis, Tokyo U. (2015).
  * Singh  _et al._ (2019) Mayank Singh, Chun Shen, Scott McDonald, Sangyong Jeon, and Charles Gale, “Hydrodynamic Fluctuations in Relativistic Heavy-Ion Collisions,” Nucl. Phys. A 982, 319–322 (2019), arXiv:1807.05451 [nucl-th] .
  * Sakai  _et al._ (2020) Azumi Sakai, Koichi Murase, and Tetsufumi Hirano, “Rapidity decorrelation of anisotropic flow caused by hydrodynamic fluctuations,” Phys. Rev. C 102, 064903 (2020), arXiv:2003.13496 [nucl-th] .
  * Sakai  _et al._ (2021) Azumi Sakai, Koichi Murase, and Tetsufumi Hirano, “Effects of hydrodynamic and initial longitudinal fluctuations on rapidity decorrelation of collective flow,” (2021), arXiv:2111.08963 [nucl-th] .
  * Akamatsu  _et al._ (2017) Yukinao Akamatsu, Aleksas Mazeliauskas, and Derek Teaney, “A kinetic regime of hydrodynamic fluctuations and long time tails for a Bjorken expansion,” Phys. Rev. C 95, 014909 (2017), arXiv:1606.07742 [nucl-th] .
  * Stephanov and Yin (2018) M. Stephanov and Y. Yin, “Hydrodynamics with parametric slowing down and fluctuations near the critical point,” Phys. Rev. D 98, 036006 (2018), arXiv:1712.10305 [nucl-th] .
  * Akamatsu  _et al._ (2019) Yukinao Akamatsu, Derek Teaney, Fanglida Yan, and Yi Yin, “Transits of the QCD critical point,” Phys. Rev. C 100, 044901 (2019), arXiv:1811.05081 [nucl-th] .
  * An  _et al._ (2019) Xin An, Gokce Basar, Mikhail Stephanov, and Ho-Ung Yee, “Relativistic Hydrodynamic Fluctuations,” Phys. Rev. C 100, 024910 (2019), arXiv:1902.09517 [hep-th] .
  * Martinez  _et al._ (2019) M. Martinez, T. Schäfer, and V. Skokov, “Critical behavior of the bulk viscosity in QCD,” Phys. Rev. D 100, 074017 (2019), arXiv:1906.11306 [hep-ph] .
  * An  _et al._ (2020) Xin An, Gökçe Başar, Mikhail Stephanov, and Ho-Ung Yee, “Fluctuation dynamics in a relativistic fluid with a critical point,” Phys. Rev. C 102, 034901 (2020), arXiv:1912.13456 [hep-th] .
  * Martinez and Schäfer (2017) Mauricio Martinez and Thomas Schäfer, “Hydrodynamic tails and a fluctuation bound on the bulk viscosity,” Phys. Rev. A 96, 063607 (2017), arXiv:1708.01548 [cond-mat.quant-gas] .
  * Bjorken (1983) J. D. Bjorken, “Highly Relativistic Nucleus-Nucleus Collisions: The Central Rapidity Region,” Phys. Rev. D 27, 140–151 (1983).
  * Kapusta and Young (2014a) J. I. Kapusta and C. Young, “Causal Baryon Diffusion and Colored Noise,” Phys. Rev. C 90, 044902 (2014a), arXiv:1404.4894 [nucl-th] .
  * Kapusta and Young (2014b) J. I. Kapusta and C. Young, “Causal Baryon Diffusion and Colored Noise for Heavy Ion Collisions,” Nucl. Phys. A 931, 1051–1055 (2014b), arXiv:1409.1957 [nucl-th] .
  * Cattaneo  _et al._ (1958) C. Cattaneo, J.K. de Fériet, and Académie des sciences (France),  _Sur une forme de l’équation de la chaleur éliminant le paradoxe d’une propagation instantanée_, Comptes rendus hebdomadaires des séances de l’Académie des sciences (Gauthier-Villars, 1958).
  * An  _et al._ (2021) Xin An, Gökçe Başar, Mikhail Stephanov, and Ho-Ung Yee, “Evolution of Non-Gaussian Hydrodynamic Fluctuations,” Phys. Rev. Lett. 127, 072301 (2021), arXiv:2009.10742 [hep-th] .
  * Gurtin and Pipkin (1968) M.E. Gurtin and A.C. Pipkin, “A general theory of heat conduction with finite wave speeds.” Arch. Rational Mech. Anal. 31, 113–126 (1968).