aainstitutetext: Department of Physics, Technion, Haifa 32000, Israel

# Stochastic gravity and turbulence

Sebastian Waeber a and Amos Yarom<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We study the ensemble average of the thermal expectation value of an energy
momentum tensor in the presence of a random external metric. In a holographic
setup this quantity can be read off of the near boundary behavior of the
metric in a stochastic theory of gravity. By numerically solving the
associated Einstein equations and mapping the result to the dual boundary
theory, we find that the non relativistic energy power spectrum exhibits a
power law behavior as expected by the theory of Kolmogorov.

##  1 Introduction

Turbulence is a ubiquitous phenomenon which manifests itself whenever fluid
flow at high Reynolds number is the appropriate physical kinematic description
of the system, be it the quark gluon plasma generated in heavy ion collisions,
fast flowing rivers, or neutron star evolution. Unfortunately, while it is a
robust phenomenon, there is much we do not know regarding turbulent behavior.
The most prominent characteristic of turbulent flow is the power law behavior
of velocity field correlators. Kolmogorov’s theory of turbulence K41
correctly predicts the power law behavior of the lower moments of these
correlators. The behavior of higher moments is currently inexplicable.

At a what seems like a completely different corner of physical theories, lies
the gauge-gravity duality Maldacena:1997re ; Witten:1998qj ; Gubser:1998bc
inspired by string theory. This duality relates, in an appropriate limit,
gravitational dynamics in asymptotically anti de Sitter (AdS) space to that of
conformal field theory. Since conformal field theories can be described by
hydrodynamics, in its regime of validity, then by the duality, this
hydrodynamic behavior must be encoded in gravitational dynamics of
asymptotically AdS black branes.

If gravity, and black holes in particular, capture turbulent dynamics then
this opens a wide range of possibilities for understanding the latter in terms
of the former. In particular, universal power law behavior of moments of the
velocity field should be captured by universal black hole behavior. While it
stands to reason that the duality implies such a relation, the chief property
of black holes responsible for turbulent flow is somewhat of a mystery.
Various attempts at identifying these key features of black holes can be found
in the literature.

Of particular relevance to our current study is the work of
Bhattacharyya:2008jc  which made the connection between the constitutive
relations of fluid dynamics and gravity in the presence of a negative
cosmological constant manifest. This work allowed the authors of Adams:2013vsa
(see also Eling:2010vr ; Carrasco:2012nf ) to obtain a conjectural relation
between the power law behavior of the energy power spectrum, predicted by
Kolmogorov, and a certain fractal dimension of the apparent horizon. The
authors then demonstrated that this relation holds in the case where the
apparent horizon was initially excited and then left to decay—corresponding to
decaying turbulent flow. In recent years there has been a significant body of
work scrutinizing, and building on, the results of Adams:2013vsa , see, e.g.,
Green:2013zba ; Balasubramanian:2013yqa ; Eling:2013sna ; Eling:2015mxa ;
Westernacher-Schneider:2015gfa ; Westernacher-Schneider:2017snn ;
Rozali:2017bll ; Westernacher-Schneider:2017xie ; Eyink:2017zfz ; Chen:2018nbh
; Andrade:2019rpn . Arguably, one of the main shortcomings of Adams:2013vsa
is that it deals with decaying turbulence where the Kolmogorov scaling
relations are notoriously difficult to capture for a prolonged time. Indeed,
steady state turbulent flow is usually achieved by introducing a random
driving force which continuously pumps energy into the fluid. In the context
of the gauge-gravity duality such a driving force should be reflected by a
random source for the bulk fields. Put differently, a random driving force is
dual to random boundary conditions for the dynamical fields, including the
metric itself. While some of the aforementioned works (Balasubramanian:2013yqa
; Westernacher-Schneider:2015gfa ; Westernacher-Schneider:2017snn ;
Westernacher-Schneider:2017xie ; Chen:2018nbh ; Andrade:2019rpn  in
particular, and also Bhattacharyya:2008ji ) have dealt with driving forces,
there is no work, to our knowledge, which properly deals with a dual
gravitational setup with random boundary conditions for the metric. The
current work aims to fill in this gap.

We have found that the main challenge in incorporating a random boundary
metric in an asymptotically AdS geometry is to obtain a well defined set of
stochastic differential equations. The crux being the non linear nature of the
Einstein equations. Not every type of random source will yield well defined
equations. In particular, as we argue in section 2, a boundary metric
characterized by white noise is incompatible with the gravitational dynamics
at hand. Instead of white noise, we have found that the metric must be
characterized by an autocorrelation function which is not a distribution and
is differentiable at least up to an appropriate, finite, order in derivatives.

After incorporating a stochastic boundary metric in a holographic setup, we
solved the stochastic Einstein equations numerically using a modest grid,
sufficient to exhibit Kolmogorov scaling behavior of the dual (non
relativistic) energy power spectrum across a decade of momenta (see figure 4).
We found, among other things, that the near AdS boundary behavior of the
metric (which determines the behavior of the stress tensor in the dual theory)
is enhanced in the deep interior due to the AdS warp factor. Thus, the horizon
provides an amplified picture of the dynamics of the fluid—inline with
qualitative notions of holography. We discuss our results in section 6.

##  2 The metric as a stochastic driving force

We often use the generating function

$Z[g_{\mu\nu}]=\int D\phi e^{-S[\phi;\,g_{\mu\nu}]}$ (1)

to characterize correlation functions of the energy momentum tensor in a
theory with some dynamical variables $\phi$ and an external metric
$g_{\alpha\beta}$, e.g.,

$\langle 0|T^{\mu\nu}|0\rangle{}_{g}=\frac{2i}{\sqrt{g}}\frac{\delta}{\delta
g_{\mu\nu}}Z[g_{\mu\nu}]\,.$ (2)

Here, and in the rest of this work, $\mu,\nu=0,\ldots d-1$ will denote the
spacetime directions of the field theory. The ensemble averaged vacuum
expectation value of a stress tensor in a quantum field theory whose metric is
characterized by a random variable, $Q_{\alpha\beta}$,

$g_{\alpha\beta}=\eta_{\alpha\beta}+Q_{\alpha\beta}\,,$ (3)

is given by

$\overline{\langle 0|T^{\mu\nu}|0\rangle}=\int
DQ_{\alpha\beta}P(Q_{\alpha\beta})\langle 0|T^{\mu\nu}|0\rangle{}_{\eta+Q}\,,$
(4)

where $P(Q)$ specifies the distribution from which $Q_{\alpha\beta}$ is drawn.
We will refer to $Q_{\alpha\beta}$ as the random, or stochastic, component of
the metric. (See Aharony:2018mjm  for a study of renormalization group flow in
a related setup.)

In this work we will be interested in thermal expectation values and not
vacuum expectation values of the energy momentum tensor in the presence of the
metric (3). However, before proceeding, it is useful to pause and consider
more broadly the properties of $T^{\mu\nu}$ in the presence of (3). Coordinate
invariance implies the Ward identity

$\nabla_{\mu}T^{\mu\nu}=0\,.$ (5)

The stochastic component of the metric, $Q_{\alpha\beta}$, contributes to (5)
through the Christoffel connection and through a possible explicit dependence
of the energy momentum tensor on the external metric (3). Thus, at the very
least, (5) will be a stochastic differential equation (SDE) where the external
random variable ($Q_{\alpha\beta}$ in our case) is differentiated and may also
appear non linearly. We refer the reader unfamiliar with the basics of
stochastic differential equations to appendix A, or to the extensive
literature on the subject specified therein.

Stochastic differential equations of the form (5) are not well defined if the
stochastic variable is characterized as white noise. Indeed, suppose
$\xi(t,\vec{x})$ denotes white noise with zero average, i.e.,
$\overline{\xi(t,\vec{x})}=0$ and
$\overline{\xi(t,\vec{x})\xi(t^{\prime},\vec{x}^{\prime})}=D(x-x^{\prime})\delta(t-t^{\prime})$
where $D$ is a smooth function of its arguments and a bar denotes an ensemble
average. It follows that rational functions of $\xi(t,\,\vec{x})$ are ill
defined, as are functions of its temporal derivatives.

To make sense of (5) in the presence of (3), we use a random force which is
determined via an Ornstein-Uhlenbeck process, or, in more involved cases, via
what we refer to as a cascade of Ornstein-Uhlenbeck processes. Recall that, in
a $0+1$ dimensional setting, an Ornstein-Uhlenbeck process, $Q(t)$, is defined
by

$\dot{Q}(t)=-\frac{Q(t)}{\tau_{1}}+\frac{\xi(t)}{\tau_{2}}$ (6)

where a dot denotes a time derivative, $\tau_{1}$ and $\tau_{2}$ are positive
numbers, and $\xi$ denotes white noise with zero average. In what follows we
use, for simplicity, $\tau_{1}=\tau_{2}=\tau$. We define a cascade of
Ornstein-Uhlenbeck processes of order $n$, $Q(t)$, $q_{1}(t)$, $\ldots$,
$q_{n}(t)$ via

$\dot{Q}(t)=-\frac{Q(t)}{\tau}+\frac{q_{1}(t)}{\tau}\,,\quad\ldots\,,\quad\dot{q}_{i}(\tau)=-\frac{q_{i}(t)}{\tau}+\frac{q_{i-1}(t)}{\tau}\,,\quad\ldots\,,\quad\dot{q}_{n}(\tau)=-\frac{q_{n}(t)}{\tau}+\frac{\xi}{\tau}\,,$
(7)

where, as before, $\xi$ denotes white noise with zero mean. (So that we may
refer to an Ornstein-Uhlenbeck process as a cascade of order $0$.) Using (7)
it is straghtforward to argue that an SDE of the form

$\dot{X}(t)=G\left(X(t),\,Q(t),\,\dot{Q}(t),\,\ddot{Q}(t),\,\ldots,\left(\frac{\partial}{\partial
t}\right)^{n}Q(t)\right)$ (8)

can be rewritten as a set of (well defined) coupled stochastic differential
equations linear in the white noise $\xi(t)$ (see e.g., (62)) as long as $Q$,
$q_{1}$, $\ldots$, $q_{n}$, are given by a cascade of at least order $n$. If
the $n$’th derivative in (8) appears linearly then a cascade of order $n-1$
will suffice.

In order to completely specify the problem it remains to state the integration
scheme associated with the stochastic noise term: as explained in appendix A,
apart from initial conditions, one needs to specify an integration scheme for
an SDE in order to obtain a unique solution to it. Which integration scheme to
use depends on the problem at hand. Often the Itô integration scheme is used
when studying formal aspects of SDE’s while the Stratonovich scheme is more
common in physical problems. In the current case, the stochastic metric is
external so that the associated integration scheme may be determined at our
discretion. After some experimentation we have decided on using the
Stratonovich scheme. Using an Itô scheme often resulted in what seemed like
unphysical features. In particular, we observed an overall decrease in the
energy density over time.

While equations (7) and (8) have been specified in a $0+1$ dimensional
setting, they can be easily generalized to include a spatial dependence. In
the remainder of this work we will focus on a $2+1$ dimensional spacetime,
where the Cartesian spatial coordinates have period $L$, $x^{i}\sim x^{i}+L$.
This implies, among other things, that momentum is discrete,
$\vec{k}_{i}=\frac{2\pi\vec{n}_{i}}{L}$ with $\vec{n}_{i}$ a pair of integers.
We choose a random metric

$Q_{\alpha,\beta}(t,\vec{x})=\left(e^{2Q(t,\vec{x})}-1\right)\sum_{i=1,2}\delta_{\alpha\,x^{i}}\delta_{\beta\,x^{i}}$
(9a) such that
$e^{2Q(t,\vec{x})}-1=q(t,\vec{x})+3\sqrt{\frac{D}{2\tau}N_{A}\left(1-e^{-\frac{2t}{\tau}}\right)}\,.$
(9b) The first term on the right hand side of (9b), $q$, is determined by a
cascade of order $0$,
$\dot{q}(t,\vec{x})=-\frac{q(t,\vec{x})}{\tau}+\frac{\xi(t,\vec{x})}{\tau}\,,$
(9c) with $\xi(t,\vec{x})$ a random variable satisfying
$\overline{\xi(t,\vec{x})}=0$ and
$\overline{\xi(t,\vec{x})\xi(t^{\prime},\vec{x}^{\prime})}=D\delta(t-t^{\prime})\sum_{{i,\,\vec{k}_{i}\in
A}}\cos\left(\vec{k}_{i}\cdot\left(\vec{x}-\vec{x}^{\prime}\right)\right)\,,$
(9d)

the sum running over all momenta $\vec{k}_{i}$ inside an annulus $A$ of inner
radius $k_{f}-\Delta k$ and outer radius $k_{f}+\Delta k$. The parameter $D$
which characterizes the strength of the random force is positive. In the
second term on the right hand side of (9b) we have defined

$N_{A}=\sum_{i,\,\vec{k}_{i}\in A}1\,,$ (10)

the ‘size’ of the annulus.

In place of (9b) and (9c), we could have chosen $Q$ itself to be specified by
a cascade of order $0$ which would have resulted in a skewed distribution
around $e^{Q}=1$. Instead, we have chosen to parameterize $e^{Q}$ as in (9b)
where the second term on the right ensures that the left hand side is
(statistically speaking) positive;

$\overline{q(t,\vec{x})q(t,\vec{x})}=\frac{D}{2\tau}N_{A}\left(1-e^{-\frac{2t}{\tau}}\right)\,.$
(11)

Correlations of the form (9d) can be obtained from the distribution

$\xi(t,\vec{x})=\sqrt{D}\sum_{i,\,\vec{k}_{i}\in
A}\left(P^{(1)}_{i}(t)\cos\left(\vec{k}_{i}\cdot\vec{x}\right)+P^{(2)}_{i}(t)\sin\left(\vec{k}_{i}\cdot\vec{x}\right)\right)$
(12)

where the $P^{(a)}_{i}(t)$ are independent random variables. At each time $t$
and for each $a=1,2$ and $\vec{k}_{i}\in A$, $P^{(a)}_{i}(t)$ is drawn from a
normal distribution, $N(0,1)$. In principle, we could have omitted the sum
over the annulus in (12) and obtained the same ensemble average as in (9d).
The sum in (12) ensures that each element of the ensemble will have an
isotropic driving force. Note that spatial derivatives of $\xi$ are well
defined—its stochastic component enters only through the $P_{i}^{(a)}$’s.

The metric specified by (3) and (9) allows us to treat conservation equations
which are at most linear in time derivatives of the random metric but may
contain arbitrary powers of it or its spatial derivatives. As we will see
shortly, (9a) has been chosen such that the resulting equations of motion are
precisely of this form.

##  3 A holographic realization

Our main goal in this work is to compute the average value of the thermal
expectation value of the energy momentum tensor, in the presence of a
stochastic boundary metric given by (3) and (9),

$\mathcal{T}^{\mu\nu}=\overline{\hbox{Tr}\left(\frac{e^{-\beta
H}}{\hbox{Tr}(e^{-\beta H})}T^{\mu\nu}\right)}\,.$ (13)

We will do so using the holographic duality Maldacena:1997re  which relates
the dynamics of field theories to the dynamics of asymptotically anti de
Sitter (AdS) spaces. In particular, suppose

$ds^{2}=\frac{\ell^{2}}{\zeta^{2}}\left(\sum_{k=0}^{\infty}g_{\mu\nu}^{(k)}(\vec{x})\zeta^{k}dx^{\mu}dx^{\nu}+d\zeta^{2}\right)$
(14)

(with latin indices $m,n=0,\ldots 3$ (in contrast to $\mu,\nu=0,\ldots,2$)
which do not run over the ‘$\zeta$’ coordinate) is a solution to the equations
of motion associated with the Einstein-Hilbert action

$S=\frac{1}{16\pi G_{N}}\int
d^{4}x\sqrt{g}\left(R+\frac{6}{\ell^{2}}\right)\,.$ (15)

This solution specifies a state, $\varrho$, in the dual field theory, such
that the expectation value of the energy momentum tensor in this state is
given by

$\hbox{Tr}\left(\varrho{{T}_{\mu\nu}}\right)\Bigg{|}_{g_{\mu\nu}=g^{(0)}_{\mu\nu}}=\frac{3\ell^{2}}{16\pi
G_{N}}g^{(3)}_{\mu\nu}\,,$ (16)

(so that the field theory is placed on a manifold with metric
$g^{(0)}_{\mu\nu}$) where $\ell^{2}/G_{N}$ is often associated with a large
number tied to the gauge group of the dual field theory. See, e..g,
deHaro:2000vlm  for details.

For instance, the black brane geometry

$ds^{2}=\frac{\ell^{2}}{\rho^{2}}\left(-\left(1-\left(\frac{\rho}{\rho_{0}}\right)^{3}\right)dt^{2}+dx_{1}^{2}+dx_{2}^{2}-2dtd\rho\right)$
(17)

is a solution to the Einstein equations following from (15) with Hawking
temperature

$T_{0}=\frac{3}{4\pi\rho_{0}}\,.$ (18)

Bringing (17) to the form (14) leads, via (16), to

$\hbox{Tr}\left(\frac{e^{-\beta H}}{\hbox{Tr}(e^{-\beta
H})}T^{\mu\nu}\right)=\frac{3\ell^{2}}{16\pi
G_{N}}\frac{64\pi^{3}T_{0}^{3}}{81}\begin{pmatrix}2&0&0\\\ 0&1&0\\\
0&0&1\end{pmatrix}\,,$ (19)

the expected expression for the energy momentum tensor of a conformal field
theory, in flat space, in a thermal state.

We propose that a holographic expression for $\mathcal{T}^{\mu\nu}$ defined in
(13) is given by

$\mathcal{T}^{\mu\nu}=\frac{3\ell^{2}}{16\pi G_{N}}\overline{g^{(3)\mu\nu}}$
(20)

where $g^{(3)\mu\nu}=g^{(0)\mu\alpha}g^{(0)\nu\beta}g^{(3)}_{\alpha\beta}$ and
$g^{(0)\mu\nu}$ is the inverse of $g^{(0)}_{\mu\nu}$. The averaging in (20) is
carried out over boundary metrics $g_{\mu\nu}^{(0)}$ specified in (3) which in
the current notation reads

$g^{(0)}_{\mu\nu}dx^{\mu}dx^{\nu}=-dt^{2}+e^{2Q(t,\vec{x})}(dx_{1}^{2}+dx_{2}^{2})\,.$
(21)

We choose a $Q$ such that $Q=0$ for $t<0$ and is a stochastic variable,
defined by (9), at later times. Put differently, in order to obtain
$\mathcal{T}^{\mu\nu}$ we must solve the Einstein equations derived from (15)
with boundary conditions given by (21) while treating the boundary value of
the metric, and therefore the metric itself, as a random variable at times
$t>0$.

To set up the problem it is convenient to use the infalling coordinate system
suggested by Chesler:2013lia  such that

$\frac{1}{\ell^{2}}ds^{2}=\Sigma(t,\vec{x},\rho)^{2}\hat{g}_{ij}(t,\vec{x},\rho)dx^{i}dx^{j}-2dt\left(F_{i}(t,\vec{x},\rho)dx^{i}+A(t,\vec{x},\rho)dt+\frac{d\rho}{\rho^{2}}\right)$
(22)

where $|\hat{g}|=1$, $t=x^{0}$ is the time coordinate, and $i=1,2$ are the
spatial coordinates. The radial coordinate $\rho$ vanishes at the asymptotic
AdS boundary. The overall factor of $\ell$ in (22) has been conveniently
chosen so that it drops out of the equations of motion:

$\displaystyle\left(\partial_{\rho}^{2}+\frac{2}{\rho}\partial_{\rho}+\frac{1}{8}\,\hat{g}^{ij}\,(\partial_{\rho}\hat{g}_{jk})\,\hat{g}^{kl}\,(\partial_{\rho}\hat{g}_{li})\right)\Sigma$
$\displaystyle=0$ (23a)
$\displaystyle\left(\delta_{i}^{j}\partial_{\rho}^{2}+\delta_{i}^{j}\frac{2}{\rho}\partial_{\rho}+P_{F}{}_{i}^{j}[\hat{g},\,\Sigma]\right)F_{j}$
$\displaystyle=S_{F\,i}[\hat{g},\,\Sigma]$ (23b)
$\displaystyle\left(\partial_{\rho}+\frac{\partial_{\rho}\Sigma}{\Sigma}\right)d_{0}\Sigma$
$\displaystyle=S_{d_{0}\Sigma}[\hat{g},\,\Sigma,\,F]$ (23c)
$\displaystyle\left(\delta^{k}_{(i}\delta^{l}_{j)}\partial_{\rho}+\frac{\partial_{\rho}\Sigma}{\Sigma}\delta^{k}_{(i}\delta^{l}_{j)}-\hat{g}^{mk}(\partial_{\rho}\hat{g}_{m(i}\delta^{l}_{j)})\right)d_{0}\hat{g}_{kl}$
$\displaystyle=S_{d_{0}\hat{g}\,}{}_{ij}[\hat{g},\,\Sigma,\,F,\,d_{0}\Sigma]$
(23d)
$\displaystyle\left(\partial_{\rho}^{2}+\frac{2}{\rho}\partial_{\rho}\right)A$
$\displaystyle=S_{A}[\hat{g},\,\Sigma,\,F,\,d_{0}\Sigma,\,d_{0}\hat{g}]$ (23e)
$\displaystyle\left(\delta^{j}_{i}\partial_{\rho}-\frac{\partial_{\rho}\Sigma}{\Sigma}\delta^{j}{}_{i}-\hat{g}^{jk}\partial_{\rho}\hat{g}_{ki}\right)d_{0}F_{j}$
$\displaystyle=S_{d_{0}F}{}_{i}[\hat{g},\,\Sigma,\,F,\,d_{0}\Sigma,\,d_{0}\hat{g},A]$
(23f) $\displaystyle d_{0}d_{0}\Sigma$
$\displaystyle=S_{d_{0}d_{0}\Sigma}[\hat{g},\,\Sigma,\,F,\,d_{0}\Sigma,\,d_{0}\hat{g},\,A,\,d_{0}F]$
(23g)

where we have defined

$d_{0}=\partial_{t}-\rho^{2}A\partial_{\rho}$ (24)

and

$\displaystyle\begin{split}P_{F}{}^{j}_{i}[\hat{g},\,\Sigma]=&-\hat{g}^{jk}\partial_{\rho}\hat{g}_{ki}+\left(-\left(\partial_{\rho^{2}}\hat{g}_{ik}+\frac{2}{\rho}\partial_{\rho}\hat{g}_{ik}\right)\hat{g}^{kj}+\partial_{\rho}\hat{g}_{ik}\hat{g}^{kl}\partial_{\rho}\hat{g}_{lm}\hat{g}^{lj}-\frac{2\partial_{\rho}\Sigma}{\Sigma}\partial_{\rho}\hat{g}_{ik}\hat{g}^{kj}\right)\\\
&+\delta_{i}^{j}\left(\frac{1}{4}\,\hat{g}^{ml}\,(\partial_{\rho}\hat{g}_{lk})\,\hat{g}^{kp}\,(\partial_{\rho}\hat{g}_{pm})-\frac{2\left(\partial_{\rho}\Sigma\right)^{2}}{\Sigma^{2}}\right)\,.\end{split}$
(25)

Expressions for the inhomogenous terms appearing on the right hand side of
(23) can be found in Chesler:2013lia . The boundary conditions for the metric
(21) imply

$\displaystyle\begin{split}\Sigma=\frac{e^{Q(t,\vec{x})}}{\rho}+\mathcal{O}(\rho^{0})\,,\quad
A=\frac{1}{2\rho^{2}}+\mathcal{O}(\rho^{-1})\,,\quad\hat{g}_{ij}=\delta_{ij}+\mathcal{O}(\rho)\,,\quad
F=\mathcal{O}(\rho)\,,\end{split}$ (26)

and we also require that the metric is regular across the apparent horizon.
Our initial conditions are such that the metric takes the form (17) at $t<0$
which ensures that the dual field theory is in thermal equilibrium at early
times.

To relate the solution to (23) to the average value of the stress tensor (20),
consider the near boundary asymptotic expansion of the dynamical fields under
the equations of motion,

$\displaystyle\begin{split}\Sigma=&\frac{e^{Q}}{\rho}+e^{Q}\lambda(t)+\mathcal{O}(\rho^{3})\\\
\hat{g}_{ij}=&\delta_{ij}+\gamma_{ij}\rho^{3}+\mathcal{O}(\rho^{4})\\\
F_{i}=&f_{i}(t,\,\vec{x})\rho+\mathcal{O}(\rho^{2})\\\
A=&\frac{1}{2\rho^{2}}+\frac{\lambda(t)-\partial_{t}{Q}}{\rho}+\left(\frac{\lambda^{2}(t)}{2}-\partial_{t}{\lambda}(t)-\lambda(t)\partial_{t}Q+\frac{1}{2}e^{-2Q}\delta^{ij}\partial_{i}\partial_{j}Q\right)+a(t,\,\vec{x})\rho\,.\end{split}$
(27)

The undetermined parameter $\lambda(t)$ is a reflection of a residual
coordinate transformation

$\rho\to\bar{\rho}(\rho)=\frac{\rho}{1+\rho\lambda(t)}\,,$ (28)

which preserves both the ansatz (22) and the boundary conditions (26). In
Chesler:2013lia  an extended version of this gauge freedom was used to tune
the parameteric location of the apparent horizon to a convenient location. In
what follows we will set $\lambda=0$.

The functions $a$, $f_{i}$ and $\gamma_{ij}$ are determined by the near
horizon boundary conditions. Together with $Q$, they specify the expectation
value of the stress tensor. Using (16) we find

$\displaystyle\begin{split}\frac{16\pi G_{N}}{3L^{2}}\hbox{Tr}\left(\varrho
T^{\mu\nu}\right)=&e^{-2Q}\begin{pmatrix}-\frac{4}{3}e^{2Q}a&f_{1}&f_{2}\\\
f_{1}&-\gamma_{22}-\frac{2}{3}a&\gamma_{12}\\\
f_{2}&\gamma_{12}&\gamma_{22}-\frac{2}{3}a\end{pmatrix}\\\
&-\frac{2}{3}e^{-2Q}\delta^{0(\mu}\delta^{\nu)i}\partial_{i}\left(e^{2Q}\delta^{jk}\partial_{j}\partial_{k}Q\right)+\frac{2}{3}e^{-2Q}\delta^{i(\mu}\delta^{\nu)j}\partial_{i}\left(e^{-2Q}\partial_{t}\partial_{j}Q\right)\\\
&+\frac{1}{3}e^{-2Q}\delta^{i\mu}\delta^{\nu
j}\delta_{ij}\delta^{kl}\partial_{k}\left(e^{-2Q}\partial_{t}\partial_{l}Q\right)\,.\end{split}$
(29)

Thus, in order to compute the ensemble averaged energy momentum tensor
$\mathcal{T}^{\mu\nu}$ in (20), it remains to compute the coefficients $a$,
$f_{i}$ and $\gamma_{ij}$ in (27). In order to do so we must resort to
numerics.

##  4 Setting up the numerical problem

In order to numerically solve (23) and the associated initial and boundary
conditions described in (17) and (26), we use the residual scaling symmetry
$\rho\to\alpha\rho$ and $x^{\mu}\to\alpha x^{\mu}$ to set $\rho_{0}=1$ in
(18). Thus, all dimensionfull parameters appearing in this section are
measured in units of the initial temperature. When we will present our results
in the next section we will go back to dimensionfull quantities by multiplying
physical quantities by appropriate powers of $\rho_{0}$.

The algorithm for solving (23) has been laid out in detail in Chesler:2013lia
: Given $\hat{g}_{ij}(t_{0},\,\vec{x})$, $a(t_{0},\vec{x})$ and
$f(t_{0},\vec{x})$ at some time $t_{0}$, we may solve the set of nested linear
differential equations (23) as follows. The second order (in $\rho$)
differential equation for $\Sigma$, (23a) can be solved for by imposing the
boundary conditions in (27) with $\lambda=0$. Once $\Sigma$ is known, we can
solve (23b) for $F_{i}$ using the asymptotic behavior in (27). With $F_{i}$
and $\Sigma$ we can solve the first order linear equation for $d_{0}\Sigma$,
(23c), as an independent variable. Using (27) we find that $d_{0}\Sigma$ is
uniquely determined (from the value of $a$). We next solve (23d) for
$d_{0}\hat{g}_{ij}$ with the asymptotic behavior given in (27) (which
translates into a vanishing $\mathcal{O}(\rho)$ term in the near boundary
expansion of (23d)). Finally, we can solve (23e) for $A$, once again with the
asymptotic behavior provided by (27) with $\lambda=0$.

At this point we can use the known value of $A$ to determine
$\partial_{t}\hat{g}_{ij}$ from $d_{0}\hat{g}_{ij}$ and then solve the
equations for $d_{0}F_{i}$ and $d_{0}d_{0}\Sigma$ to determine $\partial_{t}a$
and $\partial_{t}f_{i}$ which will allow us to step forward in time. In
practice, once the equations for $\Sigma$, $F_{i}$, $d_{0}\hat{g}$ and $A$ are
satisfied, the equations of motion for $d_{0}F_{i}$ and $d_{0}d_{0}\Sigma$
reduce to energy and momentum conservation of the boundary theory stress
tensor given by (29), $\nabla_{\mu}T^{\mu\nu}=0$. Thus, we may determine the
time derivatives of $a$ and $f_{i}$ via energy-momentum conservation.

To implement the above evolution scheme numerically, we discretize the
boundary coordinates $x^{i}$ on a periodic grid of physical size $L$ with
$N_{\parallel}$ grid points, and expand functions of $x^{i}$ in Fourier
cardinal functions. In particular the expansion coefficients of the random
force $\xi(t,\vec{x})$ are treated as stochastic functions in time, with
appropriate autocorrelation functions. Spatial derivatives of $\xi$ can be
computed from derivatives of the Fourier cardinal functions.

The radial coordinate, $\rho$, is taken to run in the range $0<\rho<2$. The
metric at early times has an apparent horizon at $\rho=1$ (with our choice of
units) and we have checked that the horizon is inside this coordinate range at
later times. We have discretized the radial coordinate on a grid with
$N_{\bot}$ grid points, and expanded functions of $\rho$ in Chebyshev cardinal
functions. Lastly, to avoid the inherently singular nature of the near
boundary behavior of some of the metric components, or of the equations of
motion themselves, we have found it convenient to switch from $\Sigma$,
$F_{i}$, $d_{0}\Sigma$, $d_{0}\hat{g}_{ij}$ and $A$ to $\widetilde{\Sigma}$,
$\widetilde{F}_{i}$, $\widetilde{d_{0}\Sigma}$, $\widetilde{d_{0}g_{ij}}$ and
$\widetilde{A}$ defined via

$\displaystyle\begin{split}\Sigma=&\frac{e^{Q}}{\rho}+\rho^{3}\widetilde{\Sigma}\\\
F_{i}=&\widetilde{F}_{i}\\\
d_{0}\Sigma=&\frac{e^{Q}}{2\rho^{2}}+\frac{1}{2}e^{-Q}\delta^{ij}\partial_{i}\partial_{j}Q+\rho\widetilde{d_{0}\Sigma}\\\
d_{0}\hat{g}_{ij}=&\rho\widetilde{d_{0}g_{ij}}\\\
A=&\frac{1}{2\rho^{2}}+\frac{-\partial_{t}{Q}}{\rho}+\frac{1}{2}e^{-2Q}\delta^{ij}\partial_{i}\partial_{j}Q+\widetilde{A}\rho\,.\end{split}$
(30)

To summarize this section let us briefly outline our evolution scheme once
more. At each time step, $t$, we start with $\hat{g}_{ij}$, $a$ and $f_{i}$.
We then determine the boundary value of the metric, characterized by
$e^{Q(t,\vec{x})}$ by solving (9b) and (9c) and drawing $P_{i}(t)$ in (12)
from a normal distribution.111In practice we have drawn $M$ values of
$\vec{k}$ from the annular region at each time step, where $M$ is roughly half
the number of grid points in $A$. The same applies to (9b).  We then solve the
nested set of linear differential equations (23a) to (23e) in terms of the
tilde’d variables in (30). To move forward in time we have used an Euler-
Maruyama algorithm (described in appendix A) using a time interval, $\Delta
t$, and, as mentioned before, treated the noise using the Stratonovich scheme.
At the end of the day, this algorithm provides us with
$\partial_{t}\hat{g}_{ij}$ from which we can obtain $\hat{g}_{ij}$ at the next
time step. Using energy momentum conservation we can then determine the time
derivatives of $a$ and $f_{i}$ and use them to obtain their value at the next
time step as well. We then use a $2/3$ filter to the updated values of
$\hat{g}_{ij}$, $a$ and $f_{i}$ and repeat the procedure. With $a$, $f_{i}$,
$\gamma_{ij}$ and $Q$ at hand, we can determine the energy momentum tensor at
all times via (29). An approximation to the ensemble averaged stress tensor
$\mathcal{T}^{\mu\nu}$ is obtained by repeating this algorithm a large,
$N_{0}$, number of times and averaging over the resulting stress tensors
obtained in each instance.

##  5 Results

The expectation value of the stress tensor obtained via the algorithm outlined
in 4 may depend, a priori, on the parameters of the noise, $\tau$, $k_{f}$,
$\Delta k$ and $D$ in (12) as well as on the length of the spatial torus $L$.
In the current text we have focused on

$\frac{\tau}{\rho_{0}}=\frac{1}{4}\,,\quad\frac{Lk_{f}}{2\pi}=21\,,\quad\frac{L\Delta
k}{2\pi}=\frac{3}{4}\,,\quad\frac{D}{\rho_{0}}=\frac{3}{1250000}\,,\quad\frac{L}{\rho_{0}}=500\,.$
(31)

We have chosen (31) since, after some exploration, these, as well as their
nearby values, allow for a stress tensor whose associated energy power
spectrum has the Kolmogorov scaling law that we’re interested in. A full study
of the dependence of the expectation value of the stress tensor on these
parameters (and other more general noise terms) is left for future work.

In order to implement the aforementioned algorithm numerically we must also
choose an appropriate spatial and temporal grid. We have found that the
somewhat modest

$N_{\bot}=70,\quad N_{\parallel}=11,\quad\frac{\Delta
t}{\rho_{0}}=\frac{1}{40}$ (32)

where sufficient to demonstrate Kolmogorov scaling with an inertial range
which spanned a decade of momenta. We plan on exploring the exact dependence
of turbulent features of the energy momentum tensor on the parameters in (31)
in future work as well. In the remainder of this section we present our
results for the stress tensor and the dual geometry associated with it given
the parameters in (31) and (32).

###  5.1 The boundary stress tensor

Once the stochastic metric is turned on, the components of the energy momentum
tensor $\hbox{Tr}\left(\varrho T^{\mu\nu}\right)$ quickly become dominated by
the driving force. A snapshot exhibiting typical behavior of the stress tensor
can be found in figure 1.

Figure 1:  A snapshot of the energy flux, $T_{01}(t,x^{1},x^{2})$ at
$t=750\rho_{0}$ where $\rho_{0}=4\pi T_{0}/3$ and $T_{0}$ is the temperature
prior to the onset of the driving force. The typical scale we observe is of
the same order as the driving scale, $1/k_{f}$.

After we ensemble averaged the expectation value of the stress tensor over 36
independent runs, we found that the magnitude of the oscillations resulting
from the driving force decreased by an order of magnitude. Carrying out a time
average over 19 snapshots, the oscillations decreased further. To demonstrate
this effect more clearly we have plotted the Fourier components of the stress
tensor and their angle averaged value in figure 2.

Figure 2:  A set of plots exhibiting the behavior of the Fourier decomposition
of the energy density (top) and energy flux (bottom), $\frac{16\pi
G_{N}}{3\ell^{2}}{T}_{00}(t,\,\vec{x})=\frac{1}{L^{2}\rho_{0}}\sum_{\vec{n}}E_{\vec{n}}e^{i\frac{2\pi\vec{n}}{L}\cdot\vec{x}}$
and $\frac{16\pi
G_{N}}{3\ell^{2}}{T}_{01}(t,\,\vec{x})=\frac{1}{L^{2}\rho_{0}}\sum_{\vec{n}}F_{\vec{n}}e^{i\frac{2\pi\vec{n}}{L}\cdot\vec{x}}$,
at $t=1125\rho_{0}$ with $\rho_{0}=4\pi T_{0}/3$ and $T_{0}$ the initial
temperature. Plots of $|E_{\vec{n}}|$ and $|F_{\vec{n}}|$ appear on the left
column. The centered ring is located at the driving scale, $k_{f}$. In the
right column we have plotted the (dimensionless) angle averaged energy density
$\hat{E}(|\vec{k}|)=\int
E_{\frac{\vec{k}L}{2\pi}}\frac{|\vec{k}|L}{2\pi}d\theta$ and flux
$\hat{F}(|\vec{k}|)=\int
F_{\frac{\vec{k}L}{2\pi}}\frac{|\vec{k}|L}{2\pi}d\theta$ (where $E_{\vec{n}}$
and $F_{\vec{n}}$ have been analytically continued to all values of $\vec{n}$)
for a single run (blue), an average over 36 independent runs (orange) and a
further average over 19 nearby time steps (pink). The dashed red line
corresponds to the driving scale $k_{f}$.

Apart from the clear signature of the driving force observed in figure 2, we
also note a possible underlying structure associated with momentum modes lower
than the driving scale. Momentum modes larger than the driving scale seem to
contain modes larger than those available on the torus and are therefore
unreliable. In future work we will treat these modes properly.

Turbulent flow in non relativistic fluids is characterized by power law
behavior of correlation functions of the velocity field, $\vec{v}$. In
particular, a hallmark of incompressible turbulent flow is that the velocity
field correlation function $\overline{v^{i}(\vec{x}+\vec{r})v^{j}(\vec{x})}$
satisfies

$\overline{v^{i}(\vec{x}+\vec{r})v^{j}(\vec{x})}=G(r)r^{i}r^{j}+H(r)\delta^{ij}$
(33a) with $G(r)\sim r^{-{\frac{4}{3}}}\qquad F(r)\sim r^{\frac{2}{3}}$ (33b)

in a certain inertial range of scales. This behavior implies that the energy
power spectrum defined via

$\hat{\epsilon}(t,\,k)=\frac{1}{(2\pi)^{2}}\int\frac{1}{2}\delta_{ij}\overline{\hat{v}^{i}\hat{v}^{j}}kd\theta$
(34)

with

$\hat{v}(t,\,\vec{k})=\int\vec{v}(t,\,\vec{x})e^{-i\vec{k}\cdot\vec{x}}d^{2}x$
(35)

satisfies the scaling law

$\hat{\epsilon}\sim k^{-\frac{5}{3}}\,,$ (36)

in the inertial range. To be slightly more precise, the power law behavior of
$\hat{\epsilon}$ behaves as in (36) for momentum scales smaller than the
forcing scale, $k_{f}$, and as $k^{-3}$ for momenta larger than the forcing
scale (assuming that the forcing scale is not too high). The former scaling
law is associated with what is usually referred to as an inverse energy
cascade and the latter to a direct enstrophy cascade. Such behavior of the
power spectrum for the velocity field is unique to $2+1$ dimensional flow and
is related to the existence of an accidentally conserved enstrophy charge. In
$3+1$ dimensional flow enstrophy is not conserved and one observers a power
law behavior of the form (36) only at momenta above the forcing scale—a direct
energy cascade.

The scaling law described in (36) is, strictly speaking, valid for non
relativistic fluid flow on a manifold with a Euclidean geometry. There are a
handful of works which attempt to generalize (36) to non relativistic systems
in such a way that the non relativistic limit of this generalization will
reduce to (36) Fouxon:2009rd ; Carrasco:2012nf ; Westernacher-
Schneider:2015gfa ; Westernacher-Schneider:2017snn ; Saridthesis . In
particular, the authors of Westernacher-Schneider:2017snn ; Saridthesis  have
numerically verified the predictions of Fouxon:2009rd ; Westernacher-
Schneider:2015gfa  for relativistic fluid flow at low Mach numbers.
Unfortunately, for the theories considered, low Mach number implies velocities
of at most $1/\sqrt{2}c$ at which relativistic corrections are still small. To
obtain truly relativistic turbulent flow, one would need to consider
supersonic flow which is notoriously challenging to simulate. It may well be
that relativistic turbulent flow is not universal Eyink:2017zfz . A study of
the relativistic entropy current can be found in Carrasco:2012nf ;
Marjieh:2020odp ; Pinzani-Fokeeva:2021klb .

Be that as it may, in order to compare our result for the expectation value of
the stress tensor to turbulent behavior, we identify the (Landau frame)
velocity field $u_{L}^{\mu}$ and proper energy density $\epsilon_{L}$ of the
fluid, via the eigenvalue equation

$\hbox{Tr}\left(\varrho
T^{\mu}{}_{\nu}\right)u_{L}^{\nu}=-\epsilon_{L}u_{L}^{\nu}\,,$ (37)

where the velocity is normalized such that $u_{L}^{\mu}u_{L\,\mu}=-1$. We
define the velocity field $\vec{v}$ via $v^{i}=u^{i}/u^{0}$. As can be
observed in figure 3 the velocities attained in a typical run are far from
being relativistic.

Figure 3:  Maximal values of the velocity field (relative to the speed of
light) as a function of time for a single run (no ensemble averaging).

Since the relativistic formalism developed in Fouxon:2009rd ; Carrasco:2012nf
; Eyink:2017zfz  reproduces the non relativistic result (36) in the non
relativistic limit, we will content ourselves with studying the non
relativistic power spectrum and its relation to the scaling behavior given in
(36).

At this point the astute reader might correctly point out that the non
relativistic energy power spectrum for the case at hand is not quite (34) due
to the stochastic component of the spatial metric. Still, the quantity
$\hat{\epsilon}$ defined in (34) will satisfy (36) due to (33). In what
follows, for lack of a better name, we will continue to refer to
$\hat{\epsilon}$ as the non relativistic energy power spectrum. A
straightforward computation of $\hat{\epsilon}$ yields a power law of the form
$k^{-\nu}$ with $\nu=1.53\pm 0.05$ in the inertial range
$5<\frac{kL}{2\pi}<14$. See figure 4.

Figure 4:  The non relativistic energy power spectrum defined in (34) as a
function of momenta for different times. The displayed data has been obtained
after ensemble averaging over 36 runs. The clearly visible maximum is due to
the driving force located at $k_{f}L/2\pi=21$, and has been clipped in order
to highlight the power law behavior at low momenta.

###  5.2 The bulk metric

As discussed in the introduction, one of the main goals of this program is to
understand how the turbulent behavior exhibited in figure 4 manifests itself
in terms of geometric quantities. As a general rule, perturbations on the
boundary of AdS are expected to inflate as one moves towards its interior due
to the AdS warp factor. In figure 5 we have plotted a snapshot of the behavior
of $F_{1}$ (the $t$, $x^{1}$ component of the bulk metric) as a function of
the spatial coordinates. We find that the stochastic behavior of the boundary
metric dominates the behavior of $F_{1}$ and its overall amplitude increases
as a function of the radial coordinate $\rho$.

Figure 5:  A stacked density plot of the $t$, $x^{1}$ component of the metric,
$F_{1}(\rho,t,\vec{x})$, as a function of the radial coordinate $\rho$ at time
$t=1125\rho_{0}$, ensemble averaged over 36 runs. As the radial coordinate
$\rho$ increases (moves farther from the boundary) the magnitude of the
perturbations induced by the driving force become larger. In each run there
exists an apparent horizon at $\rho=\rho_{h}(t,\vec{x})$ with
$0.9\lesssim\rho_{h}\lesssim 1.1$.

A more quantitative demonstration of the same feature can be found in figure 6
where we have plotted the Fourier components of $F_{1}$ as a function of
$\rho$.

Figure 6:  A plot of the Fourier decomposition of
$F_{1}(\rho,t,\vec{x})=\frac{1}{L^{2}}\sum_{\vec{n}}F_{1\,\vec{n}}(\rho,t)e^{i\frac{2\pi\vec{n}}{L}\cdot\vec{x}}$,
the $t$, $x^{1}$, component of the metric, averaged over angular directions,
$\hat{F}_{1}(\rho,t)=\frac{1}{\rho_{0}^{2}}\int
F_{1\,\frac{|\vec{k}|L}{2\pi}}(\rho,t)\frac{|\vec{k}|L}{2\pi}d\theta$, at time
$t=1125\rho_{0}$ and ensemble averaged over 36 runs. As the radial coordinate
$\rho$ increases (moves farther from the boundary) the magnitude of the
perturbations induced by the driving force become larger.

The area element of the apparent horizon follows a similar trend. Following
the construction of Chesler:2013lia , the location of the apparent horizon,
$\rho=\rho_{h}(t,\vec{x})$ may be determined by requiring that
$\rho-\rho_{h}=0$ is a null surface with zero expansion. On a typical run
$\rho_{h}$ will have features similar to those exhibited in the density plot
in figure 1 with $\rho_{h}$ within less than 10% of unity.

Once $\rho_{h}$ is known, the horizon area element can be read off of

$\mathcal{A}(t,\vec{x})=\sqrt{g_{h}(\rho_{h}(t,\vec{x}),t,\vec{x})}=\Sigma(\rho_{h}(t,\vec{x}),t,\vec{x})^{2}\,.$
(38)

As it turns out, $\Sigma$ is given to a good approximation by its boundary
value given in (27),

$\Sigma\sim\frac{e^{Q}}{\rho}\,,$ (39)

and the resulting Fourier decomposition of the horizon area element,

$\mathcal{A}=\frac{1}{L^{2}}\sum_{\vec{n}}\mathcal{A}_{\vec{n}}e^{i\frac{2\pi\vec{n}}{L}\cdot\vec{x}}\,,$
(40)

takes a form similar to that exhibited in on the left column of figure 2.

Instead of considering the Fourier components of $\mathcal{A}$, one may
consider a horizon area power spectrum defined by

$\hat{a}_{k}=\frac{1}{2\pi
L}\sum_{|\vec{n}|=\frac{kL}{2\pi}}|\Sigma_{\vec{n}}|^{2}\rightarrow\frac{1}{(2\pi)^{2}}\int\Big{|}\Sigma_{\frac{\vec{k}L}{2\pi}}\Big{|}^{2}kd\theta\,,$
(41)

where

$\Sigma=\frac{1}{L^{2}}\sum_{\vec{n}}\Sigma_{\vec{n}}e^{i\frac{2\pi\vec{n}}{L}\cdot\vec{x}}$
(42)

and the arrow denotes the continuum limit. The horizon area power spectrum is
related to the horizon area via

$A=\int\mathcal{A}d^{2}x=\frac{2\pi}{L}\sum_{n}\hat{a}_{\frac{2\pi
n}{L}}\rightarrow\int\hat{a}_{k}dk\,,$ (43)

similar to the relation between the velocity field and the non relativistic
energy power spectrum, c.f., (34). In the context of decaying turbulence, the
authors of Adams:2013vsa  observed that the horizon power spectrum follows the
non relativistic power spectrum to good accuracy. In Rozali:2017bll  it was
shown that (in the limit of a large number of dimensions) $\hat{a}$ is
proportional to $\hat{\epsilon}k^{2}$ but only for low Mach number. Here, we
make a similar observation. The apparent horizon power spectrum closely
follows $\hat{\epsilon}$ and obeys a scaling relation $\hat{a}\sim k^{-\nu}$
with $\nu=-0.28\pm 0.03$, see figure 7.

Figure 7:  A plot of the horizon area power spectrum, defined in (41), at time
$t=1250\rho_{0}$ as a function of momentum, ensemble averaged over 36 runs.
The fit to a $-5/3+2$ power law is inline with other indications for this
scaling described in the literature and discussed in the main text.

##  6 Discussion

In this work we have demonstrated that a holographic dual to a field theory
place on an external random metric is well defined provided that the driving
force is sufficiently differentiable in a sense described in detail in section
2. We then constructed a bulk dual to an initially thermalized state driven
out of equilibrium once the background metric fluctuated randomally. The
resulting energy momentum tensor exhibited a power spectrum which agrees
nicely with turbulent behavior, as exhibited in figure 4.

In figure 4 as well as figures 1 and 2 which exemplify the behavior of the
stress tensor, we observe that the dominant excitation of the fluid is the
driving force—it is several orders of magnitude larger than any other noteable
feature of the fluid. We expect that if one waits long enough then energy will
transfer from the driving scale to the lower momentum modes which will
dominate the power spectrum. Unfortunately, in the absence of large scale
friction, the energy density will eventually pile up in low momentum modes
spoiling the Kolmogorov scaling behavior observed in figure 4. A larger
inertial range, which neccessitates larger grid sizes, may ameliorate this
behavior.

In numerical simulations of $2+1$ dimensional, incompressible, flow, a
turbulent steady state is generated by adding large scale friction which
removes excess energy from low momentum modes. More precisely, energy is
pumped into the system at the driving scale, $1/k_{f}$, and gets dissipated at
large scales (by, say, friction with the boundary). Once a steady state is
reached a condensate will form. While it would be appealing to generate such a
condensate in a setup of the type discussed in this work, it is difficult to
envision how large scale friction may be incorporated in a natural way within
the framework of holography. One possibility is to induce appropriate boundary
conditions on the boundary manifold possibly following the line of reasoning
in Fujita:2011fp .

From the bulk point of view, the boundary excitations of the metric propagate
into the interior with an amplified amplitude due to the AdS warp factor. See
figure 5. Thus, if a signature of Kolmogorov scaling is apparent in the
components of the energy momentum tensor, it will also be apparent on the
event horizon. It is difficult to judge from the plots in figure 2 whether the
ensemble averaged energy density $\mathcal{T}_{00}$ exhibits an underlying
scaling behaviour. If we use the leading order hydrodynamic constitutive
relations to evaluate the energy density, and expand it around a low velocity
solution, we find

$\mathcal{T}_{00}=\overline{\epsilon
u_{0}u_{0}}=\overline{\epsilon}+\overline{\epsilon v^{2}}+\ldots$ (44)

reminiscent of the non relativistic energy density $\overline{v^{2}}$.
Further, in a conformal theory (and in the absence of additionally conserved
charges) the energy density $\epsilon$ is related to the entropy density $s$
via $s\propto\epsilon^{\frac{d-1}{d}}$ and the latter should be related to an
area power spectrum of the event horizon, similar to the one defined in (41).
Following Bhattacharyya:2008xc  (see also Adams:2013vsa ) the event horizon
and apparent horizon almost coincide as long as hydrodynamics is a good
approximation to the behavior of the energy momentum tensor. While the power
law behavior of the horizon area power spectrum is in line with the findings
of Adams:2013vsa ; Rozali:2017bll  for non relativistic theories, it is
unclear how or whether such a relation holds for relativistic flow, and
whether $\overline{\epsilon}$ captures any of its features. Unfortunately, our
current grid is too small to resolve these issues.

Increasing the grid size and the inertial range in an attempt to resolve the
behavior of the black brane in the vicinity of its horizon, is vital to
understanding the role of turbulence in black hole dynamics. However, before
doing so, we would like to step back and understand the Reynolds number, or
its analog, in the type of flows we are considering. In non relativistic,
incompressible, fluid flow the Reynolds number $Re$ is given by $Re=vX/\nu$
where $v$ is the typical scale for the velocity of the fluid, $X$ is a typical
length scale and $\nu$ is the kinematic viscosity. When the Reynolds number is
much larger than one, typically of order of thousands, the flow becomes
turbulent. The Reynolds number, or its equivalent in relativistic fluid flow
is less studied. In Fouxon:2009rd  it was suggested that in a relativistic
setting $Re=\frac{T\,X}{\eta/s}$ which is the ratio of the ideal to viscous
energy momentum tensor. Here $T$ is a typical temperature scale, $X$ a typical
length scale and $\eta/s$ the shear viscosity to entropy density ratio. Given
that $\eta/s$ is fixed in a holographic setting, we should be able to increase
the size of the inertial range by modifying $X$ or $T$. True universality
would imply that changing parameters associated with the driving force, such
as its amplitude $D$, its length scale $1/k_{f}$ or its decay time $\tau$
should not affect the flow. It would be satisfying if this feature or its
breaking can be demonstrated explicitly.

Even if the grid size is resolved and the ensemble averaged event horizon is
found to have, in some sense, a fractal dimension of the type advocated in
Adams:2013vsa  or a variant of it, c.f., Westernacher-Schneider:2017xie , its
geometric implication is still somewhat murky. It is unclear whether ensemble
averaging over classical black hole geometries gives a meaningful geometric
structure, or whether at all one should ensemble average over the event
horizon or, instead, consider the event horizon of the ensemble averaged
geometries. On a related note, perhaps it is possible to generate stochastic
noise by other means using, say, the constructions of Glorioso:2018mmw ;
deBoer:2018qqm  which generate the Schwinger Keldysh generating function,
naturally equipped with a stochastic term. We plan on dealing with these, and
other related issues in the near future.

## Acknowledgements

We thank A. Frishman and Y. Kafri for useful discussions. AY and SW are
supported in part by an Israeli Science Foundation excellence center grant
2289/18 and a Binational Science Foundation grant 2016324.

##  Appendix A Review of SDE’s

There is an extensive body of literature on stochastic differential equations
to which this extremely brief review can not do justice. In what follows we
focus on those features of SDE’s most relevant to the current work. We refer
the reader to the many textbooks on the subject, e.g., Oskendal ;
gardiner2009stochastic , for a more detailed exposition.

Consider the ordinary stochastic differential equation,

$\frac{\partial}{\partial t}X(t)=g(X(t))+h(X(t))\xi(t)\,,$ (45)

where $\xi(t)$ is a random variable. As it stands, the random nature of
$\xi(t)$, implies that (45) is ill defined. To see this, let us rewrite (45)
in integral form

$X(t)=X(t_{0})+\int_{t_{0}}^{t}g(X(t^{\prime}))dt^{\prime}+\int_{t_{0}}^{t}h(X(t^{\prime}))\xi(t^{\prime})dt^{\prime}\,.$
(46)

If we identify the integrals on the right hand side of (46) as Riemann-
Stieltjes integrals, then the second term on the right reads

$\displaystyle\begin{split}\int_{t_{0}}^{t}h(X(t^{\prime}))\xi(t^{\prime})dt^{\prime}=\lim_{\Delta
t\to 0}\sum_{n=0}^{N-1}h(X(\tau_{n}))\int_{t_{0}+n\Delta t}^{t_{0}+(n+1)\Delta
t}\xi(t^{\prime})dt^{\prime}\end{split}$ (47)

with $N\Delta t=t-t_{0}$ and $t_{0}+n\Delta t\leq\tau_{n}<t_{0}+(n+1)\Delta
t$. Intorducing $t_{n}=t_{0}+n\Delta t$ and denoting
$\tau_{n}=\cos^{2}\alpha_{n}t_{n+1}+\sin^{2}\alpha_{n}t_{n}$ we see that the
average,
$\overline{\int_{t_{0}}^{t}h(X(t^{\prime}))\xi(t^{\prime})dt^{\prime}}$, will
depend on $\alpha_{n}$ via the difference between
$\overline{X(t_{n})\xi(t^{\prime})}$ and
$\overline{X(t_{n+1})\xi(t^{\prime})}$ with $t_{n}<t^{\prime}<t_{n+1}$. Thus,
equation (45) is ill defined, and as opposed to ordinary (non stochastic)
differential equations it must be supplemented with a prescription for
carrying out the sum in (47).

One often used prescription for defining the integral in (47) is referred to
as an Itô prescription where

$\tau_{n}=t_{n}$ (48)

(or $\alpha_{n}=\pi/2$). Alternately, the Stratonovich prescription is given
by

$\displaystyle\begin{split}\int_{t_{0}}^{t}h(X(t^{\prime}))\xi(t^{\prime})dt^{\prime}=\lim_{\Delta
t\to
0}\sum_{n=0}^{N-1}h\left(\frac{X(t_{n+1})+X(t_{n})}{2}\right)\int_{t_{0}+n\Delta
t}^{t_{0}+(n+1)\Delta t}\xi(t^{\prime})dt^{\prime}\,,\end{split}$ (49)

(which is somewhat different from choosing $\alpha=0$ in (47) above.) Often,
one associates a stochastic term coming from physical processes with the
Stratonovich integral (49). In order to distinguish the two prescriptions we
will use

$\displaystyle\frac{\partial}{\partial t}X(t)$
$\displaystyle=g(X(t))+h(X(t))\xi(t)$ $\displaystyle\hbox{It\^{o}}\,,$ (50a)
$\displaystyle\frac{\partial}{\partial t}X(t)$
$\displaystyle=g(X(t))+h(X(t))\circ\xi(t)$
$\displaystyle\hbox{Stratonovich}\,,$ (50b)

in the remainder of this appendix. In the rest of this work we use the
Stratonovich scheme.

In many instances, the stochastic term $\xi(t)$ in (50) is associated with
white noise (with zero mean) where

$\overline{\xi(t)}=0\,,\qquad\overline{\xi(t^{\prime})\xi(t)}=D\delta(t-t^{\prime})\,,$
(51)

with $D>0$,222 Recall that the spectral density $S(\omega)$ is the Fourier
transform of the autocorrelation function
$R(t-t^{\prime})=\overline{\xi(t)\xi(t^{\prime})}=\int
S(\omega)e^{i\omega(t-t^{\prime})}d\omega$. (Assuming $\overline{\xi(t)}=0$.)
Thus, $\overline{\xi(t)^{2}}=\int S(\omega)d\omega$ and $S(\omega)d\omega$
specifies the contribution of frequencies in the range
$(\omega,\,\omega+d\omega)$ to the variance of $\xi(t)$. The power spectrum
for white noise is uniform. and higher moments of $\xi$ are such that

$\int_{t}^{t+T}\xi(t^{\prime})dt^{\prime}=\sqrt{DT}N(0,1)\,,$ (52)

with $N(\mu,\sigma)$ denoting a normal distribution with average $\mu$ and
standard deviation $\sigma$. In this case, one can relate the Itô and
Stratonovich differential equations via

$g+h\circ\xi=\left(g+\frac{D}{2}h^{\prime}h\right)+h\xi\,.$ (53)

To obtain (53) we made use of

$\displaystyle\int_{t}^{t+\Delta
t}h(X(t^{\prime}))\circ\xi(t^{\prime})dt^{\prime}$
$\displaystyle=h\left(\frac{X(t)+X(t+\Delta t)}{2}\right)\int_{t}^{t+\Delta
t}\xi(t^{\prime})dt^{\prime}\,,$ (54) $\displaystyle\int_{t}^{t+\Delta
t}h(X(t^{\prime}))\xi(t^{\prime})dt^{\prime}$
$\displaystyle=h\left(X(t)\right)\int_{t}^{t+\Delta
t}\xi(t^{\prime})dt^{\prime}\,,$ (55)

equations (52) and (50) and have expanded all expressions to
$\mathcal{O}\left((\Delta t)^{\frac{3}{2}}\right)$.

One somewhat unusual feature of the Itô integration scheme is that the usual
rules of calculus are modified. For instance, for a smooth function $f$, with

$f(X(t+\Delta t))=f(X(t))+f^{\prime}(X(t))(X(t+\Delta
t)-X(t))+\frac{1}{2}f^{\prime\prime}(X(t))(X(t+\Delta
t)-X(t))^{2}+\mathcal{O}\left((\Delta t)^{\frac{3}{2}}\right)\,,$ (56)

we obtain

$df(X)=f^{\prime}(X)\left(g(X)dt+h(X)\xi
dt\right)+\frac{D}{2}f^{\prime\prime}(X)h(X)^{2}dt$ (57)

where we made use of (52). In contrast, under the Stratonovich scheme,

$f(X(t+\Delta t))=f(X(t))+f^{\prime}\left(\frac{X(t+\Delta
t)+X(t)}{2}\right)(X(t+\Delta t)-X(t))+\ldots\,,$ (58)

and the usual rules of calculus hold,

$df(X)=f^{\prime}(X)\left(g(X)dt+h(X)\circ\xi dt\right)\,.$ (59)

With an integral expression for $X(t)$ one can derive a discrete time
evolution scheme by which one may evolve $X(t)$ numerically from some initial
time $t_{0}$. Expanding the Itô integral to linear order in $\Delta t$ we find
the Euler-Maruyama approximation

$X(t+\Delta t)=X(t)+h(X(t))\sqrt{\Delta tD}\phi+g(X(t))\Delta t\,,$ (60)

where $\phi$ is drawn from $N(0,1)$. A Stratonovich equation can be expanded
in a similar manner by using (53). When computing averages, the Euler Maruyama
scheme provides a good approximation to $X$ up to order $\Delta t$. But when
considering higher moments of $X$, the Euler-Maruyama scheme is good to order
$\sqrt{\Delta t}$. As in standard differential schemes there exist higher
order numerical methods which improve convergence in $\Delta t$. For instance,
the Milstein scheme is defined such that

$X(t+\Delta t)=X(t)+h(X(t))\sqrt{\Delta tD}\phi+g(X(t))\Delta
t+\frac{D}{2}h(X(t))h^{\prime}(X(t))\left(\phi^{2}-1\right)\Delta t\,,$ (61)

whose error is of order $\Delta t$, a factor of $\sqrt{\Delta t}$ higher than
the Euler-Maruyama scheme.

So far we have considered ordinary SDE’s. Our entire discussion can be
uplifted to a set of coupled ordinary SDE’s,

$\frac{\partial}{\partial t}\vec{X}=\vec{G}(\vec{X})+H(\vec{X})\vec{\xi}$ (62)

with $\vec{X}$ and $\vec{G}$, $n$ dimensional vectors, and $H$ an $n\times n$
matrix. In this work we will specialize to interactions of the form

$H_{ij}(\vec{X})=h_{i}(\vec{X})\delta_{jn}\,,$ (63)

such that a single stochastic source, $\xi$, is present. The Stratonovich
version of (62) with (63) can be recast in Itô form via an analog of (53)

$\vec{G}+\vec{h}\circ\xi=\left(\vec{G}+\frac{D}{2}\left(\vec{h}\cdot\vec{\nabla}_{X}\right)\vec{h}\right)+\vec{h}\xi\,.$
(64)

The Milstein scheme for solving the Itô SDE (62) now reads

$\vec{X}(t+\Delta t)=\vec{X}(t)+\vec{h}(t)\sqrt{D\Delta
t}\phi+\vec{G}(t)\Delta
t+\frac{D}{2}\left(\vec{h}\cdot\vec{\nabla}_{X}\right)\vec{h}(\phi^{2}-1)\Delta
t\,.$ (65)

## References

  * [1] A.N. Kolmogorov, _The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers_ , _C. R. Acad. Sci. URSS_ 30 (1941) 301.
  * [2] J.M. Maldacena, _The Large N limit of superconformal field theories and supergravity_ , _Int. J. Theor. Phys._ 38 (1999) 1113 [hep-th/9711200].
  * [3] E. Witten, _Anti-de Sitter space and holography_ , _Adv. Theor. Math. Phys._ 2 (1998) 253 [hep-th/9802150].
  * [4] S.S. Gubser, I.R. Klebanov and A.M. Polyakov, _Gauge theory correlators from noncritical string theory_ , _Phys. Lett. B_ 428 (1998) 105 [hep-th/9802109].
  * [5] S. Bhattacharyya, V.E. Hubeny, S. Minwalla and M. Rangamani, _Nonlinear Fluid Dynamics from Gravity_ , _JHEP_ 02 (2008) 045 [0712.2456].
  * [6] A. Adams, P.M. Chesler and H. Liu, _Holographic turbulence_ , _Phys. Rev. Lett._ 112 (2014) 151602 [1307.7267].
  * [7] C. Eling, I. Fouxon and Y. Oz, _Gravity and a Geometrization of Turbulence: An Intriguing Correspondence_ , 1004.2632.
  * [8] F. Carrasco, L. Lehner, R.C. Myers, O. Reula and A. Singh, _Turbulent flows for relativistic conformal fluids in 2+1 dimensions_ , _Phys. Rev. D_ 86 (2012) 126006 [1210.6702].
  * [9] S.R. Green, F. Carrasco and L. Lehner, _Holographic Path to the Turbulent Side of Gravity_ , _Phys. Rev. X_ 4 (2014) 011001 [1309.7940].
  * [10] K. Balasubramanian and C.P. Herzog, _Losing Forward Momentum Holographically_ , _Class. Quant. Grav._ 31 (2014) 125010 [1312.4953].
  * [11] C. Eling and Y. Oz, _Holographic Vorticity in the Fluid/Gravity Correspondence_ , _JHEP_ 11 (2013) 079 [1308.1651].
  * [12] C. Eling and Y. Oz, _The Anomalous Scaling Exponents of Turbulence in General Dimension from Random Geometry_ , _JHEP_ 09 (2015) 150 [1502.03069].
  * [13] J.R. Westernacher-Schneider, L. Lehner and Y. Oz, _Scaling Relations in Two-Dimensional Relativistic Hydrodynamic Turbulence_ , _JHEP_ 12 (2015) 067 [1510.00736].
  * [14] J.R. Westernacher-Schneider and L. Lehner, _Numerical Measurements of Scaling Relations in Two-Dimensional Conformal Fluid Turbulence_ , _JHEP_ 08 (2017) 027 [1706.07480].
  * [15] M. Rozali, E. Sabag and A. Yarom, _Holographic Turbulence in a Large Number of Dimensions_ , _JHEP_ 04 (2018) 065 [1707.08973].
  * [16] J.R. Westernacher-Schneider, _Fractal dimension of turbulent black holes_ , _Phys. Rev. D_ 96 (2017) 104054 [1710.04264].
  * [17] G.L. Eyink and T.D. Drivas, _Cascades and Dissipative Anomalies in Relativistic Fluid Turbulence_ , _Phys. Rev. X_ 8 (2018) 011023 [1704.03541].
  * [18] B. Chen, P.-C. Li, Y. Tian and C.-Y. Zhang, _Holographic Turbulence in Einstein-Gauss-Bonnet Gravity at Large $D$_, _JHEP_ 01 (2019) 156 [1804.05182].
  * [19] T. Andrade, C. Pantelidou, J. Sonner and B. Withers, _Driven black holes: from Kolmogorov scaling to turbulent wakes_ , 1912.00032.
  * [20] S. Bhattacharyya, R. Loganayagam, S. Minwalla, S. Nampuri, S.P. Trivedi and S.R. Wadia, _Forced Fluid Dynamics from Gravity_ , _JHEP_ 02 (2009) 018 [0806.0006].
  * [21] O. Aharony and V. Narovlansky, _Renormalization group flow in field theories with quenched disorder_ , _Phys. Rev. D_ 98 (2018) 045012 [1803.08534].
  * [22] S. de Haro, S.N. Solodukhin and K. Skenderis, _Holographic reconstruction of space-time and renormalization in the AdS / CFT correspondence_ , _Commun. Math. Phys._ 217 (2001) 595 [hep-th/0002230].
  * [23] P.M. Chesler and L.G. Yaffe, _Numerical solution of gravitational dynamics in asymptotically anti-de Sitter spacetimes_ , _JHEP_ 07 (2014) 086 [1309.1439].
  * [24] I. Fouxon and Y. Oz, _Exact Scaling Relations In Relativistic Hydrodynamic Turbulence_ , _Phys. Lett. B_ 694 (2010) 261 [0909.3574].
  * [25] A. Sarid, _Aspects of relativistic hydrodynamics_ , Master’s thesis, Technion, Israel, 2020.
  * [26] R. Marjieh, N. Pinzani-Fokeeva and A. Yarom, _Enstrophy from symmetry_ , 2009.03980.
  * [27] N. Pinzani-Fokeeva and A. Yarom, _Enstrophy without boost symmetry_ , 2103.06840.
  * [28] M. Fujita, T. Takayanagi and E. Tonni, _Aspects of AdS/BCFT_ , _JHEP_ 11 (2011) 043 [1108.5152].
  * [29] S. Bhattacharyya, V.E. Hubeny, R. Loganayagam, G. Mandal, S. Minwalla, T. Morita et al., _Local Fluid Dynamical Entropy from Gravity_ , _JHEP_ 06 (2008) 055 [0803.2526].
  * [30] P. Glorioso, M. Crossley and H. Liu, _A prescription for holographic Schwinger-Keldysh contour in non-equilibrium systems_ , 1812.08785.
  * [31] J. de Boer, M.P. Heller and N. Pinzani-Fokeeva, _Holographic Schwinger-Keldysh effective field theories_ , _JHEP_ 05 (2019) 188 [1812.06093].
  * [32] B. Oksendal, _Stochastic Differential Equations (3rd Ed.): An Introduction with Applications_ , Springer-Verlag, Berlin, Heidelberg (1992).
  * [33] C. Gardiner, _Stochastic Methods: A Handbook for the Natural and Social Sciences_ , Springer Series in Synergetics, Springer Berlin Heidelberg (2009).


2021

[]Kazuhiro Tamura

]Department of physics, Tokyo University of Science

# Analysis of shape change of droplet in dipolar Bose-Hubbard model

Shohei Watabe  Tetsuro Nikuni  [

###### Abstract

The long-range and anisotropic nature of the dipolar interaction provides the
so-called supersolid phases in Bose-Einstein condensates (BECs) in an optical
lattice. However, in a certain area of dipole interaction parameters, BECs can
form into a droplet. In this paper, in order to qualitatively understand the
droplet formations, we propose a toy model that allows us to estimate the size
and shape of droplets in dipolar Bose-Hubbard system in the optical lattice.
We compare results of the toy model with numerical solutions of the mean-field
calculation.

##  1 Introduction

Ultracold Bose gases have been studied to establish fundamental physics of
various quantum phenomena due to their flexible degrees of freedom, for
example, the controllability of the s-wave scattering length by the Feshbach
resonance. In the case of ultracold dipolar atoms, droplet formation is known
to be related to the supersolid states in which spatial translational symmetry
is brokenIlzhoefer2019 . Such supersolid states have been actively studied in
recent yearsL_onard_2017_nature ; L_onard_2017_science ; Morales_2018 . The
recently realized supersolid state in optical lattice Kovrizhin_2005  is well
described by the Bose-Hubbard model. In the earlier study Danshita_2009 ,
various ordered phases and unstable regions were found in the Bose-Hubbard
system with the dipolar interaction.

In this study, we establish a toy model describing droplets in the “unstable”
region of the dipole Bose-Hubbard system. We use our model to estimate the
size of droplets and discuss characteristic features of the shape and size of
droplets. We also compare the results of the toy model with the results of
numerical calculations and discuss the validity of the toy model. This
research provides a qualitative understanding of droplet formation and shape
control of droplets.

##  2 Ordered phases of Dipolar Bose-Hubbard model

The Hamiltonian of the dipolar Bose-Hubbard model is given by

$\displaystyle\hat{\mathcal{H}}=-\sum_{\braket{jk}}J_{jk}\hat{b}^{\dagger}_{j}\hat{b}_{k}+\frac{1}{2}\sum_{j}U_{j}\hat{b}^{\dagger}_{j}\hat{b}_{j}(\hat{b}^{\dagger}_{j}\hat{b}_{j}-1)+\frac{1}{2}\sum_{\braket{jk}}U_{jk}\hat{b}^{\dagger}_{j}\hat{b}_{j}\hat{b}^{\dagger}_{k}\hat{b}_{k},$
(1)

where $\braket{jk}$ is nearest neighbor sites. The last term is the dipolar
interaction between nearest neighbor sites.

The previous study Danshita_2009  showed the emergence of three quantum
phases, SF(superfluid), SSS(striped supersolid), and CSS(checker-board
supersolid), and an unstable region as a ground state of the two-dimensional
dipolar Bose-Hubbard model.

The energy of the dipolar Bose-Hubbard model in the mean-field approximation
is given by

$\displaystyle
E=-\sum_{\braket{jk}}J_{jk}\phi^{*}_{j}\phi_{k}+\frac{1}{2}U\sum_{j}\lvert\phi_{j}\rvert^{4}+\frac{1}{2}\sum_{\braket{jk}}U_{jk}\lvert\phi_{j}\rvert^{2}\lvert\phi_{k}\rvert^{2},$
(2)

where the order parameter is defined as an expectation value of the
annihilation operator $\phi_{j}=\braket{\hat{b}_{j}}$. The order parameter of
the three quantum phases can be written as

$\displaystyle\begin{cases}\phi_{j}^{\text{SF}}=\sqrt{n},\\\
\phi_{j}^{\text{SSS(x)}}=\sqrt{n}_{0}+\sqrt{n}_{1}\exp(i\bm{k}_{x}\cdot\bm{r}_{j}/2),\\\
\phi_{j}^{\text{CSS}}=\sqrt{n}_{0}+\sqrt{n}_{1}\exp(i\bm{k}_{xy}\cdot\bm{r}_{j}/2),\end{cases}$
(3)

where $\bm{k}_{x}=(1/a_{x},0,0),\bm{k}_{xy}=(1/a_{x},1/a_{y},0)$ with the
lattice constants $a_{x},a_{y}$, and $n$ is the average number of particles.
We can obtain a phase diagram of the dipolar Bose-Hubbard model by comparing
energies of the system in the three phases with the stability condition:

$\displaystyle\frac{\partial n}{\partial\mu}>0,$ (4)

where $n$ is given by $n_{0}+n_{1}$ for SSS and CSS. The chemical potential
$\mu$ is defined as

$\displaystyle\mu=\frac{\partial E}{\partial N},$ (5)

where $N$ is the total number of particles.

Figure 1 shows the phase diagram. The phase boundary of the phase diagram can
be determined by comparing the energy and analyzing the stability condition
(4). Previous study Danshita_2009  has shown that the particles are uniformly
distributed when the dipolar interaction is weak; the SSS state is achieved
when the dipolar interaction in one direction is stronger in the other
direction; the CSS state is achieved when the dipolar interactions in both
directions are repulsive and strong (Fig. 1).

We also determine the ground state of the system numerically using the
imaginary time evolution method without assuming the order parameter in the
form of Eq. (3). Figures 1-1 show the ground state obtained by numerical
calculations. The initial distribution of the imaginary time evolution method
is set as a random distribution with an average of 15 particles at each site.
Open boundary conditions are applied.

(a)

(b) SF

(c) SSS(x)

(d) SSS(y)

(e) CSS

(f) Unstable region

Figure 1: Density profiles of ordered phases of dipolar Bose-Hubbard model.
The hopping parameter is set as $J_{x}=0.1U,J_{y}=0.2U,J_{z}=0.05U$ and the
interaction strength is set as $U_{z}=-0.05U$. The density and size of the
system are set as $N/V=15$ and $V=16\times 16\times 16$. (1) Analytically
obtained phase diagram of ordered phases: (1) SF region $(U_{x,y}=0.025U)$,
(1) SSS(x) region $[(U_{x},U_{y})=(0.525,-0.225)U]$, (1) SSS(y) region
$[(U_{x},U_{y})=(-0.225,0.525)U]$, (1) CSS region
$[(U_{x},U_{y})=(0.575U,0.575U)]$ and (1) Unstable region $(U_{x,y}=-1.275U)$.

The SF state has a spatially uniform distribution; the SSS state has a
spatially periodic distribution in one direction; the CSS state has a
spatially periodic distribution in two directions. Actually, there are some
inhomogeneities because of the open boundary condition. The inhomogeneity is
less significant in CSS because, in this state, the occupied lattice site is
surrounded by vacant lattice sites in the $xy$ plane. In contrast, in SSS the
inhomogeneity appears along the stripe. The numerical calculations with
imaginary-time evolution show that droplet states emerge in the so-called
“unstable region”. These droplets are thin in the $z$ direction and spread
over a few sites in the $x$-$y$ plane. In later sections, we establish a toy
model of these droplet states and compare them with the numerical results.

##  3 Box-droplet model in dipolar Bose-Hubbard model

The droplets can emerge in the “unstable region” of the dipolar Bose-Hubbard
model. Here, we assume the “box-droplet” with a uniform density distribution
inside the droplet and no particles outside. The order parameter of the box
droplet can be written as

$\displaystyle\phi_{j}=\begin{cases}\displaystyle\sqrt{\frac{N}{l_{x}l_{y}l_{z}}}&\displaystyle\left(\frac{\|x_{j}\|}{\sigma_{x}},\frac{\|y_{j}\|}{\sigma_{y}},\frac{\|z_{j}\|}{\sigma_{z}}\leq
1\right)\\\ 0&\text{otherwise,}\end{cases},$ (6)

where $x_{j},y_{j},z_{j}$ are site indices and $N$ is the total number of
particles. We assumed that the center of the droplet is located at the origin
of the system $x_{j}=y_{j}=z_{j}=0$. $\sigma_{x},\sigma_{y},\sigma_{z}$ are
the “radius” of the droplet and $l_{x},l_{y},l_{z}$ are the “thickness” of the
droplet, which is defined as $l_{i}=2\sigma_{i}+1$ for $i=x,y,z$. The droplet
should have at least a single site thickness in each direction, and thus
$1\leq l_{x,y,z}$.

The mean-field energy of the box-droplet is given by

$\displaystyle E_{\mathrm{BOX}}=$
$\displaystyle-2\left(J_{x}\frac{l_{x}-1}{l_{x}}+J_{y}\frac{l_{y}-1}{l_{y}}+J_{z}\frac{l_{z}-1}{l_{z}}\right)N$
$\displaystyle\
+\frac{N^{2}}{2l_{x}l_{y}l_{z}}\left(U+2U_{x}\frac{l_{x}-1}{l_{x}}+2U_{y}\frac{l_{y}-1}{l_{y}}+2U_{z}\frac{l_{z}-1}{l_{z}}\right).$
(7)

If the system size is much greater than the size of the droplet, one can
neglect the first term of Eq.(3) as

$\displaystyle E_{\mathrm{BOX}}\simeq$
$\displaystyle\frac{N^{2}}{2l_{x}l_{y}l_{z}}\left(U+2U_{x}\frac{l_{x}-1}{l_{x}}+2U_{y}\frac{l_{y}-1}{l_{y}}+2U_{z}\frac{l_{z}-1}{l_{z}}\right).$
(8)

The width of the droplet is obtained by minimizing the energy (8). For
$\|U_{z}\|\ll U$, we obtain the following simple formulas for the droplet
width:

$\displaystyle l_{x}=\begin{dcases}1&\left(U_{y}<-U/2+2U_{x}\right),\\\
\frac{6U_{x}}{U+2(U_{x}+U_{y})}&\left(-U/2+2U_{x}<U_{y}<\frac{2U_{x}+U}{4}\right),\\\
\frac{4U_{x}}{U+2U_{x}}&\left(\frac{2U_{x}+U}{4}<U_{y}\right).\end{dcases}$
(9)

The analogous expression for $l_{y}$ can be obtained by exchanging $x$ and $y$
in Eq.(9).

The droplets described by this toy model can be roughly classified into two
types: rod-shaped droplet and rectangular droplet. The rod-shaped droplet
oriented in the $x$ direction has a single site thickness in the $y$ and $z$
directions and a width of more than one site in the $x$ direction. On the
other hand, the rectangular droplet has a single site thickness in the $z$
direction and has a width of more than one site in the $x$ and $y$ directions.
Comparing the energy of the box-droplet with that of the ordered phases, we
find that the droplet state has lower energy in the unstable region that
satisfies $\partial n/\partial\mu<0$. Figure 2 shows the phase diagram of the
dipolar Bose-Hubbard model including the droplet phase obtained by using the
box-droplet model. We also show the density profiles of droplet phases
obtained by the numerical calculation in Figs. 2 and 2.

(a)

(b) Rectangular droplet

(c) Rod-shaped droplet

Figure 2: Density profiles of droplets of dipolar Bose-Hubbard model. The
hopping parameters are set as $J_{x}=0.1U,J_{y}=0.2U,J_{z}=0.05U$ and the
interaction strength is set as $U_{z}=-0.05U$. The density and size of the
system are set as $N/V=15$ and $V=16\times 16\times 16$. (2) Phase diagram
with droplet phase, (2) Rectangular droplet with $U_{x,y}=-0.575U$ and (2)
Rod-shaped droplet with $(U_{x},U_{y})=(-0.575,0.575)U$.

##  4 Droplet width of box-droplet model

Figure 3 compares the width of the box-droplet with that of the numerical
calculation. The droplet width obtained by the box-droplet model is plotted in
Fig. 3(a) and the numerical result is shown plotted in Fig. 3(b). We also show
the droplet width estimated from the numerical results and that of our box-
droplet model in Table1.

(a) (a)Box-droplet (b)Numerical

Figure 3: Droplet width of the box-droplet model and numerical calculation. The hopping parameter and interaction strength are set as $J_{x}=0.1U,J_{y}=0.2U,J_{z}=0.05U$, and $U_{z}=-0.05U$. Droplet width of (a) the box-droplet model and (b) numerical calculation. |  | Box-droplet | Numerical | Droplet type
---|---|---|---|---
$U_{x}/U$ | $U_{y}/U$ | $l_{x}$ | $l_{y}$ | $l_{x}$ | $l_{y}$ |
-1.975 | 1.025 | 2.678 | 1.000 | 2.748 | 1.000 | Rod-shaped droplet(x)
-0.975 | 1.025 | 4.105 | 1.000 | 3.908 | 1.000 | Rod-shaped droplet(x)
-0.575 | 1.025 | 15.33 | 1.000 | 7.666 | 1.000 | Rod-shaped droplet(x)
-1.975 | -1.975 | 2.000 | 2.000 | 1.717 | 1.717 | Rectangular droplet
-0.975 | -0.975 | 2.017 | 2.017 | 2.037 | 2.639 | Rectangular droplet
-0.575 | -0.575 | 2.654 | 2.654 | 2.715 | 2.715 | Rectangular droplet
-0.275 | -0.275 | 16.50 | 16.50 | 7.367 | 7.354 | Rectangular droplet
Table 1: droplet width

In most parameter regions, the box-droplet model can accurately predict the
numerically calculated droplet width. Only in the narrow region near the
boundary between ordered and droplet phases, the numerically calculated width
of the droplet becomes very large and the box-droplet model becomes less
accurate. According to the numerical calculation, the Box-droplet model
describes the shape of droplets reasonably well when the droplet width is
smaller than 5. This condition corresponds to $U_{x}\leq-0.83U$ for the rod-
shaped droplet oriented in the $x$ direction and $U_{y}\leq-0.4U_{x}-0.5U$ for
the rectangular droplet. This difference between the numerical calculation and
the box-droplet model can be understood from the density distribution of the
droplet. The distribution is approximately uniform if the droplet is small
enough. In contrast, the density distribution cannot be considered uniform if
the droplet width is much greater than one site.

##  5 Conclusion

We have calculated the ground state of the dipolar Bose-Hubbard system by
using the imaginary time evolution method. As a result, we found that the
droplet exists in the “unstable region”, which is consistent with the previous
study. We constructed a box-droplet model that can qualitatively predict the
shape of droplets. The box-droplet model is less accurate for large droplets
where the density distribution cannot be considered uniform. In contrast, the
box-droplet model can accurately predict the droplet width for small droplets.

## Acknowledgment

This work is supported by JSPS KAKENHI Grants No. JP18K03499.

## References

  *   * (1) P. IlzhÃ¶fer, M. Sohmen, G. Durastante, C. Politi, A. Trautmann, G. Natale, G. Morpurgo, T. Giamarchi, L. Chomaz, M.J. Mark, and F. Ferlaino, Nat. Phys. 17, 356 (2021).
  * (2) J. Léonard, A. Morales, P. Zupancic, T. Esslinger, and T. Donner, Nature 543, 87 (2017).
  * (3) J. Léonard, A. Morales, P. Zupancic, T. Donner, and T. Esslinger, Science 358, 1415 (2017).
  * (4) A. Morales, P. Zupancic, J. Léonard, T. Esslinger, and T. Donner, Nat. Mater. 17, 686 (2018).
  * (5) D.L. Kovrizhin, G.V. Pai, and S. Sinha, Europhys. Lett. 72, 162 (2005).
  * (6) I. Danshita, and C.A.R. Sá de Melo, Phys. Rev. Lett. 103, 225301 (2009).