# Bayesian Linear Models: A compact general set of results.

J. Andrés Christen
Centro de Investigación en Matemáticas (CIMAT), CONACYT,
Guanajuato, GTO, Mexico<EMAIL_ADDRESS>
(Nov 2023)

###### Abstract

I present all the details in calculating the posterior distribution of the
conjugate Normal-Gamma prior in Bayesian Linear Models (BLM), including
correlated observations, prediction, model selection and comments on efficient
numeric implementations. A Python implementation is also presented. These have
been presented and available in many books and texts but, I believe, a general
compact and simple presentation is always welcome and not always simple to
find. Since correlated observations are also included, these results may also
be useful for time series analysis and spacial statistics. Other particular
cases presented include regression, Gaussian processes and Bayesian Dynamic
Models

##  1 Preliminaries

This paper is not intended to explain the details of Bayesian regression
models [1], fitting Gaussian processes [5] [4], Bayesian time series (dynamic
models) [6] and Bayessian linear models in general. The cited books should
served as a starting point to learn about each topic, and others. However,
very many of the models involved in those topics, and others related, are
based on the same basic calculations, namely, the derivation of the Normal-
Gamma conjugate posterior distribution, the Normal and t predictive and the
normalization constants for model comparisons, in the Bayesian analysis of
linear models. This paper presents, in detail, the basic algebra involved in
such derivations. These have been presented and available in many books and
texts, as already said, but, I believe, a general compact and simple
presentation is always welcome and not always simple to find.

Moreover, I include some comments on the efficient numeric evaluation for the
required hyper parameters and a compact Python class with an implementation
https://github.com/andreschristen/BLM. I also include some examples as an
illustration.

I include some preliminary results that will be required:

  * •

The Sherman-Woodbury-Morrison lemma
111https://en.wikipedia.org/wiki/Woodbury_matrix_identity. Let
$\boldsymbol{E}$ be a $n\times n$ and $\boldsymbol{D}$ a $k\times k$
invertible matrices and $\boldsymbol{B}$ $n\times k$ and $\boldsymbol{C}$
$k\times n$ matrices. Then

$(\boldsymbol{E}+\boldsymbol{B}\boldsymbol{D}\boldsymbol{C})^{-1}=\boldsymbol{E}^{-1}-\boldsymbol{E}^{-1}\boldsymbol{B}(\boldsymbol{D}^{-1}+\boldsymbol{C}\boldsymbol{E}^{-1}\boldsymbol{B})^{-1}\boldsymbol{C}\boldsymbol{E}^{-1}.$
(1)

  * •

The matrix determinant
lemma222https://en.wikipedia.org/wiki/Matrix_determinant_lemma. Let
$\boldsymbol{E}$ be a $n\times n$ invertible matrix and $\boldsymbol{u}$ and
$\boldsymbol{v}$ $n\times 1$ vectors, then

$\left|\boldsymbol{E}+\boldsymbol{u}\boldsymbol{v}^{\prime}\right|=\left|\boldsymbol{E}\right|(1+\boldsymbol{v}^{\prime}\boldsymbol{E}^{-1}\boldsymbol{u}).$
(2)

  * •

The multivariate Normal
distribution333https://en.wikipedia.org/wiki/Multivariate_normal_distribution
[2]. If the ($n\times 1$ column) vector
$\boldsymbol{Y}\sim\mathrm{N}_{n}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ then:

    1. 1.

The pdf of $\boldsymbol{Y}$ is

$f(\boldsymbol{y})=(2\pi)^{-\frac{n}{2}}\left|\boldsymbol{\Sigma}\right|^{-\frac{1}{2}}\exp\left\\{-\frac{1}{2}(\boldsymbol{y}-\boldsymbol{\mu})^{\prime}\boldsymbol{\Sigma}^{-1}(\boldsymbol{y}-\boldsymbol{\mu})\right\\}.$

$V(\boldsymbol{Y})=\boldsymbol{\Sigma}$ is the positive definite covariance
matrix and $E(\boldsymbol{Y})=\boldsymbol{\mu}$ the mean vector.

    2. 2.

Marginally, each entry of $\boldsymbol{Y}$ is univariate normal
$Y_{i}\sim\mathrm{N}(\mu_{i},\sigma_{i,i})$.

    3. 3.

If $\boldsymbol{Y}=[\boldsymbol{Y}_{1},\boldsymbol{Y}_{2}]^{\prime}$,
$\boldsymbol{Y}_{1},q\times 1$ and $\boldsymbol{Y}_{2},(n-q)\times 1$
subvectors, the conditional is

$\boldsymbol{Y}_{1}\mid\boldsymbol{Y}_{2}=\boldsymbol{y}_{2}\sim\mathrm{N}_{q}(\boldsymbol{\mu}_{1}+\boldsymbol{\Sigma}_{12}\boldsymbol{\Sigma}_{22}^{-1}(\boldsymbol{y}_{2}-\boldsymbol{\mu}_{2}),\boldsymbol{\Sigma}_{11}-\boldsymbol{\Sigma}_{12}\boldsymbol{\Sigma}_{22}^{-1}\boldsymbol{\Sigma}_{12}).$

    4. 4.

If $\boldsymbol{B}$ is a $q\times n$ matrix and $\boldsymbol{a}$ a $q\times 1$
vector then

$\boldsymbol{B}\boldsymbol{Y}+\boldsymbol{a}\sim\mathrm{N}_{q}(\boldsymbol{B}\boldsymbol{\mu}+\boldsymbol{a},\boldsymbol{B}\boldsymbol{\Sigma}\boldsymbol{B}^{\prime}).$

  * •

For the multivariate t
distribution444https://en.wikipedia.org/wiki/Multivariate_t-distribution, if
($n\times 1$ column) vector $\boldsymbol{Y}\sim
T_{n}^{\nu}(\boldsymbol{\mu},\boldsymbol{D})$ then its pdf is

$f(\boldsymbol{y})=\pi^{-\frac{n}{2}}\left|\nu\boldsymbol{D}\right|^{-\frac{1}{2}}\frac{\Gamma\left(\frac{\nu+n}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right)}\left[1+\frac{(\boldsymbol{y}-\boldsymbol{\mu})^{\prime}\boldsymbol{D}^{-1}(\boldsymbol{y}-\boldsymbol{\mu})}{\nu}\right]^{-\left(\frac{\nu+n}{2}\right)},$

$\nu$ are the degrees of freedom, $E(\boldsymbol{Y})=\boldsymbol{\mu}$ is the
mean vector, $\boldsymbol{D}$ is the dispersion matrix, and
$V(\boldsymbol{Y})=\frac{\nu}{\nu-2}\boldsymbol{D}$ if $\nu>2$, otherwise
$\boldsymbol{Y}$ does not have second moments.

  * •

If $Y\sim Ga(\alpha,\beta)$ (the Gamma distribution), its pdf, in our
preferred parametrization, is

$f(y)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}y^{\alpha-1}\exp\\{-y\beta\\}I_{(0,\infty)}(y).$

$\alpha$ is the shape parameter and $\beta$ is the rate parameter. Here
$E(Y)=\frac{\alpha}{\beta}$ and $V(Y)=\frac{\alpha}{\beta^{2}}$.

  * •

Finnaly, from the Gamma pdf one can see that

$\int_{0}^{\infty}y^{\alpha-1}\exp\\{-y\beta\\}dy=\beta^{-\alpha}\Gamma(\alpha).$

##  2 The model, the prior and the posterior

The linear model is

$\boldsymbol{Y}=\boldsymbol{X}\boldsymbol{\theta}+\boldsymbol{\epsilon};\boldsymbol{\epsilon}\sim\mathrm{N}_{n}(\boldsymbol{0},\lambda^{-1}\boldsymbol{\Sigma}).$
(3)

Where:

  * •

$\boldsymbol{Y}$ is a $n\times 1$ vector of responses.

  * •

$\boldsymbol{X}$ is a $n\times p$ design or covariate matrix

  * •

$\boldsymbol{\theta}$ is a $p\times 1$ vector of unknown parameters.

  * •

$\boldsymbol{\Sigma}$ is a $n\times n$ variance-covariance matrix, and will be
taken as known. The general case of unknown covariance is not considered here.
We let $\boldsymbol{A}=\boldsymbol{\Sigma}^{-1}$, the $n\times n$ precision
matrix.

  * •

$\lambda$ is the precision parameter, will be taken unknown and the case of
known $\lambda$ will be treated as a special case. The usual case is that
$\boldsymbol{\Sigma}$ is a correlation matrix and
$V(\boldsymbol{Y})=\lambda^{-1}$ is the common variance for all responses.
However, $\boldsymbol{\Sigma}$ may be any symmetric positive definite matrix,
and $\lambda$ is a “deflation” parameter.

Note: In the Python `pyblm` implementation I tried to use the similar names
for variables, e.g. `X` for $\boldsymbol{X}$, `th` for $\boldsymbol{\theta}$,
`la` for $\lambda$ etc.

The conjugate prior is a Normal-Gamma, that is, conditional on $\lambda$,
$\boldsymbol{\theta}$ is
$\mathrm{N}(\boldsymbol{\theta}_{0},\boldsymbol{A}_{0}^{-1})$ and $\lambda$
has a $Ga(\alpha_{0},\beta_{0})$ distribution. Namely

$\pi(\boldsymbol{\theta}|\lambda)=(2\pi)^{-\frac{p}{2}}\lambda^{\frac{p}{2}}\left|\boldsymbol{A}_{0}\right|^{\frac{1}{2}}\exp\left\\{-\frac{\lambda}{2}(\boldsymbol{\theta}-\boldsymbol{\theta}_{0})^{\prime}\boldsymbol{A}_{0}(\boldsymbol{\theta}-\boldsymbol{\theta}_{0})\right\\}$

and

$\pi(\lambda)=\frac{\beta_{0}^{\alpha_{0}}}{\Gamma(\alpha_{0})}\lambda^{\alpha_{0}-1}\exp\left\\{-\lambda\beta_{0}\right\\}I_{[0,\infty)}(\lambda),$

  * •

$\boldsymbol{\theta}_{0}$ is a $p\times 1$ mean vector for the prior.

  * •

$\boldsymbol{A}_{0}$ is a $p\times p$ precision matrix for the prior.

  * •

$\alpha_{0}>0$ and $\beta_{0}>0$ are the shape and rate parameters of the
Gamma prior for $\lambda$. In this parametrization the expected value is
$\frac{\alpha_{0}}{\beta_{0}}$.

The model is this, for some responses $\boldsymbol{Y}=\boldsymbol{y}$

$f(\boldsymbol{y}|\boldsymbol{\theta},\lambda)=(2\pi)^{-\frac{n}{2}}\lambda^{\frac{n}{2}}\left|\boldsymbol{A}\right|^{\frac{1}{2}}\exp\left\\{-\frac{\lambda}{2}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})^{\prime}\boldsymbol{A}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})\right\\}$

and the posterior is this

$\pi(\boldsymbol{\theta},\lambda|\boldsymbol{y})=\frac{f(\boldsymbol{y}|\boldsymbol{\theta},\lambda)\pi(\boldsymbol{\theta}|\lambda)\pi(\lambda)}{f(\boldsymbol{y})}$
(4)

or, copying everything,

$\displaystyle\pi(\boldsymbol{\theta},\lambda|\boldsymbol{y})=$ $\displaystyle
K(2\pi)^{-\frac{n+p}{2}}\left|\boldsymbol{A}\right|^{\frac{1}{2}}\left|\boldsymbol{A}_{0}\right|^{\frac{1}{2}}\frac{\beta_{0}^{\alpha_{0}}}{\Gamma(\alpha_{0})}\lambda^{\frac{n}{2}}\lambda^{\frac{p}{2}}\lambda^{\alpha_{0}-1}$
(5)
$\displaystyle\exp\left\\{-\lambda\left[\beta_{0}~{}+~{}\frac{1}{2}\left((\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})^{\prime}\boldsymbol{A}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta})~{}+~{}(\boldsymbol{\theta}-\boldsymbol{\theta}_{0})^{\prime}\boldsymbol{A}_{0}(\boldsymbol{\theta}-\boldsymbol{\theta}_{0})\right)\right]\right\\},$
(6)

where $K=f(y)^{-1}$, the normalizing constant. The term in the exponential is
also a quadratic form, with some constant, and therefore we should be able to
write the posterior as

$\pi(\boldsymbol{\theta},\lambda|\boldsymbol{y})\propto\lambda^{\frac{p}{2}}\exp\left\\{-\frac{\lambda}{2}(\boldsymbol{\theta}-\boldsymbol{\theta}_{n})^{\prime}\boldsymbol{A}_{n}(\boldsymbol{\theta}-\boldsymbol{\theta}_{n})\right\\}\lambda^{\alpha_{n}-1}\exp\left\\{-\lambda\beta_{n}\right\\}$

for some mean $\boldsymbol{\theta}_{n}$ vector, matrix $\boldsymbol{A}_{n}$
and parameters $\alpha_{n}=\alpha_{0}+\frac{n}{2}$ and $\beta_{n}$. If
$\boldsymbol{A}_{n}$ is positive definite and $\beta_{n}$ positive then this
is also a Normal-Gamma distribution with parameters
$\boldsymbol{\theta}_{n},\boldsymbol{A}_{n},\alpha_{n},\beta_{n}$.

##  3 The calculations

In this part I prefer to do the calculations by small steps. From (6) we
expand the two quadratic forms

$Q=\boldsymbol{y}^{\prime}\boldsymbol{A}\boldsymbol{y}~{}-~{}2\boldsymbol{y}^{\prime}\boldsymbol{A}X\boldsymbol{\theta}~{}+~{}\boldsymbol{\theta}^{\prime}\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{X}\boldsymbol{\theta}~{}+~{}\boldsymbol{\theta}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}~{}-~{}2\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}~{}+~{}\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}.$

which we need to turn into

$(\boldsymbol{\theta}-\boldsymbol{\theta}_{n})^{\prime}\boldsymbol{A}_{n}(\boldsymbol{\theta}-\boldsymbol{\theta}_{n})+c=\boldsymbol{\theta}^{\prime}\boldsymbol{A}_{n}\boldsymbol{\theta}~{}-~{}2\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{n}\boldsymbol{\theta}~{}+~{}\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{n}\boldsymbol{\theta}_{n}+c.$

  1. 1.

From terms 3 and 4 in $Q$,
$\boldsymbol{\theta}^{\prime}(\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{X}+\boldsymbol{A}_{0})\boldsymbol{\theta}=\boldsymbol{\theta}^{\prime}\boldsymbol{A}_{n}\boldsymbol{\theta}$.
Then

$\boldsymbol{A}_{n}=\boldsymbol{A}_{0}+\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{X}.$

If $\boldsymbol{X}$ is of full rank then
$\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{X}$ is positive definite and
since $\boldsymbol{A}_{0}$ is also positive definite the $\boldsymbol{A}_{n}$
is. However, even if $\boldsymbol{X}$ is not full rank, $\boldsymbol{A}_{n}$
may become positive definite. For example, a design matrix $\boldsymbol{X}$
with colinearity may still be analyzed, using an adequate prior.

  2. 2.

From terms 2 and 5 in $Q$,
$~{}-~{}2(\boldsymbol{y}^{\prime}\boldsymbol{A}\boldsymbol{X}+\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0})\boldsymbol{\theta}=-2\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{n}\boldsymbol{\theta}$.
Then

$\boldsymbol{A}_{n}\boldsymbol{\theta}_{n}=\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}+\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{y}~{}~{}\text{or}~{}~{}\boldsymbol{\theta}_{n}=\boldsymbol{A}_{n}^{-1}(\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}+\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{y}).$
(7)

  3. 3.

From terms 1 and 6 in $Q$,
$c=\boldsymbol{y}^{\prime}\boldsymbol{A}\boldsymbol{y}+\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{n}\boldsymbol{\theta}_{n}$.
Then (6) is

$\exp\left\\{-\lambda\left[\beta_{0}+\frac{1}{2}\left(Q+c\right)\right]\right\\}$

which means that

$\beta_{n}=\beta_{0}+\frac{1}{2}\left(\boldsymbol{y}^{\prime}A\boldsymbol{y}+\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{n}\boldsymbol{\theta}_{n}\right),$

which we need to prove it is always positive. See Section 3.1.

Following [4], for $\boldsymbol{A}=\boldsymbol{I}$, let
$\hat{\boldsymbol{\theta}}=(\boldsymbol{X}^{\prime}\boldsymbol{X})^{-1}\boldsymbol{X}^{\prime}\boldsymbol{y}$,
the classic MLE estimator for $\boldsymbol{\theta}$ (or the More-Penrose
inverse) and

$\boldsymbol{W}=\boldsymbol{A}_{n}^{-1}\boldsymbol{X}^{\prime}\boldsymbol{X}~{}~{}\text{and}~{}~{}\boldsymbol{I}-\boldsymbol{W}=\boldsymbol{A}_{n}^{-1}(\boldsymbol{A}_{n}-\boldsymbol{X}^{\prime}\boldsymbol{X})=\boldsymbol{A}_{n}^{-1}\boldsymbol{A}_{0},$

then

$\boldsymbol{\theta}_{n}=(\boldsymbol{I}-\boldsymbol{W})\boldsymbol{\theta}_{0}+\boldsymbol{W}\hat{\boldsymbol{\theta}}.$

That is, the posterior mean $\boldsymbol{\theta}_{n}$ is a weighted average
between the prior mean $\boldsymbol{\theta}_{0}$ and
$\hat{\boldsymbol{\theta}}$.

###  3.1 Alternative expressions for $\beta_{n}$

A very helpful shortcut for calculations is this. Since

$\displaystyle\boldsymbol{Y}=$
$\displaystyle\boldsymbol{X}\boldsymbol{\theta}+\boldsymbol{\epsilon};~{}\boldsymbol{\epsilon}\sim\mathrm{N}_{n}(\boldsymbol{0},\lambda^{-1}\boldsymbol{A}^{-1})$
$\displaystyle\boldsymbol{\theta}=$
$\displaystyle\boldsymbol{\theta}_{0}+\boldsymbol{\epsilon}_{1};~{}\boldsymbol{\epsilon}_{1}\sim\mathrm{N}_{p}(\boldsymbol{0},\lambda^{-1}\boldsymbol{A}_{0}^{-1}),$

then, using basic properties of the Normal,
$\boldsymbol{X}\boldsymbol{\theta}$ is Normal distribution and

$\boldsymbol{Y}=\boldsymbol{X}\boldsymbol{\theta}_{0}+\boldsymbol{X}\boldsymbol{\epsilon}_{1}+\boldsymbol{\epsilon}$

and

$\boldsymbol{Y}\mid\lambda\sim\mathrm{N}_{n}(\boldsymbol{X}\boldsymbol{\theta}_{0},\lambda^{-1}(\boldsymbol{A}^{-1}+\boldsymbol{X}\boldsymbol{A}_{0}^{-1}\boldsymbol{X}^{\prime})).$
(8)

The dependency on $\boldsymbol{\theta}$ has been removed. The posterior
marginal may be calculated directly, since

$\displaystyle\pi(\lambda|\boldsymbol{y})\propto$ $\displaystyle
f(\boldsymbol{y}|\lambda)\pi(\lambda)$ $\displaystyle\propto$
$\displaystyle(2\pi)^{-\frac{n}{2}}\left|\boldsymbol{A}^{-1}+\boldsymbol{X}\boldsymbol{A}_{0}^{-1}\boldsymbol{X}^{\prime}\right|^{-\frac{1}{2}}\frac{\beta_{0}^{\alpha_{0}}}{\Gamma(\alpha_{0})}\lambda^{\alpha_{0}+\frac{n}{2}-1}$
(9)
$\displaystyle\exp\left\\{-\lambda\left[\beta_{0}~{}+~{}\frac{1}{2}\left((\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{0})^{\prime}(\boldsymbol{A}^{-1}+\boldsymbol{X}\boldsymbol{A}_{0}^{-1}\boldsymbol{X}^{\prime})^{-1}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{0})\right)\right]\right\\}.$

This means that, as already seen,

$\lambda|\boldsymbol{y}\sim Ga(\alpha_{n},\beta_{n})$

with $\alpha_{n}=\alpha_{0}+\frac{n}{2}$ as above and

$\beta_{n}=\beta_{0}~{}+~{}\frac{1}{2}\left((\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{0})^{\prime}(\boldsymbol{A}^{-1}+\boldsymbol{X}\boldsymbol{A}_{0}^{-1}\boldsymbol{X}^{\prime})^{-1}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{0})\right),$
(10)

and we have proved that
$\beta_{n}=\beta_{0}+\frac{1}{2}q(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\theta}_{n},\boldsymbol{A}_{0},\boldsymbol{\theta}_{0})$
where

$\displaystyle
q(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\theta}_{n},\boldsymbol{A}_{0},\boldsymbol{\theta}_{0})=$
$\displaystyle~{}\boldsymbol{y}^{\prime}A\boldsymbol{y}+\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{n}\boldsymbol{\theta}_{n}$
$\displaystyle=$
$\displaystyle(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{0})^{\prime}(\boldsymbol{A}^{-1}+\boldsymbol{X}\boldsymbol{A}_{0}^{-1}\boldsymbol{X}^{\prime})^{-1}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{0}).$

This hints at that $\beta_{n}$ is indeed positive, if
$(\boldsymbol{A}^{-1}+\boldsymbol{X}\boldsymbol{A}_{0}^{-1}\boldsymbol{X}^{\prime})^{-1}$
is positive definite. The equivalence between the above two expressions may be
established by shear matrix algebra and using the Sherman-Woodbury-Morrison
lemma, but the equivalence is already proven above.

A third very useful expression for $\beta_{n}$ is found in the following way.
Note from (7) that
$\boldsymbol{A}_{n}\boldsymbol{\theta}_{n}=\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}+\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{y}$
then

$\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{n}\boldsymbol{\theta}_{n}=\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}+\boldsymbol{y}^{\prime}\boldsymbol{A}\boldsymbol{X}\boldsymbol{\theta}_{n}.$

Then

$\displaystyle
q(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\theta}_{n},\boldsymbol{A}_{0},\boldsymbol{\theta}_{0})=$
$\displaystyle~{}\boldsymbol{y}^{\prime}A\boldsymbol{y}+\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{n}\boldsymbol{\theta}_{n}$
$\displaystyle=$
$\displaystyle~{}\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}+\boldsymbol{y}^{\prime}A\boldsymbol{y}-\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}-\boldsymbol{y}^{\prime}\boldsymbol{A}\boldsymbol{X}\boldsymbol{\theta}_{n}$
$\displaystyle-\boldsymbol{y}^{\prime}\boldsymbol{A}\boldsymbol{X}\boldsymbol{\theta}_{n}+\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{X}\boldsymbol{\theta}_{n}+\boldsymbol{y}^{\prime}\boldsymbol{A}\boldsymbol{X}\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{X}\boldsymbol{\theta}_{n}$
$\displaystyle=$
$\displaystyle~{}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{n})^{\prime}\boldsymbol{A}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{n})~{}+$
$\displaystyle\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}+\boldsymbol{y}^{\prime}\boldsymbol{A}\boldsymbol{X}\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{X}\boldsymbol{\theta}_{n}$

For the second term we have

$\displaystyle\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}+\boldsymbol{y}^{\prime}\boldsymbol{A}\boldsymbol{X}\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{X}\boldsymbol{\theta}_{n}$
$\displaystyle=$
$\displaystyle\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{n}+\boldsymbol{y}^{\prime}\boldsymbol{A}\boldsymbol{X}\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{X}\boldsymbol{\theta}_{n}$
$\displaystyle=$
$\displaystyle\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{n}+\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{n}\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{X}\boldsymbol{\theta}_{n}$
$\displaystyle=$
$\displaystyle\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}-2\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{n}+\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{n}+\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{n}\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{X}\boldsymbol{\theta}_{n}$
$\displaystyle=$
$\displaystyle(\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n})^{\prime}\boldsymbol{A}_{0}(\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n})+\boldsymbol{\theta}_{n}^{\prime}(\boldsymbol{A}_{n}-(\boldsymbol{A}_{0}+\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{X}))\boldsymbol{\theta}_{n},$

where indeed the second term in the last equality is zero since
$\boldsymbol{A}_{n}=\boldsymbol{A}_{0}+\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{X}$.
Therefore

$q(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\theta}_{n},\boldsymbol{A}_{0},\boldsymbol{\theta}_{0})=(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{n})^{\prime}\boldsymbol{A}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{n})+(\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n})^{\prime}\boldsymbol{A}_{0}(\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n})$
(11)

or

$\beta_{n}=\beta_{0}~{}+~{}\frac{1}{2}\left[(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{n})^{\prime}\boldsymbol{A}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{n})+(\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n})^{\prime}\boldsymbol{A}_{0}(\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n})\right],$

which is indeed positive. More importantly,
$(y-\boldsymbol{X}\boldsymbol{\theta}_{n})^{\prime}\boldsymbol{A}(y-\boldsymbol{X}\boldsymbol{\theta}_{n})$
is the (weighted) sum of the squared residuals and
$(\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n})^{\prime}\boldsymbol{A}_{0}(\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n})$
is a kind of penalization for the discrepancy between the prior and posterior
means. Note that

$(E[\lambda|\boldsymbol{y}])^{-1}=\frac{\beta_{n}}{\alpha_{n}}=\frac{\beta_{0}~{}+~{}\frac{1}{2}\left[(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{n})^{\prime}\boldsymbol{A}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{n})+(\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n})^{\prime}\boldsymbol{A}_{0}(\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n})\right]}{\alpha_{0}+\frac{n}{2}}.$

##  4 The normalizing constant, model evidence and the posterior predictive
and marginal distributions

From (4) the normalizing constant, which, as a function of $\boldsymbol{y}$ is
the prior predictive model, is

$f(\boldsymbol{y})=\int f(\boldsymbol{y}|\lambda)\pi(\lambda)d\lambda.$

From (9) we have

$\displaystyle f(\boldsymbol{y})=$
$\displaystyle(2\pi)^{-\frac{n}{2}}\left|\boldsymbol{A}^{-1}+\boldsymbol{X}\boldsymbol{A}_{0}^{-1}\boldsymbol{X}^{\prime}\right|^{-\frac{1}{2}}\frac{\beta_{0}^{\alpha_{0}}}{\Gamma(\alpha_{0})}\int\lambda^{\alpha_{n}-1}\exp\left\\{-\lambda\beta_{n}\right\\}d\lambda$
$\displaystyle=$
$\displaystyle(2\pi)^{-\frac{n}{2}}|B_{0}|^{\frac{1}{2}}\frac{\beta_{0}^{\alpha_{0}}}{\Gamma(\alpha_{0})}\Gamma(\alpha_{n})\beta_{n}^{-\alpha_{n}}.$
(12)

where
$\boldsymbol{B}_{0}=(\boldsymbol{A}^{-1}+\boldsymbol{X}\boldsymbol{A}_{0}^{-1}\boldsymbol{X}^{\prime})^{-1}$.
Using the expression for $\beta_{n}$ in (10) it is simple to see that

$f(\boldsymbol{y})=\pi^{-\frac{n}{2}}\left|\nu_{0}^{-1}\boldsymbol{B}_{0}^{*}\right|^{\frac{1}{2}}\frac{\Gamma\left(\frac{\nu_{0}+n}{2}\right)}{\Gamma\left(\frac{\nu_{0}}{2}\right)}\left[1+\frac{(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{0})^{\prime}\boldsymbol{B}_{0}^{*}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{0})}{\nu_{0}}\right]^{-\left(\frac{\nu_{0}+n}{2}\right)}$

where $\nu_{0}=2\alpha_{0}$ and
$\boldsymbol{B}_{0}^{*}=\frac{\alpha_{0}}{\beta_{0}}\boldsymbol{B}_{0}$. This
is a non-central multivariate t distribution and therefore, a priori,

$\boldsymbol{Y}\sim
T_{n}^{\nu_{0}}(\boldsymbol{X}\boldsymbol{\theta}_{0},(\boldsymbol{B}_{0}^{*})^{-1}).$

For fixed $\boldsymbol{y}$ the more compact form of the normalizing constant
in (12) should be used.

###  4.1 Model comparisons, variable selection

In terms of model comparisons, one would have different design matrices
$\boldsymbol{X}_{i}$ and possible priors for $\boldsymbol{\theta}$, of
different dimensions, leading to the normalization constants for each model
$f(\boldsymbol{y}|\boldsymbol{X}_{i})$; these are called the “model evidence”
since the posterior probability of model $i$ is [3]

$\frac{f(\boldsymbol{y}|\boldsymbol{X}_{i})}{\sum_{j=1}^{M}f(\boldsymbol{y}|\boldsymbol{X}_{j})}.$

Assuming that the prior for $\lambda$ remains the same, the model evidence may
be taken as

$z(\boldsymbol{y}|\boldsymbol{X}_{i})=|\boldsymbol{B}_{0}^{i}|^{\frac{1}{2}}(\beta_{n}^{i})^{-\alpha_{n}},$
(13)

for the obvious definition of $\beta_{n}^{i}$ and
$\boldsymbol{B}_{0}^{i}=(\boldsymbol{A}^{-1}+\boldsymbol{X}_{i}\boldsymbol{A}_{0}^{-1}\boldsymbol{X}_{i}^{\prime})^{-1}$.
It should also be clear that when $\lambda$ is constant we have

$z(\boldsymbol{y}|\boldsymbol{X}_{i},\lambda)=|\boldsymbol{B}_{0}^{i}|^{\frac{1}{2}}e^{-\lambda\beta_{n}^{i}}.$

Using the Sherman-Woodbury-Morrison lemma it is easy to see that

$\boldsymbol{B}_{0}^{i}=\boldsymbol{A}-\boldsymbol{A}\boldsymbol{X}_{i}(\boldsymbol{A}_{n}^{i})^{-1}\boldsymbol{X}_{i}^{\prime}\boldsymbol{A},$

where indeed
$\boldsymbol{A}_{n}^{i}=\boldsymbol{A}_{0}+\boldsymbol{X}_{i}^{\prime}\boldsymbol{A}\boldsymbol{X}_{i}$.
When we have independent data, i.e. $\boldsymbol{A}=\boldsymbol{I}$,

$\boldsymbol{B}_{0}^{i}=\boldsymbol{I}-\boldsymbol{X}_{i}(\boldsymbol{A}_{n}^{i})^{-1}\boldsymbol{X}_{i}^{\prime},$

###  4.2 The predictive distribution

Suppose we want to predict $\boldsymbol{Z}$, a $m\times 1$ vector of
unobserved responses, corresponding to a $m\times p$ matrix of covariates
$\boldsymbol{X}_{m}$, given the data $\boldsymbol{y}$. It is more complicated
(and interesting) with correlated data. Using properties of the multivariate
Normal, and considering the joint normal distribution of
$[\boldsymbol{Y},\boldsymbol{Z}]^{\prime}$, conditional on
$\boldsymbol{Y}=\boldsymbol{y}$ we have

$\boldsymbol{Z}|\boldsymbol{\theta},\boldsymbol{Y}=\boldsymbol{y},\lambda\sim\mathrm{N}_{m}(\boldsymbol{X}_{m}\boldsymbol{\theta}+\boldsymbol{v}_{m}^{\prime}\boldsymbol{A}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}),\lambda^{-1}((\boldsymbol{A}^{m})^{-1}-\boldsymbol{v}_{m}^{\prime}\boldsymbol{A}\boldsymbol{v}_{m}))$

where $V(\boldsymbol{Z})=\lambda^{-1}(\boldsymbol{A}^{m})^{-1}$ the $m\times
m$ variance-covariance matrix of $\boldsymbol{Z}$ and
$cov(\boldsymbol{Y},\boldsymbol{Z})=\lambda^{-1}\boldsymbol{v}_{m}$ the
$n\times m$ matrix of cross covariances of
$[\boldsymbol{Y},\boldsymbol{Z}]^{\prime}$.

Indeed, we already know that

$\boldsymbol{\theta}|\boldsymbol{Y}=\boldsymbol{y},\lambda\sim\mathrm{N}_{p}(\boldsymbol{\theta}_{n},\lambda^{-1}\boldsymbol{A}_{n}^{-1})$

and therefore

$\displaystyle\boldsymbol{Z}|\boldsymbol{Y}=\boldsymbol{y},\lambda\sim\mathrm{N}_{m}(\boldsymbol{X}_{m}\boldsymbol{\theta}_{n}+\boldsymbol{v}_{m}^{\prime}\boldsymbol{A}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{n}),$
$\displaystyle\lambda^{-1}((\boldsymbol{A}^{m})^{-1}-\boldsymbol{v}_{m}^{\prime}\boldsymbol{A}\boldsymbol{v}_{m}+(\boldsymbol{X}_{m}-\boldsymbol{v}_{m}^{\prime}\boldsymbol{A}\boldsymbol{X})\boldsymbol{A}_{n}^{-1}(\boldsymbol{X}_{m}-\boldsymbol{v}_{m}^{\prime}\boldsymbol{A}\boldsymbol{X})^{\prime})).$

(For the conditional mean of $\boldsymbol{Z}$ we use the alternative
expression
$\boldsymbol{v}_{m}^{\prime}\boldsymbol{A}\boldsymbol{y}+(\boldsymbol{X}_{m}-\boldsymbol{v}_{m}^{\prime}\boldsymbol{A}\boldsymbol{X})\boldsymbol{\theta}$
to calculate the variance-covariance matrix.) Note that for independent data
$\boldsymbol{A}=\boldsymbol{I},\boldsymbol{A}_{m}=\boldsymbol{I}$ and
$\boldsymbol{v}_{m}=\boldsymbol{0}$ and the expression is far simpler giving

$\boldsymbol{Z}|\boldsymbol{Y}=\boldsymbol{y},\lambda\sim\mathrm{N}_{m}(\boldsymbol{X}_{m}\boldsymbol{\theta}_{n},\lambda^{-1}(\boldsymbol{I}+\boldsymbol{X}_{m}\boldsymbol{A}_{n}^{-1}\boldsymbol{X}_{m}^{\prime})).$
(14)

These are already the predictive distributions if $\lambda$ is known. Let

$\displaystyle\boldsymbol{\mu}_{m}=$
$\displaystyle\boldsymbol{X}_{m}\boldsymbol{\theta}_{n}+\boldsymbol{v}_{m}^{\prime}\boldsymbol{A}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{n})$
$\displaystyle\boldsymbol{B}_{m}^{-1}=$
$\displaystyle((\boldsymbol{A}^{m})^{-1}-\boldsymbol{v}_{m}^{\prime}\boldsymbol{A}\boldsymbol{v}_{m})+(\boldsymbol{X}_{m}-\boldsymbol{v}_{m}^{\prime}\boldsymbol{A}\boldsymbol{X})\boldsymbol{A}_{n}^{-1}(\boldsymbol{X}_{m}-\boldsymbol{v}_{m}^{\prime}\boldsymbol{A}\boldsymbol{X})^{\prime}$

then, following the same calculation as in Section 4, and integrating with
respect to the posterior of $\lambda$, we obtain

$\boldsymbol{Z}|\boldsymbol{Y}=\boldsymbol{y}\sim
T_{m}^{\nu_{m}}(\boldsymbol{\mu}_{m},(\boldsymbol{B}_{m}^{*})^{-1}),$ (15)

where $\nu_{m}=2\alpha_{n}$ and
$\boldsymbol{B}_{m}^{*}=\frac{\alpha_{n}}{\beta_{n}}\boldsymbol{B}_{m}$.

###  4.3 Posterior marginal distributions

It should be clear that the marginal posterior for $\lambda$ (indeed, if
unknown) is

$\lambda\mid\boldsymbol{Y}=\boldsymbol{y}\sim Ga(\alpha_{n},\beta_{n}).$

To obtain the marginal distribution of the individual parameters $\theta_{j}$,
we prefer to take a more general approach. Let $\boldsymbol{T}$ any $k\times
p$ matrix, then

$\boldsymbol{T}\boldsymbol{\theta}\mid\boldsymbol{Y}=\boldsymbol{y},\lambda\sim
N_{k}(\boldsymbol{T}\boldsymbol{\theta}_{n},\lambda^{-1}\boldsymbol{T}\boldsymbol{A}_{n}^{-1}\boldsymbol{T}^{\prime}).$

This is already the result for $\lambda$ known. Using the posterior for
$\lambda$ we obtain, as in the previous section,

$\boldsymbol{T}\boldsymbol{\theta}\mid\boldsymbol{Y}=\boldsymbol{y}\sim
T_{k}^{\nu_{m}}\left(\boldsymbol{T}\boldsymbol{\theta}_{n},\frac{\beta_{n}}{\alpha_{n}}\boldsymbol{T}\boldsymbol{A}_{n}^{-1}\boldsymbol{T}^{\prime}\right).$

The special case where
$\boldsymbol{T}=\boldsymbol{e}_{j},\boldsymbol{T}\boldsymbol{\theta}=\theta_{j}$
is

$\theta_{j}\mid\boldsymbol{Y}=\boldsymbol{y}\sim
t_{\nu_{n}}\left(\theta^{n}_{j},\frac{\beta_{n}}{\alpha_{n}}\sigma_{j,j}^{2}\right)$
(16)

where $\boldsymbol{A}_{n}^{-1}=(\sigma_{i,j}^{2})$.

###  4.4 Efficient calculations

Assuming, in the common case, that $n>>p$, then manipulating
$\boldsymbol{\Sigma}$ is the most computationally demanding process. The
Cholesky decomposition is performed on the variance-covarince matrix
$\boldsymbol{\Sigma}=\boldsymbol{L}\boldsymbol{L}^{\prime}$. One then does
Forward Substitution to solve
$\boldsymbol{L}\boldsymbol{u}_{i}=\boldsymbol{e}_{i}$ and these are used as
row vectors to construct $\boldsymbol{U}$ to obtain
$\boldsymbol{A}=\boldsymbol{\Sigma}^{-1}=\boldsymbol{U}\boldsymbol{U}^{\prime}$
the Cholesky decomposition of the precision matrix (since
$\boldsymbol{L}\boldsymbol{u}_{i}=\boldsymbol{e}_{i}$ then
$\boldsymbol{L}\boldsymbol{U}^{\prime}=\boldsymbol{I}$ or
$\boldsymbol{U}=(\boldsymbol{L}^{-1})^{\prime}$ and note that
$\boldsymbol{A}=(\boldsymbol{L}\boldsymbol{L}^{\prime})^{-1}=(\boldsymbol{L}^{-1})^{\prime}(\boldsymbol{L}^{-1})=\boldsymbol{U}\boldsymbol{U}^{\prime}$).

Note that in calculating $\boldsymbol{A}_{n}$ the unavoidable computational
burden is in calculating
$\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{X}$. One can then solve
$\boldsymbol{A}_{n}\boldsymbol{\theta}_{n}=\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}+\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{y}$
using LU decomposition, but since this is only a $p\times p$ matrix we might
as well simply calculate the inverse. The latter is in fact a good idea in
model comparison or prediction since $\boldsymbol{A}_{n}^{-1}$ is required,
see (14), for example. We do the Cholesky and inverse decomposition of
$\boldsymbol{A}_{n}$ for all calculations. Since many calculations involve
$\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{X}$ it is a good idea first
calculate $\boldsymbol{H}=\boldsymbol{U}\boldsymbol{X}$ and use it in all
calculations.

In the case of model comparisons, see (13), the determinant of
$\boldsymbol{B}_{0}^{i}=(\boldsymbol{A}^{-1}+\boldsymbol{X}_{i}\boldsymbol{A}_{0}^{-1}\boldsymbol{X}_{i}^{\prime})^{-1}$
is required. Here use the Sherman-Woodbury-Morrison result to see that

$\boldsymbol{B}_{0}^{i}=\boldsymbol{A}-\boldsymbol{A}\boldsymbol{X}_{i}(\boldsymbol{A}_{0}+\boldsymbol{X}_{i}^{\prime}\boldsymbol{A}\boldsymbol{X}_{i})^{-1}\boldsymbol{X}_{i}^{\prime}\boldsymbol{A}=\boldsymbol{A}-\boldsymbol{A}\boldsymbol{X}_{i}(\boldsymbol{A}_{n}^{i})^{-1}\boldsymbol{X}_{i}^{\prime}\boldsymbol{A}.$

From this we may use the matrix-determinant lemma, assuming the Cholesky
decomposition of
$(\boldsymbol{A}_{0}+\boldsymbol{X}_{i}^{\prime}\boldsymbol{A}\boldsymbol{X}_{i})^{-1}=(\boldsymbol{A}_{n}^{i})^{-1}=\boldsymbol{L}_{n}^{i}(\boldsymbol{L}_{n}^{i})^{\prime}$
(this is only $p_{i}\times p_{i}$)

$\displaystyle|\boldsymbol{B}_{0}^{i}|=$
$\displaystyle\left|\boldsymbol{A}-\boldsymbol{A}\boldsymbol{X}_{i}(\boldsymbol{A}_{n}^{i})^{-1}\boldsymbol{X}_{i}^{\prime}\boldsymbol{A}\right|$
$\displaystyle=$
$\displaystyle\left|\boldsymbol{A}-(\boldsymbol{A}\boldsymbol{X}_{i}\boldsymbol{L}_{n}^{i})(\boldsymbol{A}\boldsymbol{X}_{i}\boldsymbol{L}_{n}^{i})^{\prime}\right|$
$\displaystyle=$
$\displaystyle|\boldsymbol{A}|(1-(\boldsymbol{A}\boldsymbol{X}_{i}\boldsymbol{L}_{n}^{i})^{\prime}\boldsymbol{A}^{-1}(\boldsymbol{A}\boldsymbol{X}_{i}\boldsymbol{L}_{n}^{i}))$
$\displaystyle=$
$\displaystyle|\boldsymbol{A}|(1-(\boldsymbol{L}_{n}^{i})^{\prime}(\boldsymbol{X}_{i}^{\prime}\boldsymbol{A}\boldsymbol{X}_{i})(\boldsymbol{L}_{n}^{i})).$

This should save computational burden since it is no longer required to
calculate the determinant of the $n\times n$ matrix $\boldsymbol{B}_{0}^{i}$.

Finally, in the case of predictive distributions, it is requiered to calculate

$\boldsymbol{B}_{m}^{-1}=((\boldsymbol{A}^{m})^{-1}-\boldsymbol{v}_{m}^{\prime}\boldsymbol{A}\boldsymbol{v}_{m})+(\boldsymbol{X}_{m}-\boldsymbol{v}_{m}^{\prime}\boldsymbol{A}\boldsymbol{X})\boldsymbol{A}_{n}^{-1}(\boldsymbol{X}_{m}-\boldsymbol{v}_{m}^{\prime}\boldsymbol{A}\boldsymbol{X})^{\prime}$

Using the already available Cholesky decomposition of
$\boldsymbol{A}=\boldsymbol{U}^{\prime}\boldsymbol{U}$, let
$\boldsymbol{a}_{m}=\boldsymbol{U}\boldsymbol{v}_{m}$ and then
$\boldsymbol{h}_{m}=(\boldsymbol{X}_{m}-\boldsymbol{a}_{m}^{\prime}\boldsymbol{U}\boldsymbol{X})\boldsymbol{L}_{n}$

$\boldsymbol{B}_{m}^{-1}=((\boldsymbol{A}^{m})^{-1}-\boldsymbol{a}_{m}^{\prime}\boldsymbol{a}_{m})+\boldsymbol{h}_{m}\boldsymbol{h}_{m}^{\prime}.$

##  5 Particular cases

###  5.1 Regression: formulae

The general formulas for regression are outlined in Table 1. We include an
example of a regression, simulating values from

$y_{i}=1+\sin(2\pi x_{i})+\sigma e_{i};~{}e_{i}\sim N(0,1),$

for $x_{i}=i/n;i=1,2,\ldots,n$, $n=40$ and $\sigma=0.1$. We attempt the
regressions

$\boldsymbol{Y}=\boldsymbol{X}_{p}\boldsymbol{\theta}^{p}+\epsilon$

as in (3) with uncorrelated data, i.e. $\boldsymbol{\Sigma}=\boldsymbol{I}$.
The rows of the design matrix are composed with
$\boldsymbol{X}_{p}=(\phi_{p}(x_{i}))$ with the regressor function
$\phi_{p}(x)=[1,x,x^{2},\ldots,x^{p-1}]$ and
$\boldsymbol{\theta}^{p}=[\theta_{0},\theta_{1},\ldots,\theta_{p-1}]^{\prime}$.
For $p=1,2,\ldots,6$ the regressions are performed, also calculating the
posterior probability of each model; see Figures 1, 2 and 3. For the prior,
$\boldsymbol{\theta}_{0}^{p}=\boldsymbol{0}$ and
$\boldsymbol{A}_{0}^{p}=0.001\boldsymbol{I}$, that is, the prior precision is
0.1% of the data precision (diffuse prior). $\alpha_{0}=1$ and $\beta_{0}=1$.

Table 1: Prior and posterior parameters for the Normal-Gamma conjugate family. The expressions for $\boldsymbol{\theta}_{n}$ and $\boldsymbol{A}_{n}$ are the same when $\lambda$ is a known constant, just remember that in these notes $\boldsymbol{\theta}|\lambda\sim\mathrm{N}(\boldsymbol{\theta}_{0},\lambda^{-1}\boldsymbol{A}_{0}^{-1})$ and $\boldsymbol{\theta}\sim|\boldsymbol{Y}=\boldsymbol{y},\lambda\sim\mathrm{N}(\boldsymbol{\theta}_{n},\lambda^{-1}\boldsymbol{A}_{n}^{-1})$ (normally, when $\lambda$ is known, it is not included in the variance-covariance matrix of $\boldsymbol{\theta}$). The more common, but equivalent, expression $\beta_{0}+\frac{1}{2}\left(\boldsymbol{y}^{\prime}A\boldsymbol{y}+\boldsymbol{\theta}_{0}^{\prime}\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n}^{\prime}\boldsymbol{A}_{n}\boldsymbol{\theta}_{n}\right)$ is commonly used for $\beta_{n}$. The one used here is more informative with minimal added computational burden. | Prior | Posterior | Uncorrelated Data $\boldsymbol{A}=\boldsymbol{I}$
---|---|---|---
$V(\boldsymbol{\theta}|\lambda)^{-1}$ | $\boldsymbol{A}_{0}$ | $\boldsymbol{A}_{n}=\boldsymbol{A}_{0}+\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{X}$ | $\boldsymbol{A}_{0}+\boldsymbol{X}^{\prime}\boldsymbol{X}$
$E(\boldsymbol{\theta}|\lambda)$ | $\boldsymbol{\theta}_{0}$ | $\boldsymbol{\theta}_{n}=\boldsymbol{A}_{n}^{-1}(\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}+\boldsymbol{X}^{\prime}\boldsymbol{A}\boldsymbol{y})$ | $\boldsymbol{A}_{n}^{-1}(\boldsymbol{A}_{0}\boldsymbol{\theta}_{0}+\boldsymbol{X}^{\prime}\boldsymbol{y})$
Discrepancy |  | $d^{2}_{n}=(\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n})^{\prime}\boldsymbol{A}_{0}(\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n})$ | $(\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n})^{\prime}\boldsymbol{A}_{0}(\boldsymbol{\theta}_{0}-\boldsymbol{\theta}_{n})$
Res. SS |  | $s^{2}_{n}=(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{n})^{\prime}\boldsymbol{A}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{n})$ | $(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{n})^{\prime}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{n})$
$\lambda$ shape | $\alpha_{0}$ | $\alpha_{n}=\alpha_{0}+\frac{n}{2}$ | $\alpha_{0}+\frac{n}{2}$
$\lambda$ rate | $\beta_{0}$ | $\beta_{n}=\beta_{0}+\frac{1}{2}\left(s^{2}_{n}+d^{2}_{n}\right)$ | $\beta_{0}+\frac{1}{2}\left(s^{2}_{n}+d^{2}_{n}\right)$
Figure 1:  Simulated data and the true function (black curve). From left to
right, top to bottom: MAP fit (blue), that is, the resulting regression using
$\boldsymbol{\theta}_{n}^{p}$, for $p=1,2,\ldots,6$. Figure 2:  Marginal
posterior t distributions for each parameter in the regression, for
$p=1,2,\ldots,6$, using (16). The marginal posterior for $\lambda\sim
Ga(\alpha_{n},\beta_{n})$ is not included. Figure 3:  Posterior probability of
each model, using (13).

###  5.2 Gaussian processes (GP)

GP models have many applications and represent a complex topic [5]. The basics
of the Bayesian analysis of GP’s is this. With univariate response, we have a
set of observations or measurements (or computer model evaluations)
$y(\boldsymbol{x}_{i})$ at some points or locations $i=1,2,\ldots,n$ for
$\boldsymbol{x}_{i}\in D$ in some domain. Let
$\phi:D\rightarrow\mathbb{R}^{p}$ be some regressor function as in Section
5.1; commonly $\phi(x)=[1,x_{1},x_{2},\ldots,x_{p}]$, the point coordinates.
Let the design matrix be $\boldsymbol{X}=(\phi(\boldsymbol{x}_{i}))$ and the
response vector $\boldsymbol{y}=(y(\boldsymbol{x}_{i}))$

A very important aspect of GP is its covariance structure (indeed
$\boldsymbol{A}\neq\boldsymbol{I}$). For stationary isotropic covariance the
covariance is assumed as
$cov(y(\boldsymbol{x}_{i}),y(\boldsymbol{x}_{j}))=\lambda^{-1}k(d(\boldsymbol{x}_{i},\boldsymbol{x}_{j}))$,
where $k$ is a positive defined correlation function [5]; $d$ is a metric in
$D$. From this, the covariance matrix $\boldsymbol{\Sigma}$ is formed. The
idea is to use the predictive distribution $y(\boldsymbol{x})$ at any location
$\boldsymbol{x}\in D$ as a predictor of the underlying process or function.
One uses the Normal predictive (14), if $\lambda$ is known, or otherwise the t
predictive (15). Note that in any case, translating to this context, the mean
is

$\boldsymbol{\mu}(\boldsymbol{x})=\phi(\boldsymbol{x})\boldsymbol{\theta}_{n}+\boldsymbol{v}(\boldsymbol{x})^{\prime}\boldsymbol{A}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}_{n}),$

where
$\boldsymbol{v}(\boldsymbol{x})^{\prime}=[k(d(\boldsymbol{x},\boldsymbol{x}_{i}))]$.
If $k$ is continuous at zero, the predictor is an interpolator, that is
$\mu(\boldsymbol{x}_{i})=y(\boldsymbol{x}_{i})$ with zero variance; see Figure
4.

To build a smoother one includes a discontinuity of $k$ at zero (in
geostatistics this is called “the nugget”, a measuring noise for $y$). Namely
$cov(y(\boldsymbol{x}_{i}),y(\boldsymbol{x}_{j}))=\lambda^{-1}(k(\boldsymbol{x}_{i},\boldsymbol{x}_{j})+\delta_{i,j}\sigma^{2})$.

|
---|---
Figure 4:  Gaussian process interpolator; examples of the t predictive
distribution (depicted with quantiles) ant various locations. The variance
becomes zero at observations points and
$\mu(\boldsymbol{x}_{i})=y(\boldsymbol{x}_{i})$ (interpolator).

###  5.3 Bayesian Dynamic Linear Models (DLM)

[6] is the standard reference in this case. The general DLM is explained in
section 4.2, and the univariate DLM in p. 109 of [6]. The latter is

$\displaystyle Y_{t}$
$\displaystyle=\phi_{t}\boldsymbol{\theta}_{t}+\epsilon;\epsilon\sim\mathrm{N}(0,\lambda^{-1})~{}~{}\text{(observation
equation)}$ $\displaystyle\boldsymbol{\theta}_{t}$
$\displaystyle=\boldsymbol{G}_{t}\boldsymbol{\theta}_{t-1}+\boldsymbol{\epsilon}_{1};\boldsymbol{\epsilon}_{1}\sim\mathrm{N}_{p}(\boldsymbol{0},\lambda^{-1}(\boldsymbol{A}_{t-1})^{-1})~{}~{}\text{(system
equation)},$

with some prior
$\boldsymbol{\theta}_{0}|\lambda_{0}\sim\mathrm{N}(\boldsymbol{\mu}_{0},\lambda^{-1}(\boldsymbol{A}_{0})^{-1})$.

For calculating the posterior distributions the strategy is the fallowing:
when observing the first $Y_{1}=y_{1}$, the prior (given $\lambda$) for its
parameter is
$\boldsymbol{\theta}_{1}\sim\mathrm{N}_{p}(\boldsymbol{G}_{1}\boldsymbol{\mu}_{0},\boldsymbol{G}_{1}\lambda^{-1}(\boldsymbol{A}_{0})^{-1}\boldsymbol{G}_{1}^{\prime})$,
in which case we see that the posterior is
$\boldsymbol{\theta}_{1}|Y_{1}=y_{1}\sim
N(\boldsymbol{\mu}_{1},\lambda^{-1}(\boldsymbol{A}_{1})^{-1})$, using the
corresponding formulas in Table 1. We can then repeat the process, using this
a prior to obtain the posterior
$\boldsymbol{\theta}_{2}|Y_{1}=y_{1},Y_{2}=y_{2}$ etc. These become the
“updating equations” as seen in theorem 4.1 of [6]. Moreover, one can
calculate the predictive distribution of
$Y_{t+1}|\boldsymbol{Y}_{t}=\boldsymbol{y}_{t}$ with the same formulas in
Section 4.2.

##  6 Final comments

Other simpler models may be deduced from these general calculations, as for
example standard normal sampling, setting $p=1$ and $X$ a column vector of
1’s. This last model may be used a null model in regression, and use the
Bayesian model comparison in Section 4.1.

I created a Python class, BLM, to do all the calculations and an example file
with the two examples in Regression and in GP’s, precisely the examples
presented here: https://github.com/andreschristen/BLM

## References

  * [1] Lyle D Broemeling.  Bayesian analysis of linear models.  CRC Press, 2017.
  * [2] Morris L. Eaton.  Multivariate Statistics.  SPIE, January 2007.
  * [3] Jennifer A. Hoeting, David Madigan, Adrian E. Raftery, and Chris T. Volinsky.  Bayesian model averaging: A tutorial.  Statistical Science, 14(4):382–401, 1999.
  * [4] A. O’Hagan and J. Forster.  Bayesian Inference, volume 2B.  Hodder Arnold, London, 1994.
  * [5] Carl Edward Rasmussen and Christopher K. I. Williams.  Gaussian processes for machine learning. Adaptive computation and machine learning. MIT Press, 2006.
  * [6] Mike West and Jeff Harrison.  Bayesian Forecasting and Dynamic Models (Springer Series in Statistics).  Springer-Verlag, 1997.


††thanks: Corresponding author

# Heavy quark potential and jet quenching parameter in a rotating D-instanton
background

Jun-Xia Chen<EMAIL_ADDRESS>Institute of Particle Physics and
Key Laboratory of Quark and Lepton Physics (MOS), Central China Normal
University, Wuhan 430079, China  De-Fu Hou<EMAIL_ADDRESS>Institute
of Particle Physics and Key Laboratory of Quark and Lepton Physics (MOS),
Central China Normal University, Wuhan 430079, China

###### Abstract

We get the dual gravity metric of the rotating nuclear matter by performing a
standard Lorentz transformation on the static metric in the D-instanton
background. Then, we study the effects of the angular velocity, the instanton
density and the temperature on the heavy quark potential. It is shown that the
angular velocity and the temperature promote dissociation of the quark pair,
and the instanton density suppress dissociation. Similarly, according to the
result of jet quenching parameter, we found that the jet quenching parameter
increases with the increase of angular velocity, instanton density and
temperature, and the jet quenching parameter in the rotating D-instanton
background is larger than that of $\mathcal{N}=4$ SYM theory.

##  I Introduction

Many efforts have been made to investigate the properties of the quark-gluon
plasma (QGP) Shuryak (2005); et al. (STAR Collaboration, PHENIX Collaboration)
since people discovered this new form of matter created in heavy ion
collisions. Usually, the free quark and gluon can not be observed because of
color confinement. The QGP produced in the laboratory exists for a short time
so that only leptons and hadrons are measured by the detectors of the
collider. According to the information of the final state particles, we can
study the properties of quark-gluon plasma by using QGP probes, such as heavy
quarkonium Maldacena (1998), jet Casalderrey-Solana  _et al._ (2014a); Wang
and Gyulassy (1992) and so on.

Last century, Maldacena conjectured that the strongly coupled gauge theory in
four-dimensional space-time is dual to the weakly coupled string theory on
$AdS_{5}\times S^{5}$, referred to as AdS/CFT duality Aharony  _et al._
(2000); Maldacena (1999). According to AdS/CFT, many observable quantities can
be calculated. In this paper, we will study the heavy quark potential and jet
quenching parameter using AdS/CFT.

The heavy quark potential $V(L)$ increases with the increase of the distance
$L$ between the quark pair. However, when the separation distance reaches the
maximum, the heavy quark potential will no longer increase and is a constant
value, which implies that quark-antiquark interaction is screened by the QGP
between them. According to the function image between the heavy quark
potential and the separation distance, we can know the maximum value of the
heavy quark potential and the maximum separation distance of the quark pair.
The heavy quark potential for $\mathcal{N}=4$ SYM at finite temperature was
calculated in Rey  _et al._ (1998); Brandhuber  _et al._ (1998). The
subleading term of this observable has been investigated in Zhang  _et al._
(2011). Many people have studied the influences of different parameters on the
heavy quark potential. For example, the effect of hyperscaling violation on
the potential was considered in Zhang  _et al._ (2017a); Kioumarsipour and
Sadeghi (2019). The influence of magnetic fields on the potential has been
investigated in Rougemont  _et al._ (2015). The heavy quark potential as a
function of shear viscosity at strong coupling has been discussed in Noronha
and Dumitru (2009). In Jahnke and Misobuchi (2016), the author studied the
potential for different values of the Gauss-Bonnet coupling $\lambda_{GB}$.
Other important results can be seen in Andreev and Zakharov (2006); Chernicoff
_et al._ (2006); Albacete  _et al._ (2008); Boschi-Filho  _et al._ (2006);
Avramis  _et al._ (2007); Yang and Yuan (2015); Nata Atmaja  _et al._ (2015);
Wu  _et al._ (2019).

The jet quenching parameter is defined as the mean transverse momentum
acquired by the hard parton per unit distance travelled Casalderrey-Solana
_et al._ (2014b). According to this parameter, we kown how rapidly the
energetic parton loses its energy, though $\hat{q}$ is defined via transverse
momentum broadening only. The frist calculation of the jet quenching parameter
in hot $\mathcal{N}=4$ supersymmetric QCD was done by H.Liu et al Liu  _et
al._ (2006). They found $\hat{q}_{SYM}=26.69\sqrt{\alpha_{SYM}N_{c}}T^{3}$ in
the large $N_{c}$ and large $\lambda$ limite. In Armesto  _et al._ (2006), the
first correction to the parameter has been computed. After that, more and more
people study the effects of different quantities on the jet quenching
parameter. For instance, the authors investigated the influence of the
magnetic fields on the paramerter in Rougemont (2020); Zhu  _et al._ (2019);
Li  _et al._ (2016). A lattice study of the jet quenching parameter can be
seen in Panero  _et al._ (2014). The parameter in strongly coupled anisotropic
$\mathcal{N}=4$ plasma has been studied in Giataganas (2012); Chernicoff  _et
al._ (2012). Other related works can be seen in D’Eramo  _et al._ (2011);
Bitaghsir Fadafan  _et al._ (2011); Li  _et al._ (2014); Zhang and Zhu (2019);
Liu  _et al._ (2007); Buchel (2006); Nakano  _et al._ (2007); Zhu  _et al._
(2021).

An instanton is a solution to equations of motion of the classical field
theory on a Euclidean space-time. In order to study the instanton effect from
holography, we will work in a D-instanton background, since the D-instanton in
string theory is dual to the instanton in super Yang-Mills theory Park and Sin
(1998); Li and Lin (2018). The D-instanton background is first proposed in Liu
and Tseytlin (1999). They consider D-instantons homogeneously distributed over
D3-brane, which corresponds to modifying the pure D3-brane background by
adding a RR scalar charge and the dilaton charge. The type IIB string theory
in near-horizon D3/D-instanton background is dual to $\mathcal{N}=4$ SYM
theory in a constant self-dual gauge field background. The effect of
D-instanton density on different quantities has attracted much interest, such
as drag force and diffusion coefficient Zhang  _et al._ (2018), imaginary
potential Zhang  _et al._ (2017b), Schwinger effect Shahkarami  _et al._
(2018), and so on.

Recently, it was found that the fluid produced in non-central collisions have
a strong vortical structure Adamczyk  _et al._ (2017). Particles become
globally polarized along the direction of rotation because of spin-orbit
coupling. According to the results of $\varLambda$ hyperon polarization
measurements, The STAR Collaboration found there was vortex in Au+Au
collisions at low energy ($\sqrt{s_{NN}}$ = 7.7-39GeV) Adamczyk  _et al._
(2017) while no vortex at $\sqrt{s_{NN}}$ = 62.4GeV and 200GeV Abelev  _et
al._ (2007). The QGP is the most vortical system so far observed, so the
effect of vortices on the properties of QGP has attracted intense interest.
For example, the influence of rotation on the phase diagram has been studied
in Fujimoto  _et al._ (2021); Chernodub (2021); Chernodub and Gongyo (2017).
The effect of rotation on the thermodynamic quantities and equations of state
(EoS) has been investigated in Zhou  _et al._ (2021). The jet quenching
parameter in a rotating nuclear matter was calculated in Golubtsova  _et al._
(2021) using 5d Kerr-AdS black hole. Other interesting results can be seen in
McInnes (2018, 2020); Adam  _et al._ (2018); Ebihara  _et al._ (2017);
Becattini and Lisa (2020); McInnes (2016).

In order to study the effect of rotation, we will study the heavy quark
potential and jet quenching parameter in a rotating D-instanton background.
Actually, in Zhang  _et al._ (2016), the authors have calculated these two
quantities in a D-instanton background and investigated the effect of
instanton density on them. However, according to Adamczyk  _et al._ (2017),
the QGP generated in the heavy ion collision is rotational, so we want to
study the influence of the rotation on these two quantities in this paper. To
get the metric in a rotating D-instanton background, we should solve the
equation of motion of the rotating nuclear matter. However, it is very
difficult. We will use the method in Bravo Gaete  _et al._ (2017); Erices and
Martinez (2018); Chen  _et al._ (2020) to get a non-conformal rotating black
hole solution by performing a standard Lorentz transformation on the static
black hole solution. Due to the symmetry of the action, it must be a solution
of the Einstein equation. Then, we will calculate the heavy quark potential
and jet quenching parameter using the resulting metric and study the effects
of angular velocity, instanton density and temperature on these two
quantities.

The organization of the paper is as follows. In Sec.II, we get the metric in a
rotating D-instanton background by performing a standard Lorentz
transformation on the static metric in the D-instanton background. In Sec.III,
we study the effects of angular velocity, instanton density and temperature on
the heavy quark potential. In Sec.IV, similarly, we study the effects of
angular velocity, instanton density and temperature on the jet quenching
parameter. In Sec.V, we give some conclusion and discussion.

##  II Background geometry

Let’s briefly review the D-instanton background. We consider D-instantons
homogeneously distributed over D3-brane, which corresponds to modifying the
pure D3-brane background by adding a RR scalar charge and the dilaton charge.
The ten dimensional super-gravity action in Einstein frame is Gibbons  _et
al._ (1996); Kehagias and Sfetsos (1999)

$S=\frac{1}{\kappa}\int
d^{10}x\sqrt{g}(R-\frac{1}{2}(\partial\Phi)^{2}+\frac{1}{2}e^{2\Phi}(\partial\chi)^{2}-\frac{1}{6}F^{2}_{(5)})$
(1)

where $\Phi$ and $\chi$ denote the dilaton and the axion respectively, and
$F_{(5)}$ is a five-form field strength coupled to the D3-branes. If we set
$\chi=-e^{-\Phi}+\chi_{0}$, the dilaton term cancel the axion term in Eq.(1).
Then the solution in the string frame is Gwak  _et al._ (2012)

$ds^{2}=e^{\frac{\Phi}{2}}[-\frac{r^{2}}{R^{2}}f(r)dt^{2}+\frac{r^{2}}{R^{2}}(dx_{1}^{2}+dx_{2}^{2})+\frac{r^{2}}{R^{2}}d\varphi^{2}+\frac{1}{f(r)}\frac{R^{2}}{r^{2}}dr^{2}]+R^{3}d\Omega_{5}^{2}$
(2)

with

$e^{\Phi}=1+\frac{q}{r_{t}^{4}}\log\frac{1}{f(r)},\quad
f(r)=1-\frac{r_{t}^{4}}{r^{4}}$ (3)

where $(t,x_{1},x_{2},\varphi,r)$ are coordinates of $AdS_{5}$, $r_{t}$ is the
radial position of the event horizon, and gauge theory is on the boundary
$r\rightarrow\infty$. $d\Omega_{5}$ denotes the volume element of the five-
sphere $S^{5}$, $R$ is the radius of curvature, the parameter $q$ denotes the
D-instanton density, in the dual picture, it also represents the vacuum
expectation value of gluon condensation Shahkarami and Charmchi (2019).

The large vortex will be produced if the heavy-ion collision is a non-central
collision. In order to investigate the effect of the rotation on the QGP, we
should calculate the observables in the rotating background. It is very
difficult to solve the equation of motion of the rotating nuclear matter from
the action. Here, we will take an approximate approach. According to Chen  _et
al._ (2020); Bravo Gaete  _et al._ (2017); Erices and Martinez (2018), the
metric in rotating background can be obtained by performing a standard Lorentz
transformation on the static metric Eq.(2)

$t\longrightarrow\frac{1}{\sqrt{1-(\omega l)^{2}}}(t+\omega
l^{2}\varphi),\quad\varphi\longrightarrow\frac{1}{\sqrt{1-(\omega
l)^{2}}}(\varphi+\omega t)$ (4)

where $\omega$ is the angular velocity, $l$ is the radius to the rotating
axis. Since we only focus on the qualitative results, we will take
$l=1GeV^{-1}$. The resulting metric is

$\begin{split}ds^{2}=&e^{\frac{\Phi}{2}}\frac{r^{2}}{R^{2}}[\frac{1}{1-\omega^{2}}(\omega^{2}-f(r))dt^{2}+\frac{1}{1-\omega^{2}}(1-\omega^{2}f(r))d\varphi^{2}+\frac{2\omega}{1-\omega^{2}}(1-f(r))dtd\varphi\\\
&+\frac{1}{f(r)}\frac{R^{4}}{r^{4}}dr^{2}+dx_{1}^{2}+dx_{2}^{2}]\end{split}$
(5)

In order to calculate the temperature of the black hole, we need to transform
the Eq.(5) into Eq.(6) below

$ds^{2}=H(r)(d\varphi+F(r)dt)^{2}+N(r)dt^{2}+e^{\frac{\Phi}{2}}\frac{1}{f(r)}\frac{R^{2}}{r^{2}}dr^{2}+e^{\frac{\Phi}{2}}\frac{r^{2}}{R^{2}}(dx_{1}^{2}+dx_{2}^{2})$
(6)

where

$\begin{split}&H(r)=e^{\frac{\Phi}{2}}\frac{r^{2}}{R^{2}}\frac{1}{1-\omega^{2}}(1-\omega^{2}f(r)),\quad
F(r)=\omega(1-f(r))(1-\omega^{2}f(r))^{-1}\\\
&N(r)=e^{\frac{\Phi}{2}}\frac{r^{2}}{R^{2}}(\omega^{2}-1)f(r)(1-\omega^{2}f(r))^{-1}\end{split}$
(7)

The Hawking temperature of the black hole is given by

$T=|\frac{\lim\limits_{r\to
r_{t}}-\frac{1}{2}\sqrt{\frac{g^{rr}}{-\hat{g}_{00}}}\hat{g}_{00,1}}{2\pi}|=\frac{r_{t}\sqrt{1-\omega^{2}}}{\pi
R^{2}}$ (8)

here $\hat{g}_{00}$ is different from the $t-t$ component of the metric
$g_{00}$. $g^{rr}$ is also different from $g_{rr}$. $\hat{g}_{00,1}$ is the
derivative of $\hat{g}_{00}$

$\hat{g}_{00}=N(r),\quad g^{rr}=\frac{1}{g_{rr}}$ (9)

##  III Heavy quark potential

In this section, we study the effects of angular velocity, D-instanton
density, and temperature on the heavy quark potential. The quark-antiquark
pair is represented by a string whose endpoints are on the D-brane, and the
D-brane is at the boundary of AdS. The string which is not straight but goes
inside the AdS space-time is the lowest energy state. For simplicity, such a
string can be approximated by a rectangular string which can be divided into
two parts:the vertical part and the horizontal part. Only the horizontal part
contributes to the quark potential, and the vertical part gives the quark
mass.

The heavy quark potential can be extracted from the expectation value of the
Wilson loop. In the limit $\mathcal{T}\rightarrow\infty$, we have

$\langle W(C)\rangle\sim e^{-\mathcal{T}V(L)}$ (10)

where $C$ is a rectangular loop, one side represents the time $\mathcal{T}$,
the other side represents the separation distance $L$ of the quark pair in the
AdS space-time, and $V(L)$ represents the heavy quark potential.

In addition, in the large $N_{c}$ and large $\lambda$ limit, we have

$\langle W(C)\rangle\sim e^{-S_{c}}$ (11)

where $S_{c}=S-S_{0}$, $S$ represents the Nambu-Goto action of the U-shaped
string, $S_{0}$ represents the Nambu-Goto action of the straight string. The
Nambu-Goto action is

$S=-\frac{1}{2\pi\alpha^{\prime}}\int d\tau d\sigma\sqrt{-detg_{\alpha\beta}}$
(12)

where $\frac{1}{2\pi\alpha^{\prime}}$ is string tension, $g_{\alpha\beta}$ is
known as induced metric

$g_{\alpha\beta}=G_{\mu\nu}\frac{\partial
X^{\mu}}{\partial\sigma^{\alpha}}\frac{\partial
X^{\nu}}{\partial\sigma^{\beta}}$ (13)

Combining the formula (10) and (11), we can represent the heavy quark
potential in terms of the Nambu-Goto action

$V(L)=\frac{S_{c}}{\mathcal{T}}$ (14)

Now, we are going to calculate the action of the string. A particle draws a
world-line in spacetime. A string sweeps a world-sheet in spcetime. We can
take the world-sheet coordinates as $\sigma^{\alpha}=(\tau,\sigma)$. Then, the
string motion is described by $X^{\mu}(\sigma^{\alpha})$. Here, we take the
static gauge

$\tau=t,\quad\sigma=x_{1},\quad r=r(\sigma)$ (15)

Using the world-sheet coordinates $\sigma^{\alpha}$, the $AdS_{5}$ spacetime
metric is

$\begin{split}ds^{2}&=G_{\mu\nu}dX^{\mu}dX^{\nu}=g_{\alpha\beta}d\sigma^{\alpha}d\sigma^{\beta}\\\
&=e^{\frac{\Phi}{2}}\frac{r^{2}}{R^{2}}[\frac{1}{1-\omega^{2}}(\omega^{2}-f(r))d\tau^{2}+(1+\frac{1}{f(r)}\frac{R^{4}}{r^{4}}\dot{r}^{2})d\sigma^{2}]\end{split}$
(16)

Then, the action is given by

$S=\frac{\mathcal{T}}{2\pi\alpha^{\prime}}\int
d\sigma\sqrt{e^{\Phi}\frac{1}{1-\omega^{2}}(f(r)-\omega^{2})(\frac{\dot{r}^{2}}{f(r)}+\frac{r^{4}}{R^{4}})}$
(17)

where $\dot{r}=\frac{dr}{d\sigma}$, $\mathcal{T}$ is the time duration in
$t$.According to the action, the Lagrangian does not contain $\sigma$, so we
have a conserved quantity

$\frac{\partial\mathcal{L}}{\partial\dot{r}}\dot{r}-\mathcal{L}=constant$ (18)

Let us determine the constant. Assuming that the U-shaped string has the
turning point at $\sigma=0,r=r_{c}$, then, $\dot{r}\mid_{r=r_{c}}=0$, so the
left of Eq.(18) becomes a constant when we set $r=r_{c}$

$\frac{e^{\Phi}\frac{1}{1-\omega^{2}}(f(r)-\omega^{2})\frac{r^{4}}{R^{4}}}{\sqrt{e^{\Phi}\frac{1}{1-\omega^{2}}(f(r)-\omega^{2})(\frac{\dot{r}^{2}}{f(r)}+\frac{r^{4}}{R^{4}})}}=\sqrt{e^{\Phi(r_{c})}\frac{1}{1-\omega^{2}}(f(r_{c})-\omega^{2})\frac{r_{c}^{4}}{R^{4}}}$
(19)

Figure 1: The left is U-shaped string, the right is straight string.

After a simple calculation, we have

$\dot{r}=\frac{r^{2}}{R^{2}}\sqrt{f(r)\frac{g(r)-g(r_{c})}{g(r_{c})}}$ (20)

where

$\begin{split}&g(r)=e^{\Phi}(f(r)-\omega^{2})r^{4},\quad
g(r_{c})=e^{\Phi(r_{c})}(f(r_{c})-\omega^{2})r_{c}^{4}\\\
&e^{\Phi(r_{c})}=1+\frac{q}{r_{t}^{4}}\log\frac{1}{f(r_{c})},\quad
f(r_{c})=1-\frac{r_{t}^{4}}{r_{c}^{4}}\end{split}$ (21)

As shown in the Fig.1, the quark and antiquark are located at
$x_{1}=\frac{L}{2},x_{1}=-\frac{L}{2}$ respectively if the distance of the
quark pair is $L$. Integrating the above equation, we get the distance between
quark-antiquark pair

$L=2R^{2}\int_{r_{c}}^{\infty}dr\sqrt{\frac{g(r_{c})}{r^{4}f(r)[g(r)-g(r_{c})]}}$
(22)

Figure 2: Heavy quark potential $V(L)$ as a function of distance between quark
pair. (a)for different angular velocities $\omega=0.1$(red line),
$\omega=0.3$(blue line), $\omega=0.5$(black line), here, we set the
temperature $T=0.2$, the instanton density $q=0.5$. (b)for different instanton
density $q=0.2$(red line), $q=0.5$(blue line), $q=1$(black line), here, we set
temperature $T=0.2$, the angular velocities $\omega=0.1$. (c)for different
temperature $T=0.2$(red line), $T=0.25$(blue line), $T=0.3$(black line), here,
we set the angular velocities $\omega=0.1$, the instanton density $q=0.5$.The
unit for $T,\omega$ is in GeV.

Combining Eq.(17) and Eq.(20), the action of U-shaped string is

$S=\frac{\mathcal{T}}{\pi\alpha^{\prime}}\int_{r_{c}}^{\infty}dr\sqrt{\frac{e^{\Phi}(f(r)-\omega^{2})g(r)}{(1-\omega^{2})f(r)[g(r)-g(r_{c})]}}$
(23)

it represents the total energy of the quark-antiquark pair. We need to
subtract the masses of the quark and antiquark to get the potential between
the quark pair. The action of the straight string connecting the boundary and
the horizon represents the mass of the quark

$S_{0}=\frac{\mathcal{T}}{\pi\alpha^{\prime}}\int_{r_{t}}^{\infty}dr\sqrt{-g_{tt}g_{rr}}=\frac{\mathcal{T}}{\pi\alpha^{\prime}}\int_{r_{t}}^{\infty}dr\sqrt{\frac{e^{\Phi}(1-\omega^{2})}{1-\omega^{2}f(r)}}$
(24)

The heavy quark potential is

$V(L)=\frac{1}{\pi\alpha^{\prime}}\\{\int_{r_{c}}^{\infty}dr\sqrt{\frac{e^{\Phi}(f(r)-\omega^{2})g(r)}{(1-\omega^{2})f(r)[g(r)-g(r_{c})]}}-\int_{r_{t}}^{\infty}dr\sqrt{\frac{e^{\Phi}(1-\omega^{2})}{1-\omega^{2}f(r)}}\\}$
(25)

We can set $\frac{1}{\pi\alpha^{\prime}}=1,R=1$, because we just qualitatively
analyze the influences of angular velocity $\omega$, instanton density $q$ and
temperature on the heavy quark potential. As shown in Fig.2, we plot the curve
of heavy quark potential in terms of the distance between quark pair. There
are two parts for each curve:the upper part and the lower part. In fact, only
the lower part is significant. The corresponding $L=L_{m}$ at the turning
point of the curve is the maximum separation distance which is also called
screening length. At this point, the interaction between quark-antiquark has
been screened by QGP. When the distance between quark pair increases further,
the heavy quark potential will not increase, that is, when $L\geq L_{m}$, the
heavy quark potential is a constant. In the dual geometry, we believe that the
bottom of the U-shaped string has reached the horizon at this point.

From Fig.2 a), we find that the screening length decreases as the angular
velocity increases when we set the temperature $T=0.2GeV$, the instanton
density $q=0.5$, and the heavy quark potential increases with the increase of
angular velocity at the same separation distance. Both of them indicate that
the angular velocity $\omega$ promotes dissociation of the quark pair. From
Fig.2 b), we find that the screening length increases as the instanton density
increases when we set the temperature $T=0.2GeV$, the angular velocity
$\omega=0.1GeV$, and the heavy quark potential decreases with the increase of
instanton density at the same separation distance. Both of them indicate that
the instanton density $q$ suppress dissociation. From Fig.2 c), we find that
the screening length decreases as the temperature increases when we set the
angular velocity $\omega=0.1GeV$, the instanton density $q=0.2$, and the heavy
quark potential increases with the increase of temperature at the same
separation distance. Both of them indicate that high temperature promotes
dissociation.

##  IV Jet quenching parameter

In this section, we will calculate the jet quenching parameter in a rotating
D-instanton background according to AdS/CFT, and study the influences of
angular velocity, instanton density, and temperature on jet quenching
parameter. It can be extracted from the expectation value of an adjoint Wilson
loop

$\langle W^{A}(\mathcal{C})\rangle\approx
exp[-\frac{1}{4\sqrt{2}}\hat{q}L_{-}L^{2}]$ (26)

here, we assume the distance bettwen quark pair is $L$, and they propagate
along the light-cone through the thermal plasma for a distance $L^{-}$.
Contour $\mathcal{C}$ is a rectangular loop, one side is along transverse
direction with length $L$, the other side is in the $x^{-}$ direction with
length $L^{-}$.

The expectation value of lightlike Wilson loop in the fundamental
representation refers to the regularized action. According to AdS/CFT
correspondence, in the large $N_{c}$ and large $\lambda$ limite, the
expectation value is given by

$\langle W^{F}(\mathcal{C})\rangle\approx exp[-S_{c}]$ (27)

Combining Eq.(26) and Eq.(27) and using the relation $\langle
W^{A}(\mathcal{C})\rangle\sim\langle W^{F}(\mathcal{C})\rangle^{2}$, we get
the expression of the jet quenching parameter

$\hat{q}=8\sqrt{2}\frac{S_{c}}{L_{-}L^{2}}$ (28)

Using light-cone coordinates $x^{\mu}=(r,x^{+},x^{-},x_{2},\varphi)$, the
metric Eq.(5) becomes

$\begin{split}ds^{2}=&e^{\frac{\phi}{2}}\frac{r^{2}}{R^{2}}\\{\frac{1}{2}\frac{1-f(r)}{1-\omega^{2}}[(dx^{+})^{2}+(dx^{-})^{2}]+[\frac{\omega^{2}-f(r)}{1-\omega^{2}}-1]dx^{+}dx^{-}+\frac{1-\omega^{2}f(r)}{1-\omega^{2}}d\varphi^{2}\\\
&+\frac{\sqrt{2}\omega(1-f(r))}{1-\omega^{2}}(dx^{+}+dx^{-})d\varphi+dx_{2}^{2}\\}+e^{\frac{\phi}{2}}\frac{1}{f(r)}\frac{R^{2}}{r^{2}}dr^{2}\end{split}$
(29)

We still take the world-sheet coordinates as $\sigma^{\alpha}=(\tau,\sigma)$.
The Nambu-Goto action is invariant under coordinate transformation on
$\sigma^{\alpha}$, so we can set $\tau=x^{-},\sigma=x_{2}$. Here, we assume
the quark pair is along $x_{2}$. The Wilson loop is on the surface where
$\varphi$ and $x^{+}$ are constant. After gauge fixing, the metric Eq.(29)
becomes

$ds^{2}=\frac{1}{2}e^{\frac{\Phi}{2}}\frac{r^{2}}{R^{2}}[\frac{1}{1-\omega^{2}}(\omega^{2}-f(r))+1]d\tau^{2}+e^{\frac{\Phi}{2}}(\frac{1}{f(r)}\frac{R^{2}}{r^{2}}\dot{r}^{2}+\frac{r^{2}}{R^{2}})d\sigma^{2}$
(30)

where $\dot{r}=\frac{dr}{d\sigma}$, using Eq.(12) and Eq.(13), the action of
the string is

$S=\frac{\sqrt{2}L_{-}}{2\pi\alpha^{\prime}}\int_{0}^{\frac{L}{2}}d\sigma\sqrt{e^{\Phi}\frac{1}{1-\omega^{2}}(1-f(r))(\frac{\dot{r}^{2}}{f(r)}+\frac{r^{4}}{R^{4}})}$
(31)

Similarly, the Lagrangian does not contain $\sigma$, we have

$\frac{\partial\mathcal{L}}{\partial\dot{r}}\dot{r}-\mathcal{L}=C$ (32)

Solving the above equation, we get

$\dot{r}=\frac{r^{2}}{R^{2}}\sqrt{f(r)[\frac{1}{C^{2}}e^{\Phi}\frac{1}{1-\omega^{2}}(1-f(r))\frac{r^{4}}{R^{4}}-1]}$
(33)

according to the expression of $\dot{r}$, when $f(r)=0$, it is easy to get
$\dot{r}=0$. In the other words, the turning point of the U-shaped string is
at $r=r_{t}$.In the low energy limit ($C\rightarrow 0$), the action of the
U-shaped string is

$S=\frac{\sqrt{2}L_{-}}{2\pi\alpha^{\prime}}\int_{r_{t}}^{\infty}dr\sqrt{e^{\Phi}\frac{1-f(r)}{(1-\omega^{2})f(r)}}(1+\frac{1}{2\frac{1}{C^{2}}e^{\Phi}\frac{1}{1-\omega^{2}}(1-f(r))\frac{r^{4}}{R^{4}}})$
(34)

it represents the total energy of the quark pair. We need to subtract the
energy of the free quark which can be calculated by the action of the straight
string connecting the boundary and the horizon

$S_{0}=\frac{2L_{-}}{2\pi\alpha^{\prime}}\int_{r_{t}}^{\infty}dr\sqrt{g_{--}g_{rr}}=\frac{\sqrt{2}L_{-}}{2\pi\alpha^{\prime}}\int_{r_{t}}^{\infty}\sqrt{e^{\Phi}\frac{1-f(r)}{(1-\omega^{2})f(r)}}$
(35)

Figure 3: (a) $\hat{q}/\hat{q}_{SYM}$ as a function of angular velocity
$\omega$, (b) jet quenching parameter $\hat{q}$ as a function of $\omega$. The
red line represents $q=0.2,T=0.3GeV$, the blue line represents $q=1,T=0.2GeV$.
Figure 4: (a) $\hat{q}/\hat{q}_{SYM}$ as a function of the instanton density
$q$, (b) jet quenching parameter $\hat{q}$ as a function of the instanton
density $q$. The red line represents $\omega=0.2GeV,T=0.3GeV$, the blue line
represents $\omega=0.3GeV,T=0.2GeV$. Figure 5: (a) $\hat{q}/\hat{q}_{SYM}$ as
a function of the temperature $T$, (b) the jet quenching parameter $\hat{q}$
as a function of the temperature $T$. The red line represents
$\omega=0.3GeV,q=0.5$, the blue line represents $\omega=0.5GeV,q=0.3$.

Combining the Eq.(28), Eq.(34) and Eq.(35), the jet quenching parameter is

$\hat{q}=\frac{1}{\pi\alpha^{\prime}}\frac{1}{\int_{r_{t}}^{\infty}\sqrt{\frac{1-\omega^{2}}{e^{\Phi}f(r)(1-f(r))}}\frac{R^{4}}{r^{4}}dr}$
(36)

We plot the curve of $\hat{q}/\hat{q}_{SYM}$ and $\hat{q}$ in terms of the
angular velocity $\omega$, instanton density $q$ and temperature $T$ shown
separately in Fig.3, Fig.4 and Fig.5. Here we set
$\frac{1}{\pi\alpha^{\prime}}=1,R=1$, because we just qualitatively analyze
the influences of angular velocity, instanton density and temperature on the
jet quenching parameter. From Fig.3 b), we find the jet quenching parameter
increases with the increase of the angular velocity, which indicates that
rotation promotes the energy loss. From Fig.4 b), we find the jet quenching
parameter increases with the increase of the instanton density, which
indicates that the larger the instanton density is, the faster the energy
loses. From Fig.5 b), we find the jet quenching parameter increases with the
increase of the temperature, which indicates that high temperature promotes
the energy loss. From Figure a) of Fig.3, Fig.4 and Fig.5, it can be found
that the jet quenching parameter in a rotating D-instanton background is
larger than that of $\mathcal{N}=4$ SYM theory. It is consistent with the
conclusion we obtained according to Figure b) of Fig.3, Fig.4 and Fig.5.

In order to compare the results with the experimental data, we set
$\alpha^{\prime}=0.5$, which is reasonable for temperatures not far above the
QCD phase transition Liu  _et al._ (2006). According to Gubser (2007), the
value of ’t Hooft coupling constant is $5.5<\lambda<6\pi$. Here, we set
$\lambda=6\pi$, the temperature $T=250MeV$, the instanton density $q=1$. From
Eq.(36), we find $\hat{q}=5.4,6.3,9.2GeV^{2}/fm$ for $\omega=0.1,0.3,0.5GeV$,
which is consistent with the RHIC data($5\rightarrow 25GeV^{2}/fm$) Edelstein
and Salgado (2008).

##  V Conclusion and discussion

In this paper, we study the effects of vorticity, instanton density and
temperature on the properties of QGP. Here, two types of QGP probes, heavy
quarkion and jet, are used.In order to study the influence of the vorticity on
the QGP, a coordinate transformation is performed on the static black hole
solution to obtain the non-conformal rotating black hole solution.

According to AdS/CFT, we calculate the heavy quark potential in the rotating
D-instanton background and plot the curve of heavy quark potential in terms of
the distance between quark pairs. We find that the angular velocity and the
temperature promote dissociation of the quark pair, and the instanton density
suppress dissociation. Similarly, we also study the effects of angular
velocity, instanton density, and temperature on jet quenching parameter, and
find that the jet quenching parameter increases with the increase of angular
velocity, instanton density and temperature, and the jet quenching parameters
in the rotating D-instanton background are larger than that of $\mathcal{N}=4$
SYM theory. Furthermore, the values of the jet quenching parameter we obtained
are consistent with the RHIC data when we set $\alpha^{\prime}=0.5$,
$\lambda=6\pi$.

We expect the effect of vorticity on the phase transition and heavy quark
energy loss we find in this paper can help us to better understand the results
in the experiment. The effect of vorticity on the QGP is as important as the
influence of magnetic field, so it is significant to study the effect of
vorticity on the other quantities, such as imaginary potential and drag force
and so on.

## Acknowledgments

This work is supported by the National Natural Science Foundation of China by
Grant Nos. 11735007, 11890711, and 11890710. We would like to thank Hai-cang
Ren , Ziqiang Zhang and Zhou-Run Zhu for useful discussions.

## References

## References

  * Shuryak (2005) E. Shuryak, Nuclear Physics A 750, 64 (2005).
  * et al. (STAR Collaboration) J. A. et al. (STAR Collaboration), Nuclear Physics A 757, 102 (2005).
  * et al. (PHENIX Collaboration) K. A. et al. (PHENIX Collaboration), Nuclear Physics A 757, 184 (2005).
  * Maldacena (1998) J. M. Maldacena, Phys. Rev. Lett. 80, 4859 (1998), arXiv:hep-th/9803002 .
  * Casalderrey-Solana  _et al._ (2014a) J. Casalderrey-Solana, D. C. Gulhan, J. G. Milhano, D. Pablos, and K. Rajagopal, Nucl. Phys. A 931, 487 (2014a), arXiv:1408.5616 [hep-ph] .
  * Wang and Gyulassy (1992) X.-N. Wang and M. Gyulassy, Phys. Rev. Lett. 68, 1480 (1992).
  * Aharony  _et al._ (2000) O. Aharony, S. S. Gubser, J. M. Maldacena, H. Ooguri, and Y. Oz, Phys. Rept. 323, 183 (2000), arXiv:hep-th/9905111 .
  * Maldacena (1999) J. M. Maldacena, AIP Conf. Proc. 484, 51 (1999).
  * Rey  _et al._ (1998) S.-J. Rey, S. Theisen, and J.-T. Yee, Nucl. Phys. B 527, 171 (1998), arXiv:hep-th/9803135 .
  * Brandhuber  _et al._ (1998) A. Brandhuber, N. Itzhaki, J. Sonnenschein, and S. Yankielowicz, Phys. Lett. B 434, 36 (1998), arXiv:hep-th/9803137 .
  * Zhang  _et al._ (2011) Z.-q. Zhang, D. Hou, H.-c. Ren, and L. Yin, JHEP 07, 035 (2011), arXiv:1104.1344 [hep-ph] .
  * Zhang  _et al._ (2017a) Z.-q. Zhang, C. Ma, D.-f. Hou, and G. Chen, Adv. High Energy Phys. 2017, 8276534 (2017a), arXiv:1604.04349 [hep-ph] .
  * Kioumarsipour and Sadeghi (2019) M. Kioumarsipour and J. Sadeghi, Eur. Phys. J. C 79, 636 (2019).
  * Rougemont  _et al._ (2015) R. Rougemont, R. Critelli, and J. Noronha, Phys. Rev. D 91, 066001 (2015), arXiv:1409.0556 [hep-th] .
  * Noronha and Dumitru (2009) J. Noronha and A. Dumitru, Phys. Rev. D 80, 014007 (2009), arXiv:0903.2804 [hep-ph] .
  * Jahnke and Misobuchi (2016) V. Jahnke and A. S. Misobuchi, Eur. Phys. J. C 76, 309 (2016), arXiv:1510.03774 [hep-th] .
  * Andreev and Zakharov (2006) O. Andreev and V. I. Zakharov, Phys. Rev. D 74, 025023 (2006), arXiv:hep-ph/0604204 .
  * Chernicoff  _et al._ (2006) M. Chernicoff, J. A. Garcia, and A. Guijosa, JHEP 09, 068 (2006), arXiv:hep-th/0607089 .
  * Albacete  _et al._ (2008) J. L. Albacete, Y. V. Kovchegov, and A. Taliotis, Phys. Rev. D 78, 115007 (2008), arXiv:0807.4747 [hep-th] .
  * Boschi-Filho  _et al._ (2006) H. Boschi-Filho, N. R. F. Braga, and C. N. Ferreira, Phys. Rev. D 74, 086001 (2006), arXiv:hep-th/0607038 .
  * Avramis  _et al._ (2007) S. D. Avramis, K. Sfetsos, and D. Zoakos, Phys. Rev. D 75, 025009 (2007), arXiv:hep-th/0609079 .
  * Yang and Yuan (2015) Y. Yang and P.-H. Yuan, JHEP 12, 161 (2015), arXiv:1506.05930 [hep-th] .
  * Nata Atmaja  _et al._ (2015) A. Nata Atmaja, H. Abu Kassim, and N. Yusof, Eur. Phys. J. C 75, 565 (2015), arXiv:1111.7045 [hep-th] .
  * Wu  _et al._ (2019) Y. Wu, D. Hou, and H.-c. Ren, Nucl. Phys. B 938, 351 (2019), arXiv:1401.3635 [hep-ph] .
  * Casalderrey-Solana  _et al._ (2014b) J. Casalderrey-Solana, H. Liu, D. Mateos, K. Rajagopal, and U. A. Wiedemann,  _Gauge/String Duality, Hot QCD and Heavy Ion Collisions_ (Cambridge University Press, 2014) arXiv:1101.0618 [hep-th] .
  * Liu  _et al._ (2006) H. Liu, K. Rajagopal, and U. A. Wiedemann, Phys. Rev. Lett. 97, 182301 (2006), arXiv:hep-ph/0605178 .
  * Armesto  _et al._ (2006) N. Armesto, J. D. Edelstein, and J. Mas, JHEP 09, 039 (2006), arXiv:hep-ph/0606245 .
  * Rougemont (2020) R. Rougemont, Phys. Rev. D 102, 034009 (2020), arXiv:2002.06725 [hep-ph] .
  * Zhu  _et al._ (2019) Z.-R. Zhu, S.-Q. Feng, Y.-F. Shi, and Y. Zhong, Phys. Rev. D 99, 126001 (2019), arXiv:1901.09304 [hep-ph] .
  * Li  _et al._ (2016) S. Li, K. A. Mamo, and H.-U. Yee, Phys. Rev. D 94, 085016 (2016).
  * Panero  _et al._ (2014) M. Panero, K. Rummukainen, and A. Schäfer, Phys. Rev. Lett. 112, 162001 (2014), arXiv:1307.5850 [hep-ph] .
  * Giataganas (2012) D. Giataganas, JHEP 07, 031 (2012), arXiv:1202.4436 [hep-th] .
  * Chernicoff  _et al._ (2012) M. Chernicoff, D. Fernandez, D. Mateos, and D. Trancanelli, JHEP 08, 041 (2012), arXiv:1203.0561 [hep-th] .
  * D’Eramo  _et al._ (2011) F. D’Eramo, H. Liu, and K. Rajagopal, Phys. Rev. D 84, 065015 (2011), arXiv:1006.1367 [hep-ph] .
  * Bitaghsir Fadafan  _et al._ (2011) K. Bitaghsir Fadafan, B. Pourhassan, and J. Sadeghi, Eur. Phys. J. C 71, 1785 (2011), arXiv:1005.1368 [hep-th] .
  * Li  _et al._ (2014) D. Li, J. Liao, and M. Huang, Phys. Rev. D 89, 126006 (2014), arXiv:1401.2035 [hep-ph] .
  * Zhang and Zhu (2019) Z.-q. Zhang and X. Zhu, Eur. Phys. J. C 79, 107 (2019).
  * Liu  _et al._ (2007) H. Liu, K. Rajagopal, and U. A. Wiedemann, JHEP 03, 066 (2007), arXiv:hep-ph/0612168 .
  * Buchel (2006) A. Buchel, Phys. Rev. D 74, 046006 (2006), arXiv:hep-th/0605178 .
  * Nakano  _et al._ (2007) E. Nakano, S. Teraguchi, and W.-Y. Wen, Phys. Rev. D 75, 085016 (2007), arXiv:hep-ph/0608274 .
  * Zhu  _et al._ (2021) Z.-R. Zhu, J.-X. Chen, X.-M. Liu, and D. Hou, (2021), arXiv:2109.02366 [hep-ph] .
  * Park and Sin (1998) C. Park and S.-J. Sin, Phys. Lett. B 444, 156 (1998), arXiv:hep-th/9807156 .
  * Li and Lin (2018) S.-w. Li and S. Lin, Phys. Rev. D 98, 066002 (2018), arXiv:1711.06365 [hep-th] .
  * Liu and Tseytlin (1999) H. Liu and A. A. Tseytlin, Nucl. Phys. B 553, 231 (1999), arXiv:hep-th/9903091 .
  * Zhang  _et al._ (2018) Z.-q. Zhang, Z.-j. Luo, and D.-f. Hou, Nucl. Phys. A 974, 1 (2018), arXiv:1804.05517 [hep-th] .
  * Zhang  _et al._ (2017b) Z.-q. Zhang, D.-f. Hou, and G. Chen, J. Phys. G 44, 115001 (2017b), arXiv:1710.06579 [hep-th] .
  * Shahkarami  _et al._ (2018) L. Shahkarami, M. Dehghani, and P. Dehghani, Phys. Rev. D 97, 046013 (2018), arXiv:1511.07986 [hep-th] .
  * Adamczyk  _et al._ (2017) L. Adamczyk _et al._ (STAR), Nature 548, 62 (2017), arXiv:1701.06657 [nucl-ex] .
  * Abelev  _et al._ (2007) B. I. Abelev _et al._ (STAR), Phys. Rev. C 76, 024915 (2007), [Erratum: Phys.Rev.C 95, 039906 (2017)], arXiv:0705.1691 [nucl-ex] .
  * Fujimoto  _et al._ (2021) Y. Fujimoto, K. Fukushima, and Y. Hidaka, Phys. Lett. B 816, 136184 (2021), arXiv:2101.09173 [hep-ph] .
  * Chernodub (2021) M. N. Chernodub, Phys. Rev. D 103, 054027 (2021), arXiv:2012.04924 [hep-ph] .
  * Chernodub and Gongyo (2017) M. N. Chernodub and S. Gongyo, JHEP 01, 136 (2017), arXiv:1611.02598 [hep-th] .
  * Zhou  _et al._ (2021) J. Zhou, X. Chen, Y.-Q. Zhao, and J. Ping, Phys. Rev. D 102, 126029 (2021).
  * Golubtsova  _et al._ (2021) A. A. Golubtsova, E. Gourgoulhon, and M. K. Usova, (2021), arXiv:2107.11672 [hep-th] .
  * McInnes (2018) B. McInnes, (2018), arXiv:1808.00648 [hep-th] .
  * McInnes (2020) B. McInnes, Nucl. Phys. B 953, 114951 (2020), arXiv:1812.07146 [hep-th] .
  * Adam  _et al._ (2018) J. Adam _et al._ (STAR), Phys. Rev. C 98, 014910 (2018), arXiv:1805.04400 [nucl-ex] .
  * Ebihara  _et al._ (2017) S. Ebihara, K. Fukushima, and K. Mameda, Phys. Lett. B 764, 94 (2017), arXiv:1608.00336 [hep-ph] .
  * Becattini and Lisa (2020) F. Becattini and M. A. Lisa, Ann. Rev. Nucl. Part. Sci. 70, 395 (2020), arXiv:2003.03640 [nucl-ex] .
  * McInnes (2016) B. McInnes, Nucl. Phys. B 911, 173 (2016), arXiv:1604.03669 [hep-th] .
  * Zhang  _et al._ (2016) Z.-q. Zhang, D.-f. Hou, and G. Chen, Eur. Phys. J. A 52, 357 (2016), arXiv:1607.03985 [hep-ph] .
  * Bravo Gaete  _et al._ (2017) M. Bravo Gaete, L. Guajardo, and M. Hassaine, JHEP 04, 092 (2017), arXiv:1702.02416 [hep-th] .
  * Erices and Martinez (2018) C. Erices and C. Martinez, Phys. Rev. D 97, 024034 (2018), arXiv:1707.03483 [hep-th] .
  * Chen  _et al._ (2020) X. Chen, L. Zhang, D. Li, D. Hou, and M. Huang, (2020), arXiv:2010.14478 [hep-ph] .
  * Gibbons  _et al._ (1996) G. W. Gibbons, M. B. Green, and M. J. Perry, Phys. Lett. B 370, 37 (1996), arXiv:hep-th/9511080 .
  * Kehagias and Sfetsos (1999) A. Kehagias and K. Sfetsos, Phys. Lett. B 456, 22 (1999), arXiv:hep-th/9903109 .
  * Gwak  _et al._ (2012) B. Gwak, M. Kim, B.-H. Lee, Y. Seo, and S.-J. Sin, Phys. Rev. D 86, 026010 (2012), arXiv:1203.4883 [hep-th] .
  * Shahkarami and Charmchi (2019) L. Shahkarami and F. Charmchi, Eur. Phys. J. C 79, 343 (2019), arXiv:1904.09806 [hep-th] .
  * Gubser (2007) S. S. Gubser, Phys. Rev. D 76, 126003 (2007), arXiv:hep-th/0611272 .
  * Edelstein and Salgado (2008) J. D. Edelstein and C. A. Salgado, AIP Conf. Proc. 1031, 207 (2008), arXiv:0805.4515 [hep-th] .