aainstitutetext: Mandelstam Institute for Theoretical Physics, Witwatersrand
University, Johannesburg, South Africa.bbinstitutetext: International Centre
for Theoretical Sciences, Hesaraghatta, Bangalore, India.

# A study of non-linear Langevin dynamics under non-Gaussian noise with
quartic cumulant

Chandan Jana<EMAIL_ADDRESS>
###### Abstract

We consider a non-linear Langevin equation in presence of non-Gaussian noise
originating from non-linear bath. We claim, the parameters in the Langevin
equation are not physical. The physical parameters are obtained from a path-
integral description of the system, where the Langevin parameters are related
to the physical parameters by renormalisation flow equations. Then we compute
both numerically and analytically the velocity two point function and show
that it saturates to the bath temperature even in presence of non-linearity.
We also find the velocity four point function numerically and show that it
saturates to the analytically evaluated thermal velocity four point function
when the non-linear FDR Chakrabarty:2019qcp  is satisfied. When the non-linear
FDR, which is a manifestation of time reversibility of thermal bath, is
violated then the system does not seem to thermalise. Rather, its velocity
four point function settles to a steady state.

##  1 Introduction

A Brownian particle interacting with its thermal bath is successfully
described by linear Langevin equation. The particle satisfies (linear)
fluctuation-dissipation relation (FDR) which states that the variance of bath
fluctuation is directly proportional to linear drag of the bath. The linear
Langevin equation can be derived from a microscopic theory (Caldeira-Leggett
model Caldeira1982PathIA ) of the particle interacting with a harmonic bath by
a bilinear interaction. Adding a weak non-linear system-bath interaction
Chakrabarty:2019qcp ; Chakrabarty:2018dov  in the microscopic description is
expected to modify the linear Langevin equation by certain non-linear terms.

The authors of Chakrabarty:2019qcp  consider a particle interacting with a
thermal, harmonic bath at inverse temperature $\beta$ via a cubic interaction
of the form: $\sum_{i,j}\lambda_{ij}\,q\,x_{i}x_{j}$. Here, $q$ and $x_{i}$
are system and bath degrees of freedom respectively, coupled via
$\lambda_{ij}$. They integrate out the bath perturbatively assuming small
$\lambda_{ij}$ in Schwinger-Keldysh (SK) path-integral formalism
Schwinger:1960qe ; Keldysh:1964ud 111A textbook treatment of Schwinger-Keldysh
prescription and its application to non-equilibrium systems can be found in
Calzetta:1986ey ; Kamenev . and find a quartic effective action for the
particle. Along with the linear FDR, they find a non-linear FDR at the level
of non-linear effective parameters.

Identical FDRs are also derived by integrating out anharmonic, strongly
coupled bath (in contrast to harmonic bath) for a particle Chakrabarty:2019aeu
and scalar field Jana:2020vyx . More specifically, the baths in these models
are conformal field theories (CFT) at non-zero temperature,222A conformal
field theory (CFT) at zero temperature does not have any scale. But, a CFT at
non-zero temperature has only one scale which is the temperature itself.333A
strongly interacting conformal bath is in general hard to integrate out. We
circumvent this by using the duality between the anti-de Sitter (AdS) space
and conformal field theory (CFT) Aharony:1999ti . According to AdS/CFT duality
the strongly interacting CFT bath is mapped to a weakly coupled AdS bath. So,
we can integrate out the bath in AdS perturbatively using the gravitational SK
prescription Glorioso:2018mmw ; Chakrabarty:2019aeu ; Jana:2020vyx .
interacting linearly with the system. So, we notice that the FDRs do not
depend on specific details of the bath - whether a bath being harmonic or
strongly interacting CFT. The FDRs are also not restricted to appear only at
the level of quartic effective action; those are more generic. One can model
CFT baths Jana:2020vyx  such that any specific higher degree interaction terms
such as $q^{n}$ are generated in contrast to only quartic terms as found in
Chakrabarty:2019qcp . Surprisingly, the parameters corresponding to higher
degree interaction terms follow identical non-linear FDR. Hence the non-linear
FDR appears to be fairly generic; it only requires the bath to be thermal and
time reversal invariant Berges:2000ur ; Chakrabarty:2019aeu ; Jana:2020vyx ;
Chakrabarty:2019qcp . Therefore, it is worthwhile to test the implication of
the FDRs in thermalisation. In particular, whether the velocity variance for a
Brownian particle in a non-linear theory still saturates to the temperature
(like in a linear theory by assuming the Boltzmann constant $k_{B}=1$ in
suitable unit) or whether there is a systematic violation. Moreover, how are
the higher point velocity correlations affected by the FDRs?

The non-linear effective theory and FDRs can have potential application in
cosmology. Small non-Gaussian features (non-zero higher point correlation) in
the cosmic microwave background (CMB) have been estimated (Maldacena:2002vr ;
ALLEN198766 ; Gangui:1993tt ; Acquaviva:2002ud ; PhysRevD.46.4232  and
references therein) using inflationary cosmological models. The system field
in these models is in a thermal bath of field(s) which live outside the
cosmological horizon. The non-Gaussian features can be mapped to a non-linear
stochastic problem 10.1007/3-540-16452-9_6  and consequently the non-linear
FDR could play an important role for a better interpretation of the CMB data.
We leave it for future work. The FDRs may also appear in the effective theory
of gravitating (compact, astrophysical) objects Goldberger:2004jt ;
Galley:2009px . The effective theory can be thought of as an open system in a
bath of radiation. Since the radiation gravitons interact non-linearly with
the compact object, we expect non-linear FDR to emerge. Similar to the above
applications, many other weak non-linear open systems are candidates to
possess the FDRs.

In this work, we study the effect of the FDRs in correlation functions by
considering a simple quartic model of a non-linear Brownian particle. The
quartic, position-translation invariant SK effective action
Chakrabarty:2019qcp ; Chakrabarty:2019aeu  of the particle is given by

$\begin{split}S_{SK}=-\int
dt\left\\{q_{d}\ddot{q}+\gamma^{r}q_{d}\dot{q}-i\frac{f^{r}}{2}q_{d}^{2}-i\zeta_{\eta}^{r}\frac{q_{d}^{4}}{4!}-\zeta_{\gamma}^{r}q_{d}^{3}\dot{q}\right\\}\,,\end{split}$
(1)

where $q,\,q_{d}$ are the Schwinger-Keldysh classical and fluctuation degrees
of freedom respectively. The translation invariance in position allows only
the time derivative of $q$ terms in the action.
$\gamma^{r},f^{r},\zeta_{\eta}^{r},\zeta_{\gamma}^{r}$ are the effective
parameters emerged from the system-bath interaction. The significance of
$r$-label on the parameters will become clear shortly.

The SK path integral with action (1) is dual to a classical Langevin equation
following Martin-Siggia-Rose (MSR) Martin:1973zz ; DeDominicis:1977fw ;
article  prescription (discussed in §2). The dual Langevin equation is given
by Chakrabarty:2019qcp ; Chakrabarty:2019aeu

$\mathcal{E}(\dot{q},\eta)\equiv\ddot{q}+\left(\gamma+\zeta_{\gamma}\eta^{2}\right)\dot{q}-f\,\eta=0\,.$
(2)

Here $q$ is the particle’s position. The parameters $\gamma$ and $f$ are
damping coefficient and strength of the thermal noise respectively. The
damping coefficient is corrected by $\zeta_{\gamma}\eta^{2}$ \- a thermal
jitter in $\gamma$ with strength $\zeta_{\gamma}$. The duality also ensures
that the noise is drawn from a non-Gaussian distribution as the following.

$P(\eta)\sim e^{-\delta
t\left\\{\frac{f}{2}\eta^{2}+\frac{\zeta_{\eta}}{4!}\eta^{4}\right\\}}\,,$ (3)

$\zeta_{\eta}$ being the strength of quartic noise. The above distribution is
solely determined by the system-bath interaction; adding potential terms to
(2) does not affect the distribution. This means, e.g., if we study the
Langevin equation (2) numerically with a Gaussian noise, the particle should
not thermalise (see figure 5). So, the non-Gaussian noise plays a role in
particle’s correlations functions.

Note that the parameters in (2) and (3) are not same as those appearing in
(1). We show in §2 that the Langevin parameters are related to the respective
$r$-labelled SK parameters by a renormalisation flow, where the SK parameters
are the renormalised/physical parameters and the Langevin parameters are the
bare parameters. The renormalised parameters in the SK description follow
certain relations among themselves as the following:

$\gamma^{r}=\frac{\beta}{2}f^{r}\qquad\text{and}\qquad\zeta_{\gamma}^{r}=-\frac{\beta}{12}\zeta_{\eta}^{r}\,.$
(4)

The first relation can be identified as the linear FDR between the
renormalised damping coefficient $\gamma^{r}$ and the renormalised strength of
the noise $f^{r}$. $\beta$ is the inverse temperature of the bath. The second
relation is a non-linear generalisation of the linear FDR, thus can be called
a non-linear FDR Wang:1998wg ; Chakrabarty:2019qcp ; Chakrabarty:2019aeu ;
Jana:2020vyx . It states that the renormalised coefficient
$\zeta_{\gamma}^{r}$ is linearly related to the renormalised non-Gaussian
parameter $\zeta_{\eta}^{r}$. The linear FDR requires the bath to be thermal
whereas the non-linear FDR also needs the thermal bath to be time reversal
invariant Chakrabarty:2019qcp ; Chakrabarty:2018dov . The explicit numerical
value of $\zeta_{\eta}^{r}$ as derived in Chakrabarty:2019aeu ; Jana:2020vyx
turns out to be negative which makes the noise distribution ill defined. So,
we regulate the noise distribution (§4.1) by a cut-off.

In the first part of this work, we use the SK path integral to compute the
velocity two and four point functions analytically in §3. Since (1) is an
interacting SK action, we choose $\zeta^{r}_{\eta},\zeta^{r}_{\gamma}$ to be
sufficiently small to implement perturbation theory. However, due to lack of
available analytic tools, we evaluate only thermal correlations analytically
and match against numeric.

On the other side, we consider the Langevin description for a numerical study
of velocity correlations in §4. In the Langevin description, first we prepare
the non-Gaussian noise (3) using a numerical tool called the rejection
sampling vonNeumann1951 . According to this tool one can prepare an arbitrary
noise distribution from a known distribution. Given the required non-Gaussian
noise distribution, we solve (2) numerically and find that the velocity
variance for a non-linear Brownian particle saturates to $\frac{1}{\beta}$ \-
the bath temperature. We expect that the connected velocity four point
function (equal time) to be non-zero due to presence of non-linearity. The
numerically evaluated four point function is indeed non-zero and it saturates
after some characteristic time. We show that it thermalises (by checking
against analytic result) only if the non-linear FDR is valid. Otherwise, the
particle reaches a steady state2002cond.mat..2501M ; Harris_2007 ;
Parrondo_2009 ; Andrieux_2009 ; PhysRevLett.101.090602 . This means, thermal
nature of the bath is not sufficient to guarantee thermalisation of a non-
linear Brownian particle; it also requires time-reversal
invarianceBerges:2000ur .

##  2 Renormalisation of parameters in the discretised problem

The relations between the SK and the Langevin parameters will be central in
matching analytic correlations against the numeric ones. So, we begin by
deriving those relations following the Martin-Siggia-Rose (MSR) prescription.
The MSR prescription establishes a duality between SK path-integral and
stochastic description of a system. This means, each Langevin equation has a
SK path integral interpretation. The parameters appearing in the SK
description are physical, thus we call those as renormalised parameters. The
parameters appearing in the Langevin equation are related to renormalised ones
by renormalisation flow equations. We call those as the bare parameters. In
this section, we find relations among the bare and renormalised parameters.

The duality between the Langevin and path integral description is established
by the following steps

  * •

First we write an identity for the noise path-integral:
$1=\mathcal{Z}=\int\mathcal{D}\eta\,P[\eta]$ where $\eta$ is the noise and
$P[\eta]$ is the distribution functional for $\eta$.

  * •

Then we introduce the Langevin equation by a Dirac-Delta function as

$\mathcal{Z}=\int\mathcal{D}\eta\mathcal{D}q\,P[\eta]\,\text{det}\left[\frac{\delta\mathcal{E}}{\delta
q}\right]\,\delta(\mathcal{E}(q,\eta)),$

where $\mathcal{E}(q,\eta)=0$ is the Langevin equation.

  * •

We write the exponential form of the Dirac-Delta function by introducing an
auxiliary variable $q_{d}$ as the following.

$\delta(\mathcal{E})=\int\mathcal{D}q_{d}\,e^{-i\int dt\,\mathcal{E}q_{d}}.$

  * •

The noise $\eta$ is then integrated out preceded by a shift:
$\eta\rightarrow\eta+iq_{d}$.

  * •

Finally we identify $q$ and $q_{d}$ to the SK classical and fluctuation
degrees of freedom respectively. The final path-integral is identified to the
SK path-integral which is unity in absence of any sources.

Following the above steps of MSR rules, the path integral of the non-linear
Langevin particle (2) with noise distribution (3) is given by

$\begin{split}1=\mathcal{Z}=&\int\mathcal{D}q\,\mathcal{D}\eta\,\det\left[\frac{\delta\mathcal{E}}{\delta
q}\right]\,P[\eta]\,\delta\left(\mathcal{E}(t)\right)\\\
=&\int\mathcal{D}q\,\mathcal{D}\eta\,\mathcal{D}q_{d}\,\det\left[\frac{\delta\mathcal{E}}{\delta
q}\right]\,e^{-\int
dt\left\\{\left[\frac{f}{2}\eta^{2}+\frac{\zeta_{\eta}}{4!}\eta^{4}\right]+iq_{d}\mathcal{E}(t)\right\\}}\,.\end{split}$
(5)

Here $q_{d}$ is the difference/fluctuation d.o.f. in SK prescription. We
follow the Stratonovich conventioncite-key1 ; Stratonovich1992NonlinearNT  and
write the determinant in (5) as a path integral over two Grassmannian
variables ($\bar{c},c$) as the following.444One could have also used the Ito
convention cite-key ; Kamenev  and set the determinant to unity. In doing so
we will end up making the $\zeta_{\gamma}$ term non-local Chaudhuri:2020sc  by
one time step $\delta t$. Correspondingly, the path integral description
becomes non-local which is cumbersome to handle.

$\begin{split}\det\left[\frac{\delta\mathcal{E}}{\delta
q}\right]=&\int\mathcal{D}\bar{c}\,\mathcal{D}c\,\,e^{-i\int
dt\,\bar{c}\left[\frac{d}{dt}+\gamma+\zeta_{\gamma}\eta^{2}\right]c}\,.\end{split}$
(6)

We give a shift to the noise by $\eta\rightarrow\eta+iq_{d}$. The $q_{d}$ can
be thought of as a dual variable corresponding to the noise and is identified
to the fluctuation d.o.f. in SK description. The shifting does not essentially
change the path integral:

$\begin{split}\mathcal{Z}=&\int\mathcal{D}q\,\mathcal{D}\eta\,\mathcal{D}q_{d}\,\mathcal{D}\bar{c}\mathcal{D}c\,\,e^{-\int
dt\left\\{\left[\frac{f}{2}(\eta+iq_{d})^{2}+\frac{\zeta_{\eta}}{4!}(\eta+iq_{d})^{4}\right]+iq_{d}\mathcal{E}(t)+i\,\bar{c}\left[\frac{d}{dt}+\gamma+\zeta_{\gamma}(\eta+iq_{d})^{2}\right]c\right\\}}\,.\end{split}$
(7)

Then we integrate out the noise by computing noise-noise correlation
functions. Since there is non-Gaussianity in the noise distribution, we can
only integrate out the noise perturbatively. The path integral in terms of
noise correlation functions takes the following form in perturbation theory.

$\begin{split}\mathcal{Z}=&\int\mathcal{D}q\mathcal{D}q_{d}\mathcal{D}\bar{c}\mathcal{D}c\,e^{-\int
dt\left[\frac{f}{2}q_{d}^{2}+\frac{\zeta_{\eta}}{4!}q_{d}^{4}+iq_{d}\left(\ddot{q}+(\gamma-\zeta_{\gamma}q_{d}^{2})\dot{q}\right)+i\,\bar{c}\left(\frac{d}{dt}+\gamma-\zeta_{\gamma}q_{d}^{2}\right)c\right]}\\\
&\Bigg{\\{}1+\int dt\left[\frac{\zeta_{\eta}}{2!2!}\langle\eta^{2}\rangle
q_{d}^{2}-i\zeta_{\gamma}\langle\,\eta^{2}\rangle
q_{d}\dot{q}-i\zeta_{\gamma}\langle\eta^{2}\rangle\bar{c}c\right]+\frac{i\zeta_{\gamma}\zeta_{\eta}}{4!}\int
dt\,dt^{\prime}\,\langle\eta^{4}(t^{\prime})\eta^{2}(t)\rangle
q_{d}(t)\dot{q}(t)\\\ &-\frac{\zeta_{\eta}^{2}}{4!}\int
dt\,dt^{\prime}\left[\frac{1}{2^{2}}\langle\eta^{2}(t)\eta^{4}(t^{\prime})\rangle\,q_{d}^{2}(t)+\frac{1}{3}\langle\eta^{3}(t)\eta^{3}(t^{\prime})\rangle\,q_{d}(t)q_{d}(t^{\prime})\right]\\\
&+\frac{\zeta_{\eta}^{2}}{4!}\int
dt\,dt^{\prime}\left[\frac{2}{3}\langle\eta(t)\eta^{3}(t^{\prime})\rangle\,q_{d}^{3}(t)q_{d}(t^{\prime})+\frac{3}{4}\langle\eta^{2}(t)\eta^{2}(t^{\prime})\rangle\,q_{d}^{2}(t)q_{d}^{2}(t^{\prime})\right]\\\
&-i\zeta_{\eta}\zeta_{\gamma}\int
dt\,dt^{\prime}\left[\frac{1}{3}\langle\eta(t)\eta^{3}(t^{\prime})\rangle\,q_{d}^{2}(t)\dot{q}(t)q_{d}(t^{\prime})+\frac{1}{4}\langle\eta^{2}(t)\eta^{2}(t^{\prime})\rangle\,q_{d}^{2}(t)q_{d}(t^{\prime})\dot{q}(t^{\prime})\right]\\\
&-i\zeta_{\gamma}\zeta_{\eta}\int
dt\,dt^{\prime}\left[\frac{1}{4}\langle\eta^{2}(t)\eta^{2}(t^{\prime})\rangle\,\bar{c}(t)c(t)q_{d}^{2}(t^{\prime})+\frac{1}{3}\langle\eta(t)\eta^{3}(t^{\prime})\rangle\,\bar{c}(t)c(t)q_{d}(t)q_{d}(t^{\prime})\right]+\cdots\Bigg{\\}}\,.\end{split}$
(8)

Note that, in the Taylor expanded form, we have kept terms till quartic
order(in $q,q_{d},\bar{c},c$). The higher order terms are suppressed due to
small non-linear couplings. The noise correlation functions are UV divergent;
thus we put a UV cut-off $\delta t$ in time which also acts as the smallest
time step in the discretised problem discussed through out this paper.

We evaluate various noise correlation functions of (8) in appendix A. With the
noise correlation function, we exponentiate the RHS of (8) back and find the
following effective path integral.

$\begin{split}\mathcal{Z}=\int\mathcal{D}q\,\mathcal{D}q_{d}\,\mathcal{D}\bar{c}\,\mathcal{D}c\times\exp[-\int
dt\Bigg{\\{}iq_{d}\ddot{q}+i\gamma^{r}q_{d}\dot{q}+\frac{f^{r}}{2}q_{d}^{2}+i\bar{c}\left(\frac{d}{dt}+\gamma-\zeta_{\gamma}q_{d}^{2}\right)c\\\
+\zeta_{\eta}^{r}\frac{q_{d}^{4}}{4!}-i\zeta_{\gamma}^{r}q_{d}^{3}\dot{q}+\cdots\Bigg{\\}}\Bigg{]}\,.\end{split}$
(9)

The $r$ labelled parameters are renormalised parameters (evaluated in appendix
A) and are given by

$\displaystyle f^{r}$ $\displaystyle=f-\frac{\zeta_{\eta}}{2f\delta
t}+\frac{2\zeta_{\eta}^{2}}{3f^{3}(\delta t)^{2}}+\cdots\,,$ (10a)
$\displaystyle\gamma^{r}$ $\displaystyle=\gamma+\frac{\zeta_{\gamma}}{f\delta
t}-\frac{\zeta_{\gamma}\zeta_{\eta}}{2f^{3}(\delta t)^{2}}+\cdots\,,$ (10b)
$\displaystyle\zeta_{\eta}^{r}$
$\displaystyle=\zeta_{\eta}-\frac{7}{2}\frac{\zeta_{\eta}^{2}}{f^{2}\delta
t}+\frac{149}{12}\frac{\zeta_{\eta}^{3}}{f^{4}(\delta t)^{2}}+\cdots\,,$ (10c)
$\displaystyle\zeta_{\gamma}^{r}$
$\displaystyle=\zeta_{\gamma}-\frac{3\zeta_{\eta}\zeta_{\gamma}}{2f^{2}\delta
t}+\frac{10}{3}\frac{\zeta_{\eta}^{2}\zeta_{\gamma}}{f^{4}(\delta
t)^{2}}+\cdots\,.$ (10d)

Here $\delta t$ is the UV cut-off and the smallest time step and the
perturbation parameter is $\sim\frac{\zeta_{\eta}}{f^{2}\delta t}$. Note that
the front factors, in RHS of the above relations, are larger down the series.
So, we choose the bare parameters such that perturbation theory is valid. All
of the renormalised parameters are not independent. The linear renormalised
parameters $\\{f^{r},\gamma^{r}\\}$, as discussed in the introduction, are
related by the linear FDR:

$\gamma^{r}=\frac{\beta}{2}f^{r}\,.$ (11)

If the bath also possesses a time reversal symmetry, the non-linear
renormalised parameters are related by the non-linear FDR:

$\zeta_{\gamma}^{r}=-\frac{\beta}{12}\zeta_{\eta}^{r}\,.$ (12)

The parameters corresponding to the Grassmann variables $\bar{c},c$ also
obtain correction from noise correlation function, but those can be neglected
due to the following reason. Since, $\bar{c}$ and $c$ do not constitute any
physical observable due to their Grassmannian nature, they do not have any
direct contribution in any correlation function; they only have sub-leading
loop contribution. So, the leading parameters corresponding to Grassmann
vertices in (9) make sub-leading contribution to correlation function. Hence
noise loop correction to those parameters will be sub-sub-leading and thus are
ignored (see appendix B).

One might be concerned about LHS of (10) being called as renormalised
parameters, since they may receive correction due to quartic terms in (9) and
be further renormalised. We will argue in the next section that the
renormalised parameters are not further corrected by $\zeta_{\gamma}^{r}$ or
$\zeta_{\eta}^{r}$ interactions. Although we get $\bar{c}$-$c$ loop
corrections but those are sub-leading as discussed above. Therefore we can
safely call LHS of (10) as the renormalised parameters.

The bare parameters in (10), according to the notion of renormalisation,
should flow with $\delta t$ whereas the renormalised parameters are fixed by
the system itself. This means that the bare parameters can take different
values along the flow depending on $\delta t$. One can find the flow equations
by writing the bare parameters in terms of the renormalised ones, which can be
found easily by making suitable approximation.

##  3 Thermal velocity correlations in Schwinger-Keldysh formalism

Correlation function for a generic initial condition is hard to compute
analytically. In this section we rather compute thermal correlations such as
velocity variance and equal time velocity four point correlation analytically
using the SK path integral.

###  3.1 Schwinger-Keldysh propagators and velocity variance

We evaluate the thermal SK propagators (explicitly in appendix B) to compute
the velocity correlation functions. The SK propagators are given by

$\displaystyle iG^{K}(t_{2},t_{1})\equiv$
$\displaystyle\langle\dot{q}(t_{2})\dot{q}(t_{1})\rangle=\frac{f^{r}}{2\gamma^{r}}e^{-\gamma^{r}|t_{2}-t_{1}|}\,,$
(13) $\displaystyle iG^{A}(t_{2},t_{1})\equiv$ $\displaystyle\langle
q_{d}(t_{2})\dot{q}(t_{1})\rangle=-i\,\Theta(t_{1}-t_{2})e^{-\gamma^{r}(t_{1}-t_{2})}\,,$
(14) $\displaystyle iG^{R}(t_{2},t_{1})\equiv$
$\displaystyle\langle\dot{q}(t_{2})q_{d}(t_{1})\rangle=-i\,\Theta(t_{2}-t_{1})e^{-\gamma^{r}(t_{2}-t_{1})}\,,$
(15) $\displaystyle iG^{d}(t_{2},t_{1})\equiv$ $\displaystyle\langle
q_{d}(t_{2})q_{d}(t_{1})\rangle=0\,.$ (16)

Here, $\dot{q}$ is the particle’s velocity. The above non-zero propagators are
the Keldysh, advanced and the retarded propagators respectively. We depict the
propagators in figure 1.

$iG^{K}(t_{2},t_{1})\,:$$t_{1}$$t_{2}$$iG^{A}(t_{2},t_{1})\,:$$t_{1}$$t_{2}$$iG^{R}(t_{2},t_{1})\,:$$t_{1}$$t_{2}$$iG^{d}(t_{2},t_{1})\,:$$t_{1}$$t_{2}$
Figure 1: Diagrammatic representation of Schwinger-Keldysh propagators: The
explicit expression for the propagators are the following.
$iG^{K}(t_{2},t_{1})=\langle\dot{q}(t_{2})\dot{q}(t_{1})\rangle\,,iG^{A}(t_{2},t_{1})=\langle
q_{d}(t_{2})\dot{q}(t_{1})\rangle\,,iG^{R}(t_{2},t_{1})=\langle\dot{q}(t_{2})q_{d}(t_{1})\rangle$
and $iG^{d}(t_{2},t_{1})=\langle q_{d}(t_{2})q_{d}(t_{1})\rangle$. The solid
and dashed lines correspond to $\dot{q}$ and $q_{d}$ respectively.

Note that the propagator in (16) is zero. In absence of this propagator and
with the available interaction terms in (9), there will not be any loop
correction consisting of $q_{d},\dot{q}$ to velocity variance. Although there
will be Grassmann loop corrections, those are negligibly small, thus can be
ignored (see appendix B). Therefore, in absence of any non-linear correction,
the Keldysh propagator (13) is the velocity two point function and at equal
time is the velocity variance. Thus the velocity variance, by using the linear
FDR, is given by

$\langle\dot{q}^{2}(t)\rangle=\frac{f^{r}}{2\gamma^{r}}=\frac{1}{\beta}\,.$
(17)

So, we find that the velocity variance is equal to the temperature of the
bath, like in linear Langevin dynamics. Weak non-linearity does not bring in
correction to velocity variance in thermal equilibrium.

###  3.2 Microscopic time reversibility of bath and non-linear FDR

The non-linear FDR, as discussed in the introduction, seems to be generic from
the discussions of Chakrabarty:2018dov ; Chakrabarty:2019qcp  and
Chakrabarty:2019aeu . In Chakrabarty:2018dov ; Chakrabarty:2019qcp , non-
linear FDR was derived assuming weak non-linear system-bath interaction. On
the other hand, Chakrabarty:2019aeu  found an identical FDR with linear
system-bath interaction where the bath being non-linear and strongly coupled.
Thus non-linear FDR does not seem to depend on the details of bath, rather it
depends on a generic principle: microscopic time reversal invariance of the
bath.

Microscopic time reversibility of bath, in general, imposes certain
constraints on the system’s effective theory. These relations were first
discovered by Onsager Onsager:1931jfa ; Onsager:1931kxm  and extended by
Casimir RevModPhys.17.343  in the context of a system with multiple degrees of
freedom. According to this observation, if $a)$ the bath Hamiltonian $H_{B}$
commutes with the time reversal operator T and $b)$ the bath operator
$\mathcal{O}$, through which the bath interacts with the system, is even under
time reversal:

$\begin{split}\left[\textbf{T},H_{B}\right]=0\,,\qquad\textbf{T}^{-1}\,\mathcal{O}(t)\,\textbf{T}=\mathcal{O}(-t)\,,\end{split}$
(18)

then certain parameters in the system’s effective theory are related.

Following Chakrabarty:2018dov ; Chakrabarty:2019qcp , we motivate the above
assertion in the context of non-linear FDR by considering the following bath
four point function:

$\langle\mathcal{O}(t_{1})\mathcal{O}(t_{2})\mathcal{O}(t_{3})\mathcal{O}(t_{4})\rangle\,.$
(19)

This kind of four point functions would appear in the influence phase of the
effective theory when we integrate out the bath perturbatively. Due to time
reversal invariance we have

$\langle\mathcal{O}(t_{1})\mathcal{O}(t_{2})\mathcal{O}(t_{3})\mathcal{O}(t_{4})\rangle=\langle\textbf{T}^{\dagger}\mathcal{O}(t_{1})\mathcal{O}(t_{2})\mathcal{O}(t_{3})\mathcal{O}(t_{4})\textbf{T}\rangle^{\ast}$
(20)

Using time reversal (18) and hermiticity of $\mathcal{O}$, we get

$\langle\mathcal{O}(t_{1})\mathcal{O}(t_{2})\mathcal{O}(t_{3})\mathcal{O}(t_{4})\rangle=\langle\mathcal{O}(-t_{1})\mathcal{O}(-t_{2})\mathcal{O}(-t_{3})\mathcal{O}(-t_{4})\rangle^{\ast}=\langle\mathcal{O}(-t_{4})\mathcal{O}(-t_{3})\mathcal{O}(-t_{2})\mathcal{O}(-t_{1})\rangle$
(21)

Let us choose $t_{1}=t_{4}=0$ and $t_{2}=t_{3}=t$. Then we get

$\langle\mathcal{O}(0)\mathcal{O}(t)\mathcal{O}(t)\mathcal{O}(0)\rangle=\langle\mathcal{O}(0)\mathcal{O}(-t)\mathcal{O}(-t)\mathcal{O}(0)\rangle\,.$
(22)

Note that the LHS is a SK correlation function for $t>0$ since it can be
captured in SK contourKamenev ; Haehl . But the RHS does not belong to SK
correlation - it is rather an out-of-time-ordered (OTO) correlation function.
(A review on generalised OTO correlation can be found in Haehl:2017qfl ;
Chaudhuri:2018ymp ; Chaudhuri:2018ihk .) Now, we impose thermality on the RHS
via the Kubo-Martin-Sschwinger (KMS) Kubo:1957mj ; Martin:1959jp  relation and
find

$\langle\mathcal{O}(0)\mathcal{O}(t)\mathcal{O}(t)\mathcal{O}(0)\rangle=\langle\mathcal{O}(-i\beta)\mathcal{O}(0)\mathcal{O}(-t)\mathcal{O}(-t)\rangle\,.$
(23)

Thus we find that two different SK correlation functions are related by both
time-reversal invariance and thermality. When we integrate out the bath then
bath correlations at the system time scale appears as parameters in the
effective theory. Thus from (23), we expect to get certain relation between
corresponding parameters. This motivates the idea of why should there be a
non-linear FDR under time reversal invariance and thermal nature of the bath.

Note that, for our system, there is no Onsager like relations at the level of
two point correlations. The effective path integral (9) describes only one
particle with two parameters ($f^{r},\gamma^{r}$) at quadratic order.
Corresponding to these parameters there are two two-point bath correlations
Chakrabarty:2019qcp  which are related by the KMS condition. Consequently,
$f^{r}$ and $\gamma^{r}$ are related by the linear FDR (11). There are no
other parameters so that time reversal invariance would act and give non-
trivial relation(s) among those parameters. Therefore, time reversal
invariance acts trivially at the level of two point correlations for our
system.

###  3.3 Equal-time velocity four point function

We compute equal-time velocity four point function to study the implication of
the non-linear FDR (12). The connected part of the equal-time velocity four
point function is defined as

$\langle\dot{q}^{4}(t)\rangle_{c}=\langle\dot{q}(t)^{4}\rangle-3\langle\dot{q}^{2}(t)\rangle^{2}\,.$
(24)

We compute it perturbatively from the path-integral (9) using Feynman
diagrams. As Feynman rules, we need the propagators and the vertex factors.
The propagators are evaluated in appendix B and the vertex factors for
$q_{d}^{4}$ and $\dot{q}q_{d}^{3}$ vertices are given by $(-\zeta_{\eta}^{r})$
and $(i4!\zeta_{\gamma}^{r})$ respectively.

The SK diagrams contributing to $\langle\dot{q}^{4}(t)\rangle_{c}$ at leading
order, are drawn in figure 2. We find that both $\zeta_{\eta}^{r}$ and
$\zeta_{\gamma}^{r}$ contribute at the tree level. One can show that all loop
diagrams consist of only sub-leading Grassmann loops, provided the propagator
$G^{d}$ in (16) is zero. Hence, tree level contribution to (24) is the most
significant contribution.

$\left\langle(\dot{q}(t))^{4}\right\rangle_{c}\,=$$(-\zeta_{\eta}^{r})\,\times$$\dot{q}(t)$$\dot{q}(t)$$\dot{q}(t)$$\dot{q}(t)$$t^{\prime}$$+\,\,\,(i4!\zeta_{\gamma}^{r})\,\times$$\dot{q}(t)$$\dot{q}(t)$$\dot{q}(t)$$\dot{q}(t)$$t^{\prime}$
Figure 2: Leading order contribution to $\langle\dot{q}^{4}(t)\rangle_{c}$:
The vertex factors of the SK diagrams are $(-\zeta_{\eta}^{r})$ and
$(4!i\zeta_{\gamma}^{r})$ respectively. We assume an out-going time convention
which means time flows from $t^{\prime}$ to $t$. In this convention, all of
the propagators (except the Keldysh propagator on the right figure) are
retarded propagators.

To evaluate the SK diagrams, we choose an ’out-going time’ convention for the
diagrams in figure 2. According to this convention, all propagators (except
the Keldysh propagator) are retarded. The integral representation of the SK
diagrams in figure 2 is given by

$\begin{split}\langle\dot{q}^{4}(t)\rangle_{c}=&\,-\zeta_{\eta}^{r}\int_{-\infty}^{\infty}dt^{\prime}\langle\dot{q}(t)q_{d}(t^{\prime})\rangle^{4}+(4!i\zeta_{\gamma}^{r})\int_{-\infty}^{\infty}dt^{\prime}\langle\dot{q}(t)q_{d}(t^{\prime})\rangle^{3}\langle\dot{q}(t^{\prime})\dot{q}(t)\rangle\,.\end{split}$
(25)

Substituting (13) and (15) we evaluate the integrals. Then we use the linear
and non-linear FDRs (4) to get

$\langle\dot{q}^{4}(t)\rangle_{c}=-\frac{\zeta^{r}_{\eta}}{4\gamma^{r}}-4!\zeta_{\gamma}^{r}\frac{f^{r}}{2\gamma^{r}}\frac{1}{4\gamma^{r}}=\frac{\zeta^{r}_{\eta}}{4\gamma^{r}}\,.$
(26)

Note that the above quantity is a negative number since $\zeta_{\eta}^{r}$ is
negative as pointed out in the introduction. We define a quartic deviation
(analogous to standard deviation) of velocity as

$v_{QD}\equiv-|\langle\dot{q}^{4}(t)\rangle_{c}|^{1/4}\,$ (27)

and match with the numeric in figure 6.

##  4 A numerical analysis of velocity correlations

###  4.1 Preparation of noise

Figure 3: The target distribution $P$ is found from the known reference
distribution $P_{g}$. $P_{g}$ is chosen to be a scaled Gaussian distribution
whose exact form is not relevant; it only has to contain the full target
distribution. The target distribution $P(\eta)\sim e^{-\delta
t\left(\frac{f}{2}\eta^{2}+\frac{\zeta_{\eta}}{4!}\eta^{4}\right)}$. The green
dots generated from the rejection sampling, plot $u\times P_{g}$. Note that
these dots do not go out of the red curve.

We prepare the non-Gaussian noise distribution $P(\eta)$ in (3) from a
reference Gaussian distribution by a method called rejection sampling555I
thank Anugu Sumith Reddy for bringing the rejection sampling method to my
notice.vonNeumann1951  originally proposed by von Neumann. (One could use
different numerical methods such as Metropolis Hastings
doi:10.1080/00031305.1995.10476177 ). According to rejection sampling, we
choose a reference noise distribution (Gaussian) function $P_{g}(\eta)$. The
reference distribution $P_{g}(\eta)$ need not be normalised and should be such
that it contains the full target distribution $P(\eta)$ within, as shown in
figure 3. Then we choose a random number $u$ between $[0,1]$. If
$u<\frac{P(\eta)}{P_{g}(\eta)}$ for a generated $\eta$ from $P_{g}$, we keep
that $\eta$ for $P$; else we discard. By repeating this process many times we
can generate $P$. Out of 100 iteration figure 3 plots those data which follow
the above inequality. An infinite number of iteration should populate the
entire target distribution.

It is worthwhile to point out a subtlety in the target noise distribution
$P(\eta)$. Derivation of the non-linear Langevin equation from microscopic
theories Chakrabarty:2019aeu ; Jana:2020vyx  finds that the actual numerical
value of $\zeta_{\eta}$ is negative which makes $P(\eta)$ ill defined since it
diverges as $\eta\rightarrow\infty$. To stabilise it, we add a small regulator
such that the regulated distribution is given by

$P_{\text{regulated}}(\eta)\sim e^{-\delta
t\left(\frac{f}{2}\eta^{2}+\frac{\zeta_{\eta}}{4!}\eta^{4}+\epsilon\eta^{6}\right)}\,,\,\,\,1>>\epsilon>0\,.$
(28)

The justification of the above regulator could be the following. All
derivations of the Langevin equation consider only the leading order influence
phase Feynman:1963fq  from the bath. The sub-leading correction should
stabilise the distribution to make it a sensible theory. We choose $\epsilon$
such that it nullifies any sixth cumulant generated during the noise
preparation. We generate sufficient number ($10^{8}$) of $\eta$ using
rejection sampling and plot two, four and six-point connected noise
correlations. Since $\eta$ is a white noise, all of its correlations can have
peak (non-zero values) when all $\eta$ are evaluated at the same time. So, we
choose to vary one time-argument of each correlation, keeping the rest fixed
at equal times as the following:

$\langle\eta(t)\eta(t^{\prime})\rangle\,,\qquad\left\langle\eta(t)(\eta(t^{\prime}))^{3}\right\rangle_{c}\qquad\text{and}\qquad\left\langle\eta(t)(\eta(t^{\prime}))^{5}\right\rangle_{c}\,.$
(29)

The connected four and six point correlations are defined as

$\displaystyle\left\langle\eta(t)(\eta(t^{\prime}))^{3}\right\rangle_{c}$
$\displaystyle=$
$\displaystyle\left\langle\eta(t)(\eta(t^{\prime}))^{3}\right\rangle-3\left\langle(\eta(t^{\prime}))^{2}\right\rangle\langle\eta(t)\eta(t^{\prime})\rangle$
(30)
$\displaystyle\qquad\left\langle\eta(t)(\eta(t^{\prime}))^{5}\right\rangle_{c}$
$\displaystyle=$
$\displaystyle\qquad\left\langle\eta(t)(\eta(t^{\prime}))^{5}\right\rangle-5\left\langle(\eta(t^{\prime}))^{2}\right\rangle\left\langle\eta(t)(\eta(t^{\prime}))^{3}\right\rangle_{c}$
(31)
$\displaystyle-10\left\langle(\eta(t^{\prime}))^{4}\right\rangle_{c}\left\langle\eta(t)\eta(t^{\prime})\right\rangle-15\left\langle(\eta(t^{\prime}))^{2}\right\rangle^{2}\left\langle\eta(t)\eta(t^{\prime})\right\rangle$

where correlations are evaluated from the below explicit formula.

$\displaystyle\langle\eta(t)(\eta(t^{\prime}))^{n}\rangle=\frac{1}{N}\sum_{i=1}^{N}\eta_{i}(t)(\eta_{i}(t^{\prime}))^{n}\,,\qquad
N\text{ being the number of iteration.}$ (32)

We depict the above two point and connected four & six point correlations in
figure 4. We find that two and four point functions show sharp peak while the
six point function shows only fluctuation. This ensures that the noise
distribution have insignificant sixth cumulant as required. Since we have
successfully prepared the noise, we are ready to study the particle’s dynamics
numerically.

$(a)$$(b)$$(c)$ Figure 4: We compute $\langle\eta(t)\eta(1.25)\rangle$,
$\langle\eta(t)\eta(1.25)^{3}\rangle$ and
$\langle\eta(t)\eta(1.25)^{5}\rangle$ in figure $(a),\,(b)$ and $(c)$
respectively. Here, the horizontal axis is in some unit of time. So, both
$t,1.25$ are numerical value in that unit of time. Since the noise does not
have any dynamics, we get a peak (non-zero correlation) only when time of all
noise variable are equal. The parameters for non-Gaussian noise distribution
are chosen as $f=5.11$ and $\zeta_{\eta}=-0.086$ which are also the bare
parameters discussed in section 2. The units of these parameters are given in
footnote 6. The six point function in $(c)$ is essentially indistinguishable
from the fluctuation. The number of iteration for each graph is $10^{7}$.

#### Comment on the regime of parameters:

The bare parameters given in (2) are used in the numerical computation whereas
the renormalised ones are used in the analytics. The renormalised parameters
are related to the bare parameters by an infinite series as written in (10).
We evaluate till second order noise loop corrections in this series. This
means our choice of bare parameters should be such that the convergence is
fast enough to ensure matching of numeric and analytic data. We choose the
bare parameters and the time step as

$f=5.11\,,\qquad\zeta_{\eta}=-0.086\qquad\text{and}\qquad\delta t=0.05$ (33)

respectively so that the perturbation parameter
$\frac{\zeta_{\eta}}{f^{2}\delta t}(=0.066)$ in (10) makes the second order
term small (compared to the first order term). We also choose the numerical
value of temperature to be $\frac{1}{\beta}=2$. The corresponding renormalised
parameters666 We have not explicitly written the units of the parameters in
the main text. If we assume the dimensions of $\dot{q}$ and $t$ as
[length][time]-1 and [time] respectively, then the dimensions of $f$ and
$\gamma$ are found to be [length]2[time]-3 and [time]-1 respectively.
Similarly, the dimensions of $\zeta_{\eta},\,\zeta_{\gamma}$ are respectively
[length]4 [time]-5, [length]2[time]-3. The dimension of $\beta$, using the
FDRs, is given by [time]2 [length]-2 for unit mass and unit Boltzmann
constant. But, most importantly, the perturbation parameter
$\zeta_{\eta}/(f^{2}\delta t)$ is a dimensionless number. are given by

$f^{r}=5.29\,,\qquad\gamma^{r}=1.31\,,\qquad\zeta_{\eta}^{r}=-0.112\,,\qquad\zeta_{\gamma}^{r}=0.0046\,\,.$
(34)

Let us now ask the following question. How much can we deviate from (33) - the
values of individual bare parameters, keeping the perturbation parameter at
same order? If we make each parameter larger then the Euler-Maruyama error
grows due to large values of $f$ and $\delta t$. Thus we will find deviation
from analytic results. On the other hand, if we make each parameter smaller
than in (33), we observe the following. We may get rid of the Euler-Maruyama
error, but small $f$ and $\delta t$ bring large fluctuation in higher
cumulants. For example, for $f=1,\,\zeta_{\eta}=-0.001$ and $\delta t=0.5$,
$5\times 10^{8}$ number of ensemble is insufficient to tame the fluctuation
and get a definitive plot for connected velocity four point function defined
in (24). Thus we conclude that we should not make the parameters too large or
too small. However, if one gets improvements on either side by better
numerical methods, then the domain of validity of parameters can be made
bigger.

###  4.2 Velocity variance

In this section, we numerically compute the velocity variance of the non-
linear Langevin particle and show that its velocity variance matches to a high
accuracy with that of linear Langevin particle at same temperature.

In order to study the dynamics of (2), we discretise it as the following.

$\begin{split}q(t+\delta t)=&\,q(t)+\delta t\,\dot{q}\,,\\\ \dot{q}(t+\delta
t)=&\,\dot{q}(t)-\delta t[\gamma+\zeta_{\gamma}\eta^{2}(t)]\dot{q}(t)+f\delta
t\,\eta(t)\,.\end{split}$ (35)

The key point is that the parameters in the above equation are bare
parameters. With these parameters we plot variance of velocity of the particle
in figure 5. In this figure, we also simulate the linear Langevin equation
with renormalised linear couplings ($f^{r},\gamma^{r}$) and find that the
velocity variances in both cases match. One might argue, discrepancy between
them in figure 5 are not observed due to _weak_ system-bath coupling and thus
there is no significant distinction between bare and renormalised parameters.
This is not correct. Had we used bare linear parameters ($f,\gamma$) and made
$\zeta_{\eta},\zeta_{\gamma}=0$ in (2), we would have got the dashed curve.
So, the weak non-linear couplings do make a significant contribution.

Figure 5: The numerical value of bare parameters are given by
$f=5.11,\,\gamma=1.31,\,\zeta_{\eta}=-0.086,\,\zeta_{\gamma}=0.0042$ in
respective units 6. The corresponding renormalised parameters are evaluated
from (10) as $f^{r}=5.29,\,\gamma^{r}=1.32,\,\zeta_{\eta}^{r}=-0.112$ and
$\zeta_{\gamma}^{r}=0.0046$ for $\delta t=0.05$ . We assume both initial
position and momentum as zero and find that velocity variance exactly matches
in both linear and non-linear cases. Notice that saturation value does not
occur exactly at temperature $\frac{1}{\beta}=\frac{f^{r}}{2\gamma^{2}}=2$ due
to the Euler-Maruyama error. We have chosen $f,\delta t$ large enough that the
error is noticeable. However, since we are matching with linear Langevin
result which also have identical error, the apparent deviation is irrelevant.
The dashed line plots velocity variance with bare linear parameters but zero
non-linear parameters. The number of iteration for each curve is $10^{7}$.

###  4.3 Velocity four point function and non-linear FDR

We find numerically the connected velocity four point function defined in (24)
for $5\times 10^{7}$ iterations. However, due to our choice of parameters,
velocity four point function is comparable to typical fluctuation for the
above number of iterations. We can reduce the fluctuation substantially by
taking a fourth root of the velocity four point function as the following.

$\displaystyle
v_{QD}(t)\big{|}_{\text{numeric}}\equiv\begin{cases}\left(\langle\dot{q}^{4}(t)\rangle_{c}\right)^{1/4}\,\,,\qquad\text{if}\qquad\langle\dot{q}^{4}(t)\rangle_{c}>0\\\
-\left|\langle\dot{q}^{4}(t)\rangle_{c}\right|^{1/4}\,\,,\qquad\text{if}\qquad\langle\dot{q}^{4}(t)\rangle_{c}<0\end{cases}$
(36)

The above definition for numerical analysis is different from (27) due to our
prior ignorance about the sign of $\langle\dot{q}^{4}(t)\rangle_{c}$ in
numerical analysis. We plot the above quantity in figure 6 and find that the
saturation value of the numerical data matches quite well with the analytic
result when the non-linear FDR (12) is satisfied.

In figure 6, there is a bump during the initial time evolution. To understand
this feature, let us consider the discretised non-linear Langevin equation
(35) at $t=0$. Since the initial velocity is zero, there is no
$\zeta_{\gamma}$ contribution at $t=0$. But we still have a $\zeta_{\eta}$
contribution from the noise which has a positive cumulant. Thus we get a
positive bump at $t=\delta t$. After a few more time-steps $\zeta_{\gamma}$
term takes over since it has a larger negative contribution as pointed out in
(25).

Figure 6: We plot the quartic deviation $v_{QD}(t)|_{\text{numeric}}$ (36) for
$\alpha=1,2,3,4$ as defined in (37). We also plot the analytic, thermal
$v_{QD}$ as found in (26) and find that $v_{QD}$ agrees well with
$v_{QD}(t)|_{\text{numeric}}$ when $\alpha=1$, i.e., non-linear FDR is
satisfied. The numerical value of the bare parameters for $\alpha=1$ are
chosen as $f=5.11,\,\gamma=1.31,\,\zeta_{\eta}=-0.086,\,\zeta_{\gamma}=0.0042$
in respective units. For $\alpha\neq 1$, only $\zeta_{\gamma}$s are different
and are chosen as $0.0084,\,0.0126,\,0.0167$ to get $\alpha=2,\,3,\,4$
respectively. We find that fluctuation in $v_{QD}(t)|_{\text{numeric}}$
reduces as we make $\alpha$ bigger. So, we consider different number of
iteration for each plot. The numbers of iteration are $5\times 10^{7}$ for
$\alpha=1$, $3\times 10^{7}$ for $\alpha=2$ and $10^{7}$ for $\alpha=3,4$
each.

We also check how the quartic deviation behaves if we violate the non-linear
FDR. We choose different values of the bare parameter $\zeta_{\gamma}$ such
that the FDR differs by a factor of $\alpha$ defined as the following:

${\zeta_{\gamma}^{r}}_{\text{dev}}\equiv\alpha\frac{-\zeta_{\eta}^{r}}{12}\,,$
(37)

where $\alpha=1$ defines the FDR. We find that the numeric quartic deviations
differ from the analytically evaluated thermal quartic deviations (27) when
$\alpha\neq 1$, as shown in figure 6. (Note that $\alpha=1$ for the analytic
data in figure 6 as corresponds to thermal quartic deviation.) The
interpretation of this disagreement is the following. If the time reversal
invariance of bath is broken ($\alpha\neq 1$), then we do not expect the
particle to thermalise Berges:2000ur ; 2002cond.mat..2501M ; Harris_2007 ;
Parrondo_2009 ; Andrieux_2009 ; PhysRevLett.101.090602 . Rather, in our case,
it reaches a steady state.

##  5 Conclusion

This work deals with systematic analysis of non-linear Langevin equation in
presence of non-Gaussian noise. We show how the non-Gaussianity introduces the
notion of renormalisation such that the bare/unphysical parameters are the
parameters appearing in Langevin equation. On the other hand, the
renormalised/physical parameters appear in the Schwinger-Keldysh effective
action. We show, the bare and renormalised parameters are related by noise
loop correction.

In the second half of this work, we do a numerical analysis to find the
velocity variance and velocity four point function. We show that the velocity
variance for both linear and non-linear Langevin dynamics, saturates to the
temperature of the bath. We compute velocity four point function numerically
and show that it agrees with thermal velocity four point function (calculated
analytically) when non-linear FDR is satisfied. If we deviate from the non-
linear FDR, the system does not thermalise but settles to a steady state.

###### Acknowledgements.

I would like to thank Soumyadeep Chaudhuri for initial collaboration and
guidance. I thank Arghya Das, Anupam Kundu, R. Loganayagam, Ajit Mehta, Avaya
Pal, Archak Purkayastha, Anugu Sumith Reddy, Arnab Rudra, Prashant Singh,
Rahul Singh and Akhil Sivakumar for useful suggestions.

##  Appendix A Noise correlation functions in corrections of parameters

We consider the path integral in (8) which evaluates correction to all
parameters. To evaluate the loop diagrams we need the noise propagator as a
Feynman rule. Since the noise does not have any dynamics, the propagator can
only be Dirac-Delta function of time. But, in our numerical analysis, since
the smallest available time scale is the time step $\delta t$, we write the
following expression for noise propagator.

$\begin{split}\langle\eta(t_{1})\eta(t_{2})\rangle=\begin{cases}\frac{1}{f\delta
t}\qquad\text{if}\qquad t_{1}=t_{2}\\\
0\qquad\text{else}\,.\end{cases}\end{split}$ (38)

The corresponding diagrammatic representation is given in figure 7.

$t$$t$$:$$\langle\eta(t)\eta(t)\rangle=\frac{1}{f\delta t}$ Figure 7:
Diagrammatic representation of noise propagator

### Correction to $q_{d}^{2}$ vertex:

The correction to $f$ are obtained from the noise two loop diagrams drawn in
figure 8. The correction to $f$ is also pointed out below each diagram.

$\times$$-\frac{f^{r}}{2}$$=$$\times$$-\frac{f}{2}$+$\times$$\frac{1}{4}\frac{\zeta_{\eta}}{f\delta
t}$+$\times$$\times$$-\frac{1}{12}\frac{\zeta_{\eta}^{2}}{f^{3}\delta
t^{2}}$+$\times$$\times$$-\frac{1}{8}\frac{\zeta_{\eta}^{2}}{f^{3}\delta
t^{2}}$+$\times$$-\frac{1}{8}\frac{\zeta_{\eta}^{2}}{f^{3}\delta t^{2}}$
Figure 8: Loop correction to $q_{d}^{2}$ vertex in SK path integral (8). The
corrected parameter is the renormalised parameter.

Adding all the contributions we get

$f^{r}=f-\frac{\zeta_{\eta}}{2f\delta
t}+\frac{2\zeta_{\eta}^{2}}{3f^{3}(\delta t)^{2}}$ (39)

### Correction to $q_{d}\dot{q}$ vertex:

The diagrams which contribute to $q_{d}\dot{q}$ vertex are drawn in figure 9.

$\times$$-i\gamma^{r}$$=$$\times$$-i\gamma$+$\times$$\frac{-i\zeta_{\gamma}}{f\delta
t}$+$\times$$\frac{i}{2}\frac{\zeta_{\gamma}\zeta_{\eta}}{f^{3}\delta t^{2}}$
Figure 9: Loop correction to $q_{d}\dot{q}$ vertex in SK path integral (8).

Adding all contribution we get the following expression for $\gamma^{r}$.

$\gamma^{r}=\gamma+\frac{\zeta_{\gamma}}{f\delta
t}-\frac{\zeta_{\gamma}\zeta_{\eta}}{2f^{3}(\delta t)^{2}}\,.$ (40)

### Correction to $q_{d}^{4}$ vertex:

The correction to $q_{d}^{4}$ vertex gets contribution from the diagrams shown
in figure 10.

$\times$$-\frac{1}{4!}\zeta_{\eta}^{r}$$=$$\times$$-\frac{1}{4!}\zeta_{\eta}$+$\times$$\times$$\frac{3}{2\times
4!}\frac{\zeta_{\eta}^{2}}{f^{2}\delta
t}$+$\times$$\times$$\frac{2}{4!}\frac{\zeta_{\eta}^{2}}{f^{2}\delta
t}$+$\times$$\times$$-\frac{3}{2\times 4!}\frac{\zeta_{\eta}^{3}}{f^{4}\delta
t^{2}}$+$\times$$\times$$-\frac{3}{4\times
4!}\frac{\zeta_{\eta}^{3}}{f^{4}\delta
t^{2}}$+$\times$$\times$$-\frac{2}{3\times
4!}\frac{\zeta_{\eta}^{3}}{f^{4}\delta
t^{2}}$+$\times$$\times$$-\frac{1}{4!}\frac{\zeta_{\eta}^{3}}{f^{4}\delta
t^{2}}$+$\times$$\times$$-\frac{1}{4!}\frac{\zeta_{\eta}^{3}}{f^{4}\delta
t^{2}}$+$\times$$\times$$\times$$-\frac{3}{4!}\frac{\zeta_{\eta}^{3}}{f^{4}\delta
t^{2}}$+$\times$$\times$$\times$$-\frac{3}{4!}\frac{\zeta_{\eta}^{3}}{f^{4}\delta
t^{2}}$+$\times$$\times$$\times$$-\frac{3}{2\times
4!}\frac{\zeta_{\eta}^{3}}{f^{4}\delta t^{2}}$ Figure 10: Loop correction to
$q_{d}^{4}$ vertex in SK path integral (8). The sixth 2-loop diagram contains
two topologically distinct diagrams which have reflection symmetry.

Combining all of these correction we get the renormalised $\zeta_{\eta}^{r}$
as

$\zeta_{\eta}^{r}=\zeta_{\eta}-\frac{7}{2}\frac{\zeta_{\eta}^{2}}{f^{2}\delta
t}+\frac{149}{12}\frac{\zeta_{\eta}^{3}}{f^{4}(\delta t)^{2}}\,.$ (41)

### Correction to $q_{d}^{3}\dot{q}$ vertex:

The correction to $q_{d}^{3}\dot{q}$ vertex gets contribution from the
diagrams shown in figure 11.

$\times$$i\zeta_{\gamma}^{r}$$=$$\times$$i\zeta_{\gamma}$+$\times$$\times$$-\frac{i}{2}\frac{\zeta_{\gamma}\zeta_{\eta}}{f^{2}\delta
t}$+$\times$$\times$$-i\frac{\zeta_{\gamma}\zeta_{\eta}}{f^{2}\delta
t}$+$\times$$\times$$\frac{i}{2}\frac{\zeta_{\eta}\zeta_{\gamma}^{2}}{f^{4}\delta
t^{2}}$+$\times$$\times$$\frac{i}{4}\frac{\zeta_{\eta}\zeta_{\gamma}^{2}}{f^{4}\delta
t^{2}}$+$\times$$\times$$\frac{i}{3}\frac{\zeta_{\eta}\zeta_{\gamma}^{2}}{f^{4}\delta
t^{2}}$+$\times$$\times$$\frac{i}{2}\frac{\zeta_{\eta}\zeta_{\gamma}^{2}}{f^{4}\delta
t^{2}}$+$\times$$\times$$\frac{i}{2}\frac{\zeta_{\eta}\zeta_{\gamma}^{2}}{f^{4}\delta
t^{2}}$+$\times$$\times$$\times$$\frac{i}{2}\frac{\zeta_{\eta}\zeta_{\gamma}^{2}}{f^{4}\delta
t^{2}}$+$\times$$\times$$\times$$\frac{i}{2}\frac{\zeta_{\eta}\zeta_{\gamma}^{2}}{f^{4}\delta
t^{2}}$+$\times$$\times$$\times$$\frac{i}{4}\frac{\zeta_{\eta}\zeta_{\gamma}^{2}}{f^{4}\delta
t^{2}}$ Figure 11: Loop correction to $q_{d}^{3}\dot{q}$ vertex in SK path
integral (8). The sixth 2-loop diagram contains two topologically distinct
diagrams which have reflection symmetry, as done in figure 10.

Combining all of these correction we get the renormalised $\zeta_{\gamma}^{r}$
as

$\zeta_{\gamma}^{r}=\zeta_{\gamma}-\frac{3\zeta_{\eta}\zeta_{\gamma}}{2f^{2}\delta
t}+\frac{10}{3}\frac{\zeta_{\eta}^{2}\zeta_{\gamma}}{f^{4}(\delta t)^{2}}\,.$
(42)

##  Appendix B Functional derivation of Schwinger-Keldysh propagators

The propagators are obtained by considering the quadratic part of path-
integral (9). Then we introduce sources corresponding to all variables. The
propagators for a generic initial condition is hard to compute. However, we
can compute the thermal propagators. We shift the variables and finally take
functional derivative of the path integral w.r.t. the sources and obtain the
thermal propagators. The quadratic path integral with the sources is given by

$\begin{split}\mathcal{Z}_{0}[J_{\dot{q}},J_{d},J_{c},J_{\bar{c}}]&=\int\mathcal{D}\dot{q}\,\mathcal{D}q_{d}\,\mathcal{D}\bar{c}\,\mathcal{D}c\,\,\exp[-\int
dt\left\\{J_{\dot{q}}\,q_{d}+J_{d}\,\dot{q}+J_{c}\,c+\bar{c}\,J_{\bar{c}}\right\\}\Big{]}\\\
&\times\exp[-\int
dt\Bigg{\\{}iq_{d}\ddot{q}+i\gamma^{r}q_{d}\dot{q}+\frac{f^{r}}{2}q_{d}^{2}+i\bar{c}\left(\frac{d}{dt}+\gamma\right)c\Bigg{]}\,.\end{split}$
(43)

Here $J_{\dot{q}},\,J_{d},\,J_{c}$ and $J_{\bar{c}}$ are sources of
$q_{d},\,\dot{q},\,c$ and $\bar{c}$ respectively. It may seem that $q_{d}$ and
$\dot{q}$ sources are flipped. But, it is a convention followed in SK
literature.

The path integral in Fourier space

$\begin{split}\mathcal{Z}_{0}=&\int\mathcal{D}\dot{q}\,\mathcal{D}q_{d}\,\mathcal{D}\bar{c}\,\mathcal{D}c\times
e^{-\int\frac{d\omega}{2\pi}\left\\{\tilde{J}_{\dot{q}}(-\omega)\,\tilde{q}_{d}(\omega)+\tilde{J}_{d}(-\omega)\,\tilde{\dot{q}}(\omega)+\tilde{J}_{c}(-\omega)\,\tilde{c}(\omega)+\tilde{\bar{c}}(-\omega)\,\tilde{J}_{\bar{c}}(\omega)\right\\}}\\\
&\qquad\qquad\qquad\times
e^{-\int\frac{d\omega}{2\pi}\left\\{\tilde{q}_{d}(-\omega)\left(\omega+i\gamma^{r}\right)\tilde{\dot{q}}+\frac{f^{r}}{2}\tilde{q}_{d}(-\omega)\tilde{q}_{d}(\omega)+\tilde{\bar{c}}(-\omega)\left(\omega+i\gamma\right)\tilde{c}(\omega)\right\\}}\,.\end{split}$
(44)

We give a shift to the variables as the following such that the variables and
their corresponding sources are decoupled.

$\displaystyle\tilde{q}_{d}^{\prime}(\omega)$ $\displaystyle=$
$\displaystyle\tilde{q}_{d}(\omega)+\frac{\tilde{J_{d}}(\omega)}{-\omega+i\gamma^{r}}\,,$
(45) $\displaystyle\tilde{\dot{q}}^{\prime}(\omega)$ $\displaystyle=$
$\displaystyle\tilde{\dot{q}}(\omega)-\frac{f^{r}\,\tilde{J}_{d}(\omega)}{(-\omega+i\gamma^{r})(\omega+i\gamma^{r})}+\frac{\tilde{J}_{\dot{q}}(\omega)}{\omega+i\gamma^{r}}\,,$
(46)
$\displaystyle\tilde{c}^{\prime}(\omega)=\tilde{c}(\omega)+\frac{\tilde{J}_{\bar{c}}(\omega)}{\omega+i\gamma}\,,\qquad\tilde{\bar{c}}^{\prime}(\omega)=\tilde{\bar{c}}(\omega)+\frac{\tilde{J}_{c}(\omega)}{-\omega+i\gamma}\,.$
(47)

The corresponding path-integral is given by

$\begin{split}\mathcal{Z}_{0}=&\int\mathcal{D}\dot{q}^{\prime}\mathcal{D}q^{\prime}_{d}\mathcal{D}\bar{c}^{\prime}\mathcal{D}c^{\prime}\,\,e^{\int\frac{d\omega}{2\pi}\left\\{\frac{f^{r}}{2}\frac{\tilde{J}_{d}(-\omega)\,\tilde{J}_{d}(\omega)}{\omega^{2}+(\gamma^{r})^{2}}-\frac{\tilde{J}_{\dot{q}}(-\omega)\,\tilde{J}_{d}(\omega)}{\omega-i\gamma^{r}}+\frac{\tilde{J}_{c}(-\omega)\tilde{J}_{\bar{c}}(\omega)}{\omega+i\gamma}\right\\}}\\\
&\qquad\qquad\times
e^{-\int\frac{d\omega}{2\pi}\left\\{\frac{f^{r}}{2}\tilde{q}^{\prime}_{d}(-\omega)\tilde{q}^{\prime}_{d}(\omega)+(\omega+i\gamma^{r})\tilde{q}^{\prime}_{d}(-\omega)\tilde{\dot{q}}^{\prime}(\omega)+(\omega+i\gamma)\tilde{\bar{c}}^{\prime}(-\omega)\tilde{c}^{\prime}(\omega)\right\\}}\,.\end{split}$
(48)

The first line in the above expression consists of only sources and the second
line does not involve any source. Thus, integrating out the second line gives
a constant which can be absorbed into the definition of $\mathcal{Z}_{0}$.
Thus we find

$\mathcal{Z}_{0}=\exp\left[\int\frac{d\omega}{2\pi}\left\\{\frac{f^{r}}{2}\frac{\tilde{J}_{d}(-\omega)\,\tilde{J}_{d}(\omega)}{\omega^{2}+(\gamma^{r})^{2}}-\frac{\tilde{J}_{\dot{q}}(-\omega)\,\tilde{J}_{d}(\omega)}{\omega-i\gamma^{r}}+\frac{\tilde{J}_{c}(-\omega)\tilde{J}_{\bar{c}}(\omega)}{\omega+i\gamma}\right\\}\right]\,.$
(49)

Fourier transforming back to real time, we get

$\mathcal{Z}_{0}=\exp\left[\int
dt_{1}dt_{2}\int\frac{d\omega}{2\pi}e^{-i\omega(t_{2}-t_{1})}\left\\{\frac{f^{r}}{2}\frac{J_{d}(t_{2})\,J_{d}(t_{1})}{\omega^{2}+(\gamma^{r})^{2}}-\frac{J_{\dot{q}}(t_{2})\,J_{d}(t_{1})}{\omega-i\gamma^{r}}+\frac{J_{c}(t_{2})J_{\bar{c}}(t_{1})}{\omega+i\gamma}\right\\}\right]$
(50)

The propagators are obtained by taking functional derivative of
$\mathcal{Z}_{0}$ w.r.t. the sources. The SK propagators are accordingly given
by

$\displaystyle\langle\dot{q}(t_{2})\dot{q}(t_{1})\rangle$ $\displaystyle=$
$\displaystyle(-1)^{2}\frac{\partial^{2}\mathcal{Z}_{0}}{\partial{J_{q}}(t_{2})\,\partial{J_{d}}(t_{1})}\Big{|}_{J_{\dot{q}},J_{d},J_{c},J_{\bar{c}}=0}=\int\frac{d\omega}{2\pi}\,\frac{f^{r}\,e^{-i\omega(t_{2}-t_{1})}}{\omega^{2}+(\gamma^{r})^{2}}\,,$
(51) $\displaystyle\langle q_{d}(t_{2})\dot{q}(t_{1})\rangle$ $\displaystyle=$
$\displaystyle(-1)^{2}\frac{\partial^{2}\mathcal{Z}_{0}}{\partial
J_{\dot{q}}(t_{2})\,\partial
J_{d}(t_{1})}\Big{|}_{J_{\dot{q}},J_{d},J_{c},J_{\bar{c}}=0}=-\int\frac{d\omega}{2\pi}\,\frac{e^{-i\omega(t_{2}-t_{1})}}{\omega-i\gamma^{r}}\,,$
(52) $\displaystyle\langle\dot{q}(t_{2})q_{d}(t_{1})\rangle$ $\displaystyle=$
$\displaystyle(-1)^{2}\frac{\partial^{2}\mathcal{Z}_{0}}{\partial
J_{d}(t_{2})\,\partial
J_{\dot{q}}(t_{1})}\Big{|}_{J_{\dot{q}},J_{d},J_{c},J_{\bar{c}}=0}=\int\frac{d\omega}{2\pi}\,\frac{e^{-i\omega(t_{2}-t_{1})}}{\omega+i\gamma^{r}}\,,$
(53) $\displaystyle\langle q_{d}(t_{2})q_{d}(t_{1})\rangle$ $\displaystyle=$
$\displaystyle(-1)^{2}\frac{\partial^{2}\mathcal{Z}_{0}}{\partial
J_{\dot{q}}(t_{2})\,\partial
J_{\dot{q}}(t_{1})}\Big{|}_{J_{\dot{q}},J_{d},J_{c},J_{\bar{c}}=0}=0\,.$ (54)

The first three propagators are called the Keldysh, advanced and retarded
propagator respectively. Performing the $\omega$ integral we get

$\displaystyle\langle\dot{q}(t_{2})\dot{q}(t_{1})\rangle$ $\displaystyle=$
$\displaystyle\frac{f^{r}}{2\gamma^{r}}e^{-\gamma^{r}|t_{2}-t_{1}|}\,,$ (55)
$\displaystyle\langle q_{d}(t_{2})\dot{q}(t_{1})\rangle$ $\displaystyle=$
$\displaystyle-i\,\Theta(t_{1}-t_{2})e^{-\gamma^{r}(t_{1}-t_{2})}\,,$ (56)
$\displaystyle\langle\dot{q}(t_{2})q_{d}(t_{1})\rangle$ $\displaystyle=$
$\displaystyle-i\,\Theta(t_{2}-t_{1})e^{-\gamma^{r}(t_{2}-t_{1})}\,,$ (57)
$\displaystyle\langle q_{d}(t_{2})q_{d}(t_{1})\rangle$ $\displaystyle=$
$\displaystyle 0\,.$ (58)

We implement the above propagators in computation of real time correlation
function in §3.1 and in §3.3. Similarly the Grassmannian propagators are
obtained as the following.

$\displaystyle\langle\bar{c}(t_{2})\,c(t_{1})\rangle$ $\displaystyle=$
$\displaystyle(-1)^{2}\frac{\partial^{2}\mathcal{Z}_{0}}{\partial
J_{\bar{c}}(t_{2})\,\partial
J_{c}(t_{1})}\Big{|}_{J_{\dot{q}},J_{d},J_{c},J_{\bar{c}}=0}$ (59)
$\displaystyle=$
$\displaystyle-\int\frac{d\omega}{2\pi}\,\frac{e^{-i\omega(t_{2}-t_{1})}}{\omega-i\gamma}\,.$

### Sub-leading correction to $f^{r},\zeta_{\eta}^{r}$ due to $\bar{c}$$c$
loop:

The correction to $f^{r}$ is obtained from a Grassmannian loop integral given
by

$\displaystyle(f^{r})^{\prime}$ $\displaystyle=$ $\displaystyle
f^{r}+(-1)(2i\zeta_{\gamma})(-1)\int\frac{d\omega}{2\pi}\frac{1}{\omega-i\gamma}$
(60) $\displaystyle=$ $\displaystyle f^{r}-\zeta_{\gamma}\,.$

The first $(-1)$ appears due to the Grassmann loop, the factor $(2i\gamma)$ is
the Feynman rules obtained from the effective path-integral (9). In our case
$f^{r}=5.29$ and $\zeta_{\gamma}=0.0042$. So, the correction is smaller than
the second order noise loop correction, thus can be ignored.

In a similar way we obtain the correction to $\zeta_{\gamma}^{r}$, given by

$\displaystyle(\zeta_{\eta}^{r})^{\prime}$ $\displaystyle=$
$\displaystyle\zeta_{\eta}^{r}+(-1)(2i\zeta_{\gamma})^{2}\int\frac{d\omega}{2\pi}\frac{-1}{(\omega-i\gamma)(\omega+i\gamma)}$
(61) $\displaystyle=$
$\displaystyle\zeta_{\eta}^{r}-\frac{2\zeta_{\gamma}^{2}}{\gamma}\,.$

The correction is again sufficiently small. Thus we conclude that the
Grassmannian correction, for the domain of our parameters, does not
significantly contribute to correlation function.

## References

  * (1) B. Chakrabarty and S. Chaudhuri, “Out of time ordered effective dynamics of a quartic oscillator,” SciPost Phys. 7 (2019) 013, arXiv:1905.08307 [hep-th].
  * (2) A. Caldeira and A. J. Leggett, “Path integral approach to quantum brownian motion,” Physica A-statistical Mechanics and Its Applications 121 (1982) 587–616.
  * (3) B. Chakrabarty, S. Chaudhuri, and R. Loganayagam, “Out of Time Ordered Quantum Dissipation,” JHEP 07 (2019) 102, arXiv:1811.01513 [cond-mat.stat-mech].
  * (4) J. S. Schwinger, “Brownian motion of a quantum oscillator,” J. Math. Phys. 2 (1961) 407–432.
  * (5) L. V. Keldysh, “Diagram technique for nonequilibrium processes,” Zh. Eksp. Teor. Fiz. 47 (1964) 1515–1527. [Sov. Phys. JETP20,1018(1965)].
  * (6) E. Calzetta and B. L. Hu, “Closed Time Path Functional Formalism in Curved Space-Time: Application to Cosmological Back Reaction Problems,” Phys. Rev. D35 (1987) 495.
  * (7) A. Kamenev, Field Theory of Non-Equilibrium Systems.  Cambridge University Press, Cambridge, 1985.  http://www.cambridge.org/us/academic/subjects/physics/condensed-matter-physics-nanoscience-and-mesoscopic-physics/field-theory-non-equilibrium-systems?format=HB&isbn=9780521760829.
  * (8) B. Chakrabarty, J. Chakravarty, S. Chaudhuri, C. Jana, R. Loganayagam, and A. Sivakumar, “Nonlinear Langevin dynamics via holography,” JHEP 01 (2020) 165, arXiv:1906.07762 [hep-th].
  * (9) C. Jana, R. Loganayagam, and M. Rangamani, “Open quantum systems and Schwinger-Keldysh holograms,” JHEP 07 (2020) 242, arXiv:2004.02888 [hep-th].
  * (10) O. Aharony, S. S. Gubser, J. M. Maldacena, H. Ooguri, and Y. Oz, “Large N field theories, string theory and gravity,” Phys. Rept. 323 (2000) 183–386, arXiv:hep-th/9905111.
  * (11) P. Glorioso, M. Crossley, and H. Liu, “A prescription for holographic Schwinger-Keldysh contour in non-equilibrium systems,” arXiv:1812.08785 [hep-th].
  * (12) J. Berges and J. Cox, “Thermalization of quantum fields from time reversal invariant evolution equations,” Phys. Lett. B 517 (2001) 369–374, arXiv:hep-ph/0006160.
  * (13) J. M. Maldacena, “Non-Gaussian features of primordial fluctuations in single field inflationary models,” JHEP 05 (2003) 013, arXiv:astro-ph/0210603.
  * (14) T. Allen, B. Grinstein, and M. B. Wise, “Non-gaussian density perturbations in inflationary cosmologies,” Physics Letters B 197 no. 1, (1987) 66–70. https://www.sciencedirect.com/science/article/pii/0370269387903431.
  * (15) A. Gangui, F. Lucchin, S. Matarrese, and S. Mollerach, “The Three point correlation function of the cosmic microwave background in inflationary models,” Astrophys. J. 430 (1994) 447–457, arXiv:astro-ph/9312033.
  * (16) V. Acquaviva, N. Bartolo, S. Matarrese, and A. Riotto, “Second order cosmological perturbations from inflation,” Nucl. Phys. B 667 (2003) 119–148, arXiv:astro-ph/0209156.
  * (17) T. Falk, R. Rangarajan, and M. Srednicki, “Dependence of density perturbations on the coupling constant in a simple model of inflation,” Phys. Rev. D 46 (Nov, 1992) 4232–4234. https://link.aps.org/doi/10.1103/PhysRevD.46.4232.
  * (18) A. A. Starobinsky, “Stochastic de sitter (inflationary) stage in the early universe,” in Field Theory, Quantum Gravity and Strings, H. J. de Vega and N. Sánchez, eds., pp. 107–126.  Springer Berlin Heidelberg, Berlin, Heidelberg, 1986.
  * (19) W. D. Goldberger and I. Z. Rothstein, “An Effective field theory of gravity for extended objects,” Phys. Rev. D 73 (2006) 104029, arXiv:hep-th/0409156.
  * (20) C. R. Galley and M. Tiglio, “Radiation reaction and gravitational waves in the effective field theory approach,” Phys. Rev. D 79 (2009) 124027, arXiv:0903.1122 [gr-qc].
  * (21) P. C. Martin, E. D. Siggia, and H. A. Rose, “Statistical Dynamics of Classical Systems,” Phys. Rev. A 8 (1973) 423–437.
  * (22) C. De Dominicis and L. Peliti, “Field Theory Renormalization and Critical Dynamics Above t(c): Helium, Antiferromagnets and Liquid Gas Systems,” Phys. Rev. B 18 (1978) 353–376.
  * (23) H. Janssen, “On a lagrangean for classical field dynamics and renormalization group calculations of dynamical critical properties,” Zeitschrift für Physik B Condensed Matter and Quanta 23 (01, 1976) 377–380.
  * (24) E. Wang and U. W. Heinz, “A Generalized fluctuation dissipation theorem for nonlinear response functions,” Phys. Rev. D 66 (2002) 025008, arXiv:hep-th/9809016.
  * (25) J. von Neumann, “Various techniques used in connection with random digits,” in Monte Carlo Method, A. S. Householder, G. E. Forsythe, and H. H. Germond, eds., vol. 12 of National Bureau of Standards Applied Mathematics Series, ch. 13, pp. 36–38.  US Government Printing Office, Washington, DC, 1951.
  * (26) C. Maes and K. Netocny, “Time-Reversal and Entropy,” arXiv e-prints (Feb., 2002) cond–mat/0202501, arXiv:cond-mat/0202501 [cond-mat.stat-mech].
  * (27) R. J. Harris and G. M. Schütz, “Fluctuation theorems for stochastic dynamics,” Journal of Statistical Mechanics: Theory and Experiment 2007 no. 07, (Jul, 2007) P07020–P07020. https://doi.org/10.1088/1742-5468/2007/07/p07020.
  * (28) J. M. R. Parrondo, C. V. den Broeck, and R. Kawai, “Entropy production and the arrow of time,” New Journal of Physics 11 no. 7, (Jul, 2009) 073008. https://doi.org/10.1088/1367-2630/11/7/073008.
  * (29) D. Andrieux, P. Gaspard, T. Monnai, and S. Tasaki, “The fluctuation theorem for currents in open quantum systems,” New Journal of Physics 11 no. 4, (Apr, 2009) 043014. https://doi.org/10.1088/1367-2630/11/4/043014.
  * (30) E. H. Feng and G. E. Crooks, “Length of time’s arrow,” Phys. Rev. Lett. 101 (Aug, 2008) 090602. https://link.aps.org/doi/10.1103/PhysRevLett.101.090602.
  * (31) I. Oppenheim, “Nonlinear nonequilibrium thermodynamics i. linear and nonlinear fluctuation-dissipation theorems,” Journal of Statistical Physics 77 no. 5, (1994) 1109–1110. https://doi.org/10.1007/BF02183157.
  * (32) R. L. Stratonovich, “Nonlinear nonequilibrium thermodynamics ii,”  1992\.
  * (33) . Øksendal, B. K. (Bernt Karsten), Stochastic differential equations : an introduction with applications.  Sixth edition. Berlin ; New York : Springer, $[$2003$]$ ©2003, $[$2003$]$.  https://search.library.wisc.edu/catalog/999949285502121.
  * (34) S. Chaudhuri, “ unpublished notes ,”.
  * (35) L. Onsager, “Reciprocal Relations in Irreversible Processes. I.,” Phys. Rev. 37 no. 4, (1931) 405–426.
  * (36) L. Onsager, “Reciprocal Relations in Irreversible Processes. II.,” Phys. Rev. 38 no. 12, (1931) 2265–2279.
  * (37) H. B. G. Casimir, “On onsager’s principle of microscopic reversibility,” Rev. Mod. Phys. 17 (Apr, 1945) 343–350. https://link.aps.org/doi/10.1103/RevModPhys.17.343.
  * (38) F. Haehl, R. Loganayagam, and M. Rangamani, “Schwinger-keldysh formalism i: Brst symmetries and superspace,” Journal of High Energy Physics 2017 (10, 2016) .
  * (39) F. M. Haehl, R. Loganayagam, P. Narayan, and M. Rangamani, “Classification of out-of-time-order correlators,” SciPost Phys. 6 no. 1, (2019) 001, arXiv:1701.02820 [hep-th].
  * (40) S. Chaudhuri, C. Chowdhury, and R. Loganayagam, “Spectral Representation of Thermal OTO Correlators,” JHEP 02 (2019) 018, arXiv:1810.03118 [hep-th].
  * (41) S. Chaudhuri and R. Loganayagam, “Probing Out-of-Time-Order Correlators,” JHEP 07 (2019) 006, arXiv:1807.09731 [hep-th].
  * (42) R. Kubo, “Statistical mechanical theory of irreversible processes. 1. General theory and simple applications in magnetic and conduction problems,” J. Phys. Soc. Jap. 12 (1957) 570–586.
  * (43) P. C. Martin and J. S. Schwinger, “Theory of many particle systems. 1.,” Phys. Rev. 115 (1959) 1342–1373.
  * (44) S. Chib and E. Greenberg, “Understanding the metropolis-hastings algorithm,” The American Statistician 49 no. 4, (1995) 327–335.
  * (45) R. P. Feynman and F. L. Vernon, Jr., “The Theory of a general quantum system interacting with a linear dissipative system,” Annals Phys. 24 (1963) 118–173. [Annals Phys.281,547(2000)].


# AI-Empowered RIS-Assisted Networks: CV-Enabled RIS Selection and DNN-Enabled
Transmission

Conggang Hu, Yang Lu, , Hongyang Du, , Mi Yang ,
Bo Ai, , Dusit Niyato Conggang Hu and Yang Lu are with the School of Computer
Science and Technology, Beijing Jiaotong University, Beijing 100044, China
(e-mail: 23120352@bjtu.edu.cn,yanglu@bjtu.edu.cn).Hongyang Du and Dusit Niyato
are with the School of Computer Science and Engineering, Nanyang Technological
University, Singapore 639798 (e-mail:
hongyang001@e.ntu.edu.sg,dniyato@ntu.edu.sg).Mi Yang and Bo Ai are with the
School of Electronics and Information Engineering, Beijing Jiaotong
University, Beijing 100044, China (e-mail<EMAIL_ADDRESS>boai@bjtu.edu.cn).

###### Abstract

This paper investigates artificial intelligence (AI) empowered schemes for
reconfigurable intelligent surface (RIS) assisted networks from the
perspective of fast implementation. We formulate a weighted sum-rate
maximization problem for a multi-RIS-assisted network. To avoid huge channel
estimation overhead due to activate all RISs, we propose a computer vision
(CV) enabled RIS selection scheme based on a single shot multi-box detector.
To realize real-time resource allocation, a deep neural network (DNN) enabled
transmit design is developed to learn the optimal mapping from channel
information to transmit beamformers and phase shift matrix. Numerical results
illustrate that the CV module is able to select of RIS with the best
propagation condition. The well-trained DNN achieves similar sum-rate
performance to the existing alternative optimization method but with much
smaller inference time.

###### Index Terms:

AI, RIS, CV, DNN.

##  I Introduction

Providing high-quality wireless coverage for mobile intelligent applications
has become the central consideration of B5G wireless communication systems. As
a wireless coverage enabler, reconfigurable intelligent surface (RIS) draws
great attention due to its low energy consumption, easy deployment, and
flexible programming[1]. In particular, by jointly designing the transmit
beamformer at the transmitter and tuning the reflecting elements on the RIS,
the propagation condition is improved by mitigating the impact of obstacles on
the radio frequency (RF) signals. Nevertheless, the RIS complicates the
channel estimation and resource allocation, and the complexity is proportional
to the number of reflecting elements, which challenges the large-scale
deployment of RIS[2].

Recently, there has been increasing interest in leveraging artificial
intelligence (AI) to empower wireless networks [3]. On the one hand, the
computer vision (CV) methods have been employed to distinguish line-of-sight
(LoS) link and Non-LoS (NLoS) link to facilitate precise channel estimation
[4]. On the other hand, the “learning-to-optimize” paradigm shows great
advantages in realizing near-optimal and real-time computation when handling
the resource allocation problems in wireless networks [5]. Consequently,
applying AI methods to tackle the challenges in designing RIS-assisted
networks is a promising approach.

Thus far, there have been some existing works on AI-empowered RIS-assisted
networks. In [6], a CV-based approach to aid RIS for dynamic beam tracking was
proposed, which validates the effectiveness of using visual information to
control the RIS. In [7], a CV module was utilized to compute the distances
between RISs and mobile users for reconfigurable radiation pattern selection,
and the robustness of CV against the blockages was demonstrated numerically.
Nevertheless, the existing CV-assisted RIS networks activated all reflecting
elements. Most existing works show that the performance gain tends to
saturation with the number of reflecting elements [8]. Therefore, activating
partial reflecting elements not only mitigates the heavy overhead for channel
estimation but also induces limited performance loss. Besides, in [9], the
deep neural network (DNN)-enabled transmit design was proposed for RIS-
assisted networks to realize end-to-end mapping and joint beamforming of
transmitter and RIS. In [10], the DNN was employed to optimize the phase shift
coefficients with regard to bit error rate for RIS-assisted cooperative
communications. However, the supervised learning was adopted in [9, 10] which
burdens the generation of the training set. In [11], the DNN-enabled transmit
design was shown to achieve good performance via unsupervised learning.

In this paper, we consider a multi-RIS-assisted multi-user multiple-input-
multiple-output (MISO) network and formulate a weighted sum-rate maximization
problem. To reduce channel estimation overhead, a CV-enabled RIS selection
scheme is proposed to activate the RIS with the best propagation condition
among a set of RISs. Moreover, to realize real-time optimization, a DNN-
enabled transmit design trained via unsupervised learning is proposed to
jointly optimize the transmit beamformers and the phase shift matrix.
Numerical results validate the proposed schemes. It is shown that both the CV
module and the DNN-based resource allocation are able to realize fast
implementation with exceptional performance.

The rest of this paper is organized as follows. Section II introduces the
system model. The proposed CV-enabled RIS selection and DNN-enabled transmit
design are presented in Section III and Section IV, respectively. Numerical
results are presented in Sections V. Finally, Section VI concludes this paper.

##  II System Model

Figure 1: (i) A top view of a multi-RIS-assisted network. (ii) Graph
representation via CV. (iii) RIS selection.

Considered a multi-RIS-assisted multi-user MISO network as shown in Figure
1(i) where a $N_{\rm T}$-antenna transmitter serves $K$ mobile users with the
help of $2L$ passive RISs. Each RIS is with $N$ reflecting elements. The
transmitter intends to activate one111As each RIS is equipped with a large
number of passive reflecting elements, a huge overhead for channel estimation
of all transmitter-RIS links shall be induced. Therefore, we activate only one
RIS at a time. RIS of $2L$ RISs to enhance the information transmission. Due
to obstacles, the LoS links may be unavailable between the transmitter/RISs
and users, and propagation attenuation is induced on NLoS links. It is
difficult for traditional channel estimation methods to distinguish LoS link
and NLoS link in complex wireless environment [12], which however, can be
easily realized by CV. By input a top view of the considered system into a CV
module, it is capability of detecting and tracking the users and obstacles.
For clarity, we use ${\cal K}\triangleq\\{1,...,K\\}$, ${\cal
L}\triangleq\\{1,...,2L\\}$ and ${\cal N}\triangleq\\{1,...,N\\}$ to denote
the sets of users, RISs and reflecting elements of each RIS, respectively.

Denote the channel between the transmitter and the RIS $l$ by ${\bf
H}_{l}\in{\mathbb{C}}^{{N}\times{N_{\rm T}}}$, the channel between the
transmitter and the user $k$ by ${\bf g}_{k}\in{\mathbb{C}}^{N_{\rm T}}$ and
the channel between the RIS $l$ and the user $k$ by ${\bf
h}_{l,k}\in{\mathbb{C}}^{N}$. The received signal at the user $k$ when the RIS
$l$ is activated is given by

$\displaystyle{y_{l,k}}=\left({{\bf{g}}_{k}^{H}+{\bf{h}}_{l,k}^{H}{\bm{\Phi}}_{l}{{\bf{H}}}_{l}}\right)\sum\nolimits_{i\in\mathcal{K}}{{{\bf{w}}_{i}}}{s_{i}}+{n_{k}},$
(1)

where $s_{i}$ denotes the transmit signal for user $i$ with
${\mathbb{E}}\\{|s_{i}|^{2}\\}=1$ and ${\bf w}_{i}$ denotes the corresponding
beamformer, ${\bm{\Phi}}_{l}\triangleq{\rm
diag}\\{[e^{j\phi^{(l)}_{1}},e^{j\phi^{(l)}_{2}},...,e^{j\phi^{(l)}_{N}}]\\}\in\mathbb{C}^{N\times
N}$ denotes the phase shift matrix of the RIS $l$ with
$\phi^{(l)}_{n}\in[0,2\pi)$, and $n_{k}\sim\mathcal{CN}({0,{\sigma_{k}^{2}}})$
denotes the additive white Gaussian noises (AWGN) at the user $k$.

When the RIS $l$ is activated, the received information rate at the user $k$
is given by

$\displaystyle{{R_{k}}\left({\left\\{{{\bf{w}}_{i}}\right\\},{{\bm{\Phi}}_{l}}}\right)}=$
(2)
$\displaystyle{\log_{2}}\left({1+\frac{{{{\left|{\left({{\bf{g}}_{k}^{H}+{\bf{h}}_{k}^{H}{\bm{\Phi}}_{l}{\bf{H}}_{l}}\right){{\bf{w}}_{k}}}\right|}^{2}}}}{{\sum\nolimits_{i\in\mathcal{K}\backslash
k}{{{\left|{\left({{\bf{g}}_{k}^{H}+{\bf{h}}_{k}^{H}{\bm{\Phi}}_{l}{\bf{H}}_{l}}\right){{\bf{w}}_{i}}}\right|}^{2}}}+\sigma_{k}^{2}}}}\right)$

Our goal is to select the RIS with the best cascade channel condition, and
then, maximize the weighted sum rate under constraints of the power budget at
the transmitter and the phase shift coefficient of each reflecting element.
The corresponding optimization problem is mathematically expressed as

$\displaystyle{{\bf P}_{0}:}$
$\displaystyle\mathop{\max}\limits_{{\left\\{{{{\bf{w}}_{i}}}\right\\},l,{\bm{\Phi}}_{l}}}\sum\nolimits_{k\in{\cal
K}}\omega_{k}\mathop{\max}\limits_{l}\left\\{{{R_{k}}\left({\left\\{{{\bf{w}}_{i}}\right\\},{{\bm{\Phi}}_{l}}}\right)}\right\\}$
(3a) $\displaystyle{\rm s.t.}~{}$ $\displaystyle\sum\nolimits_{i\in{\cal
K}}{\left\|{{{\bf{w}}_{i}}}\right\|}_{2}^{2}\leq{P_{\rm max}},$ (3b)
$\displaystyle\phi^{\left(l\right)}_{n}\in\left[0,2\pi\right),$ (3c)
$\displaystyle\forall i,k\in{\mathcal{K}},\forall l\in\mathcal{L},\forall
n\in{\mathcal{N}},$ (3d)

where $\omega_{k}$ denotes the weight factor of the user $k$ and $P_{\rm max}$
denotes the power budget at the transmitter. To achieve fast implementation, a
CV-enabled RIS selection and DNN-enabled transmit design are respectively
given in the subsequent two sections.

##  III CV-enabled RIS Selection

A Cartesian coordinate system is defined with the location of transmitter
being fixed by $(0,0)$ and that of the RIS $l$ being fixed by $(a_{l},b_{l})$.
Then, the effective pathloss of the $l$-th transmitter-RIS link is denoted by

$e^{({\rm T-R})}_{l}\triangleq{\alpha_{l}}{\sqrt{a_{l}^{2}+b_{l}^{2}}},$

where $\alpha_{l}\in[1,\infty)$ represents the attenuation factor, where
$\alpha_{l}=1$ denotes the LoS link, and $\alpha_{l}>1$ denotes the NLoS link
and the value of $\alpha_{l}$ is determined by the impact of the obstacles on
the NLoS link.

Denote the coordinate of the user $k$ by $(x_{k},y_{k})$. Then, the effective
pathloss of the $(l,k)$-th RIS-user link is given by

$e^{({\rm
R-U})}_{l,k}\triangleq{\beta_{l,k}}\sqrt{{{({a_{l}}-{x_{k}})}^{2}}+{{({b_{l}}-{y_{k}})}^{2}}},$

where $\beta_{l,k}\in[1,\infty)$ represents the obstacle-aware attenuation
factor.

The RIS selection problem is expressed as

$\displaystyle{l^{\star}}=\mathop{\arg\min}\limits_{l\in{\cal
L}}\sum\nolimits_{k\in{\cal
K}}\underbrace{e_{l}^{({\rm{T-R}})}e_{l,k}^{({\rm{R-U}})}}_{{\rm{Effective~{}pathloss~{}of~{}cascade~{}channel}}},$
(4)

where ${l^{\star}}$ denotes the index of an RIS to be activated.

The main idea of the CV-enabled RIS selection is to utilize the visual
information via a CV module to compute the effective pathloss of cascade
channel (defined in (4)) associated with each RIS to select the RIS with the
best propagation condition. Therefore, the key step to solve the problem (4)
is to obtain $\\{(x_{k},y_{k})\\}$ and $\\{\alpha_{l},\beta_{l,k}\\}$, which
is realized via a single shot multi-box detector (SSD) illustrated in Figure 2
and described as follows.

###  III-1 SSD-based user and obstacle detection

With a top view of the considered system, the users and obstacles are detected
by a one-stage detector in order to realize fast detection.

In particular, the VGG16 is adopted as the backbone network to accept a frame
of the top view and yield its feature maps, with which a SSD [13] directly
predicts category (i.e., user or obstacle) scores of all objects required to
detect and offsets of multiple-scale default boxes. The category of one object
is determined by comparing its category score and a pre-defined threshold
(e.g., $0.5$). Besides, the offsets of default boxes can be utilized to
compute the central coordinate, height and width of the bounding box of each
object.

###  III-2 User and obstacle location

With the location information of bounding boxes of objects, more precise
coordinates of users and obstacles can be obtained by the edge detection.

In particular, the location information of bounding boxes enables each object
to be roughly cropped from the frame to reduce the background noise during the
edge detection. Then, the edge information of each object is detected via the
Canny algorithm [14] and shaped into a polygon via the polygon approximation
algorithm [15]. The coordinates of vertexes of polygon is also obtained, which
are used to calculate the coordinates of users, i.e., $(x_{k},y_{k})$, or the
coordinates of vertexes of obstacles.

With the location information of users and obstacles, the penetration distance
of each transmitter-RIS or RIS-user link over obstacles is obtained, which can
be utilized to calculate the corresponding attenuation factor, i.e.,
$\alpha_{l}$ or $\beta_{l,k}$.

Once $\\{(x_{k},y_{k})\\}$ and $\\{\alpha_{l},\beta_{l,k}\\}$ are available,
the problem (4) can be directly solved and the RIS with the best propagation
condition is selected with its index being the optimal solution to the problem
(4).

Figure 2: SSD-based user and obstacle detection in the CV module.

##  IV DNN-enabled Transmit Design

Once (4) is solved, ${l^{\star}}$ is obtained, with which the problem (3) is
rewritten as

$\displaystyle{{\bf P}_{1}:}$
$\displaystyle\mathop{\max}\limits_{{\left\\{{{{\bf{w}}_{i}}}\right\\},{\bm{\Phi}}_{{l^{\star}}}}}\sum\nolimits_{k\in{\cal
K}}\alpha_{k}{{R_{k}}\left({\left\\{{{\bf{w}}_{i}}\right\\},{{\bm{\Phi}}_{l^{\star}}}}\right)}$
(5a) $\displaystyle{\rm s.t.}~{}$
$\displaystyle\phi^{\left({l^{\star}}\right)}_{n}\in\left[0,2\pi\right),$ (5b)
$\displaystyle{\rm(\ref{p0:b})},\forall i,k\in{\mathcal{K}},\forall
n\in{\mathcal{N}}.$ (5c)

For notation simplicity, we use $\phi_{n}$ instead of
$\phi_{n}^{\left({{l^{\star}}}\right)}$ and ${\bf h}_{k}$ instead of ${\bf
h}_{{l^{\star}},k}$.

Figure 3: DNN-enabled transmit design.

Each beamformer can be represented by

$\displaystyle{\bf w}_{k}={\sqrt{p}_{k}}{{\overline{\bf w}}_{k}},$ (6)

where $p_{k}\in{\mathbb{R}}^{+}$ and ${{\overline{\bf w}}_{k}}$ satisfying
$\|{{\overline{\bf w}}_{k}}\|_{2}=1$ denote the power component and the
direction component, respectively.

Define

$\displaystyle{\bm{\phi}}$
$\displaystyle\triangleq{\left[{\phi_{1},\phi_{2},...,\phi_{N}}\right]^{T}}\in{\mathbb{R}}^{N\times
1},$ (7) $\displaystyle{\bf f}_{k}\left({\bm{\phi}}\right)$
$\displaystyle\triangleq{{{{{{\bf{g}}_{k}^{H}+{\bf{h}}_{k}^{H}{\bm{\Phi}}_{l}{\bf{H}}_{l}}}}}}\in{\mathbb{C}}^{1\times
N_{\rm T}},$ (8)

and

$\displaystyle{\bf G}\left({\bm{\phi}}\right)\triangleq\left[{\bf
f}_{1}\left({\bm{\phi}}\right);{\bf f}_{2}\left({\bm{\phi}}\right);...;{\bf
f}_{K}\left({\bm{\phi}}\right)\right]\in{\mathbb{C}}^{K\times{N_{\rm T}}}.$
(9)

With the minimum mean square error (MMSE) scheme, the direction of the $k$-th
beamformer is set by

${{\overline{\bf w}}_{k}}=\frac{{\bf v}_{k}({\bm{\phi}})}{\|{\bf
v}_{k}({\bm{\phi}})\|_{2}},$

where ${{\bf v}_{k}}({\bm{\phi}})$ is the $k$-th column of

${\bf V}({\bm{\phi}})\triangleq{\bf G}({\bm{\phi}})^{H}({\bf
G}({\bm{\phi}}){\bf G}({\bm{\phi}})^{H}+\sigma^{2}{\bf I}_{K})^{-1}.$

Then, the problem (5) is approximated by the following the problem (10), i.e.,

$\displaystyle{{\bf P}_{2}:}$
$\displaystyle\mathop{\max}\limits_{\left\\{p_{k},\phi_{n}\right\\}}\sum\nolimits_{k\in{\cal
K}}\alpha_{k}{{{\widetilde{R}}_{k}}\left(\left\\{p_{k},\phi_{n}\right\\}\right)}$
(10a) $\displaystyle{\rm s.t.}~{}$
$\displaystyle\phi_{n}\in\left[0,2\pi\right),$ (10b)
$\displaystyle\sum\nolimits_{i\in{\cal K}}p_{i}\leq{P_{\rm max}},$ (10c)
$\displaystyle\forall i,k\in{\mathcal{K}},\forall n\in{\mathcal{N}},$ (10d)

where

$\displaystyle{{\widetilde{R}}_{k}}$
$\displaystyle\left(\left\\{p_{k},\phi_{n}\right\\}\right)$ (11)
$\displaystyle={\log_{2}}\left({1+\frac{{{p_{k}}{{\left|{{{\bf{f}}_{k}}\left({\bm{\phi}}\right){{\overline{\bf
w}}_{k}}}\right|}^{2}}}}{\sum\nolimits_{i\in{{\cal K}\backslash
k}}{{p_{i}}{{\left|{{{\bf{f}}_{k}}\left({\bm{\phi}}\right){{\overline{\bf
w}}_{k}}}\right|}^{2}}+\sigma_{k}^{2}}}}\right).$

Our goal is to train a DNN to learn the mapping from $\\{{\bf H},{\bf
h}_{k},{\bf g}_{k}\\}$ to the near-optimal solution to the problem (10) as
shown in Figure 3. Particularly, a fully connected neural network is
constructed which includes one input layer, $G$ hidden layers and one output
layer.

In particular, the input layer accepts an real-value input with the size of

$\displaystyle a_{1}\triangleq 2\left(N\times N_{\rm T}+K\times\left(N+N_{\rm
T}\right)\right).$ (12)

For the $g$-th ($l\in\\{1,2,...,G\\}$) hidden layer, denote the learnable
weights by ${\bm{\Omega}}_{g}\in{\mathbb{R}}^{{a_{g+1}}\times{a_{g}}}$ and
${\bm{\xi}}_{g}\in{\mathbb{R}}^{{a_{g+1}}}$, where $a_{g}/a_{g+1}$ represents
the input/output dimension of the $g$-th hidden layer. Then, the output of the
$g$-th hidden layer is given by

$\displaystyle{\bf d}_{g+1}={\rm ReLU}\left({\bm{\Omega}}_{g}{\bf
d}_{g}+{\bm{\xi}}_{g}\right),~{}g\in\\{1,2,...,G\\},$ (13)

where ${\bf d}_{g}\in{\mathbb{R}}^{a_{g}}$/${\bf
d}_{g+1}\in{\mathbb{R}}^{a_{g+1}}$ represents the input/output of the $g$-th
hidden layer and ${\rm ReLU}(\cdot)$ represents the ReLU activation function
for nonlinear mapping. The output layer is to shape ${\bf d}_{G+1}$ with the
same dimension of $\\{p_{k},\phi_{n}\\}$, i.e., $(K+N)$, which is given by

$\displaystyle\left[{\widetilde{p}}_{1},{\widetilde{p}}_{2},...,{\widetilde{p}}_{K}\right]^{T}$
$\displaystyle={\rm Sig}\left({\bm{\Theta}}_{1}{\bf
d}_{G+1}+{\bm{\mu}}_{1}\right),$ (14)
$\displaystyle\left[{\widetilde{\phi}}_{1},{\widetilde{\phi}}_{2},...,{\widetilde{\phi}}_{N}\right]^{T}$
$\displaystyle={\rm Sig}\left({\bm{\Theta}}_{2}{\bf
d}_{G+1}+{\bm{\mu}}_{2}\right),$ (15)

where ${\bm{\Theta}}_{1}\in{\mathbb{R}}^{K\times a_{G+1}}$,
${\bm{\Theta}}_{2}\in{\mathbb{R}}^{N\times a_{G+1}}$,
${\bm{\mu}}_{1}\in{\mathbb{R}}^{K}$ and ${\bm{\mu}}_{2}\in{\mathbb{R}}^{N}$
denote the learnable weights of the output layer, and ${\rm Sig}(\cdot)$
represents the Sigmoid activation function. Further, to guarantee the output
of the DNN satisfies (10b) and (10c),
$\\{{\widetilde{p}}_{k},\widetilde{\phi}_{n}\\}$ are respectively updated222As
the problem (10) aims to maximize the sum rates, the transmitter intends to
use all power budget, and (16) guarantees that $\sum\nolimits_{i\in{\cal
K}}{\widetilde{p}}_{i}={P_{\rm max}}$.  by

$\displaystyle{\widetilde{p}}_{k}:=\frac{{{e^{{{\widetilde{p}}_{k}}}}}}{{\sum\nolimits_{i\in{\cal
K}}{{e^{{{\widetilde{p}}_{i}}}}}}}{P_{{\rm{max}}}},~{}\forall
k\in{\mathcal{K}},$ (16)
$\displaystyle\widetilde{\phi}_{n}:=2\pi\widetilde{\phi}_{n},~{}\forall
n\in{\mathcal{N}}.$ (17)

Additionally, to deal with the disappearance of the gradient due to the
Sigmoid function, the batch normalization operation is adopted to process the
output data after each hidden layer.

Note that the proposed DNN is able to generate a feasible solution to the
problem (10) and thus, we can train the DNN to be competent to solve problem
(10) via unsupervised learning. Define
${\bm{\theta}}\triangleq\\{{\bm{\Omega}}_{g},{\bm{\xi}}_{g},{\bm{\Theta}}_{g},{\bm{\Theta}}_{2},{\bm{\mu}}_{1},{\bm{\mu}}_{2}\\}$
to represent all the learnable weights of the DNN. With a batch of $M$
training samples, i.e., $\\{{\bf H}^{(m)},{\bf h}_{k}^{(m)},{\bf
g}_{k}^{(m)}\\}_{m=1}^{M}$, the corresponding outputs with ${\bm{\theta}}$ are
given by
$\\{{\widetilde{p}}_{k}^{(m)},{\widetilde{\phi}}_{n}^{(m)}|{\bm{\theta}}\\}_{m=1}^{M}$.
Then, the loss function is given by

$\displaystyle{\cal
L}_{M}\left({\bm{\theta}}\right)=-\frac{1}{M}\sum\nolimits_{k\in{\cal
K}}\alpha_{k}{{R_{k}}\left(\left\\{{\widetilde{p}}_{k}^{(m)},{\widetilde{\phi}}_{n}^{(m)}|{\bm{\theta}}\right\\}\right)}$
(18)

with which ${\bm{\theta}}$ is updated via the Adam method by

${\bm{\theta}}:={\bm{\theta}}-\frac{\omega{{{\bm{v}}_{t}}\left(\nabla{\cal
L}_{M}\left({\bm{\theta}}\right)\right)}}{\sqrt{{{\bm{s}}_{t}}\left(\nabla{\cal
L}_{M}\left({\bm{\theta}}\right)\right)+\varepsilon}},$ (19)

where $\omega$ denotes the learning rate, $\varepsilon$ is a hyperparameter,
${{{\bm{v}}_{t}}(\nabla{\cal L}_{M}({\bm{\theta}}))}$ and
${{\bm{s}}_{t}}(\nabla{\cal L}_{M}({\bm{\theta}}))$ are two functions defined
by Adam to modify gradient and learning rate, respectively.

1 Training set: $\\{{\bf H}^{(d)},{\bf h}_{k}^{(d)},{\bf g}_{k}^{(d)}\\}_{d}$;

2 Initialize $\bm{\theta}$;

3 for  _epoch $e\in[0,1,\dots,E]$_ do

4 for  _minibatch $b\in[0,1,\dots,B]$_ do

5 Sample $\\{{\bf H}^{(m)},{\bf h}_{k}^{(m)},{\bf g}_{k}^{(m)}\\}_{m=1}^{M}$
from the training set;

6 Obtain the output of DNN
$\\{{\widetilde{p}}_{k}^{(m)},{\widetilde{\phi}}_{n}^{(m)}|{\bm{\theta}}\\}_{m=1}^{M}$;

7 Calculate the loss function ${\cal L}_{M}\left({\bm{\theta}}\right)$ via
(18);

8 Update ${\bm{\theta}}$ via (19);

9

10

Return $\bm{\theta}$.

Algorithm 1 Training DNN for solving problem (10)

The training process of DNN is summarized in Algorithm 1. With the output of
DNN, i.e., $\\{{\widetilde{p}}_{k},\widetilde{\phi}_{n}\\}$, the beamformers
can be recovered via (6). Besides, the computational complexity of DNN-enabled
transmit design is on the order of
${\mathcal{O}}({\sum\nolimits_{g=1}^{G}{({a_{g+1}}{a_{g}^{2}}+{a_{g+1}}})+(K+N)(1+{a_{G+1}^{2}})})$,
where the problem dimension is relevant to $a_{1}$ defined in (12).

##  V Numerical Results

This section provides numerical results to evaluate the CV-enabled RIS
selection and DNN-enabled transmit design.

Simulation scenario: The numbers of RISs, reflecting elements of each RIS,
transmit antennas and mobile users are respectively set by $2L=6$, $N=8$,
$N_{\rm T}=8$, and $K=4$. The power budget of the transmitter is set by
$P_{\rm max}=1$, the average ratio between path loss and noise power is set to
$20$ dB.

Datasets: The CV module is trained on a dataset of top views of the considered
system. In each sample, the coordinates of users and obstacles as well as the
number of obstacles are randomly generated.

The DNN is trained on a set of $100,000$ unlabeled samples. Besides, the test
set includes $1,000$ labeled samples and each sample has two labels, i.e., the
result of $25$ iterations via the AO-based algorithm (denoted by AO 25) and
the result of $50$ iterations via the AO-based algorithm (denoted by AO 50).

Our implementation is developed using Python 3.9.18 with Pytorch 2.2.0 on a
computer with Intel(R) Core(TM) i5-9300H CPU and one NVIDIA RTX 1650 GPU (4 GB
of memory).

###  V-1 Test performance of CV module

Figure 4 illustrates the RIS selection by the CV module. It is observed that
the selected RIS suffers from the least penetrate loss and multiplicative
fading, which validates the effectiveness of the CV-enabled RIS selection.
Besides, the CV module achieves $11$ FPS for $300\times 300$ input (i.e., the
top view). Note that a CV selection result is applicable for several frames
due to the low-speed mobility of users and obstacles.

Figure 4: Illustration of the CV-enabled RIS selection.

###  V-2 Test performance of DNN

In Table I, the average sum rate over $1,000$ test samples yield by the well-
trained DNN is very close to that by the AO-based algorithm (i.e., the
performance loss is lower than $10\%$) while the inference speed of DNN is one
million faster than the AO-based algorithm which is able to realize real-time
optimization.

TABLE I: Evaluation of DNN-enabled transmit design based on average sum rate and average inference time.  Method | Sum rate [bit/s/Hz] | Inference time
---|---|---
AO 25 | $37.83$ |  $69.8$ s
AO 50 | $38.02$ |  $118.9$ s
DNN | $34.44$ |  $0.02$ ms

##  VI Conclusion

This paper applied two AI methods to design multi-RIS-assisted networks to
realize fast implementation. A CV module was designed to select the RIS with
the best RF signal propagation condition among a set of RISs to reduce the
channel estimation overhead. With the selected RIS, a DNN was trained to
jointly optimize the transmit beamformers and phase shift matrix. The result
of the CV-enabled RIS selection scheme was illustrated. The well-trained DNN
was shown to be close to the AO-based algorithm but with much faster inference
speed. These results demonstrate the advantage of applying AI to wireless
networks.

## References

  * [1] R. Liu et al, “A path to smart radio environments: An industrial viewpoint on reconfigurable intelligent surfaces,” _IEEE Wireless Commun._ , vol. 29, no. 1, pp. 202-208, Feb. 2022.
  * [2] X. Wei, D. Shen, and L. Dai, “Channel estimation for RIS assisted wireless communications-Part I: Fundamentals, solutions, and future Opportunities,” _IEEE Wireless Commun. Lett._ , vol. 25, no. 5, pp. 1398-1402, May 2021.
  * [3] Y. Lu, W. Mao, H. Du, O. A. Dobre, D. Niyato, and Z. Ding, “Semantic-aware vision-assisted integrated sensing and communication: Architecture and resource allocation,” early accessed by _IEEE Wireless Commun._ , 2024.
  * [4] Z. Hua, Y. Lu, G. Pan, K. Gao, D. B. d. Costa, and S. Chen, “Computer vision-aided mmWave UAV communication systems,” _IEEE Internet Things J._ , vol. 10, no. 14, pp. 12548-12561, July 2023.
  * [5] Y. Li, Y. Lu, R. Zhang, B. Ai, and Z. Zhong, “Deep learning for energy efficient beamforming in MU-MISO networks: A GAT-based approach,” _IEEE Wireless Commun. Lett._ , vol. 12, no. 7, pp. 1264-1268, July 2023.
  * [6] M. Ouyang, F. Gao, Y. Wang, S. Zhang, P. Li, J. Ren, “Computer vision-aided reconfigurable intelligent surface-based beam tracking: Prototyping and experimental results,” _IEEE Trans. Wireless Commun._ , vol. 22, no. 12, pp. 8681-8693, Dec. 2023.
  * [7] V. Loscrí, C. Rizza, A. Benslimane, A. M. Vegni, E. Innocenti, and R. Giuliano, “BEST-RIM: A mmWave beam steering approach Bbased on computer vision-enhanced reconfigurable intelligent metasurfaces,” _IEEE Trans. Veh. Technol._ , vol. 72, no. 6, pp. 7613-7626, June 2023.
  * [8] C. Huang, A. Zappone, G. C. Alexandropoulos, M. Debbah, and C. Yuen, “Reconfigurable intelligent surfaces for energy efficiency in wireless communication,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 8, pp. 4157-4170, Aug. 2019.
  * [9] H. Jiang, L. Dai, M. Hao, and R. MacKenzie, “End-to-end learning for RIS-aided communication systems,” _IEEE Trans. Veh. Technol._ , vol. 71, no. 6, pp. 6778-6783, June 2022.
  * [10] B. Sagir, E. Aydin, and H. Ilhan, “Deep-learning-assisted IoT-based RIS for cooperative communications,” _IEEE Internet Things J._ , vol. 10, no. 12, pp. 10471-10483, June 2023.
  * [11] Y. Li, Y. Lu, B. Ai, O. A. Dobre, Z. Ding and D. Niyato, “GNN-based beamforming for sum-rate maximization in MU-MISO networks,” early accessed by _IEEE Trans. Wireless Commun._ , 2024.
  * [12] M. Yang et al., “Machine-learning-based scenario identification using channel characteristics in intelligent vehicular communications,” _IEEE Trans. Intell. Transp. Syst._ , vol. 22, no. 7, pp. 3961-3974, July 2021.
  * [13] W. Liu et al, “SSD: Single shot multibox detector,” in _Proc. ECCV_ , Sept. 2016.
  * [14] J. Canny, “A computational approach to edge detection,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. PAMI-8, no. 6, pp. 679-698, Nov. 1986.
  * [15] M. Visvalingam and J. D. Whyatt, “The Douglas-Peucker algorithm for line simplification: Re-evaluation through visualization,” _Comput. Graphics Forum_ , no. 9, pp. 213–225, Sept. 1990.