# Emergence of hydrodynamic spatial long-range correlations in nonequilibrium
many-body systems

Benjamin Doyon  Department of Mathematics, King’s College London, Strand,
London WC2R 2LS, U.K.  Gabriele Perfetto  Institut für Theoretische Physik,
Eberhard Karls Universität Tübingen, Auf der Morgenstelle 14, 72076 Tübingen,
Germany.  Tomohiro Sasamoto  Department of Physics, Tokyo Institute of
Technology, Ookayama 2-12-1, Tokyo 152-8551, Japan  Takato Yoshimura  All
Souls College, Oxford OX1 4AL, U.K.  Rudolf Peierls Centre for Theoretical
Physics, University of Oxford, 1 Keble Road, Oxford OX1 3NP, U.K.

###### Abstract

At large scales of space and time, the nonequilibrium dynamics of local
observables in extensive many-body systems is well described by hydrodynamics.
At the Euler scale, one assumes that each mesoscopic region independently
reaches a state of maximal entropy under the constraints given by the
available conservation laws. Away from phase transitions, maximal entropy
states show exponential correlation decay, and independence of fluid cells
might be assumed to subsist during the course of time evolution. We show that
this picture is incorrect: under ballistic scaling, regions separated by
macroscopic distances develop long-range correlations as time passes. These
correlations take a universal form that only depends on the Euler
hydrodynamics of the model. They are rooted in the large-scale motion of
interacting fluid modes, and are the dominant long-range correlations
developing in time from long-wavelength, entropy-maximised states. They
require the presence of interaction and at least two different fluid modes,
and are of a fundamentally different nature from well-known long-range
correlations occurring from diffusive spreading or from quasi-particle
excitations produced in far-from-equilibrium quenches. We provide a universal
theoretical framework to exactly evaluate them, an adaptation of the
macroscopic fluctuation theory to the Euler scale. We verify our exact
predictions in the hard-rod gas, by comparing with numerical simulations and
finding excellent agreement.

Introduction.— Finding universal laws that govern many-body extended systems
Bloch  _et al._ (2008); Eisert  _et al._ (2015); Gogolin and Eisert (2016);
D'Alessio  _et al._ (2016); Langen  _et al._ (2016) at large scales away from
equilibrium is a fundamental problem in physics. Hydrodynamics is arguably the
most far-reaching and successful such set of laws Spohn (1991). The largest
scale of hydrodynamics, the Euler scale, exists and is nontrivial as soon as
the system admits ballistic transport and interactions are on short enough
distances; in particular, the system must possess at least a few extensive
conserved quantities, and hydrodynamic modes are in one-to-one correspondence
with these. Euler hydrodynamics applies to a wide array of many-body systems,
including gases and fluids of interacting particles. A prominent example,
which came to the fore recently, is the theory of generalised hydrodynamics
(GHD) Castro-Alvaredo  _et al._ (2016); Bertini  _et al._ (2016); De Nardis
_et al._ (2018); Doyon and Yoshimura (2017); Doyon  _et al._ (2018a);
Bastianello  _et al._ (2022), which proposes a universal structure for the
Euler hydrodynamics of many-body integrable systems, and which has been shown
to correctly describe cold atomic gases constrained to one dimension Schemmer
_et al._ (2019); Møller  _et al._ (2021); Malvania  _et al._ (2021) and
experimentally accessible gases of solitons El (2021); Bonnemain  _et al._
(2022). The Euler hydrodynamics of a microscopic model only requires the
knowledge of basic aspects of the emergent degrees of freedom, such as their
fluid velocities and their static correlations in stationary states. From
these data, it makes a range of nontrivial physical predictions, including the
large-scale motion of local observables, correlations at large separations in
space-time Doyon and Spohn (2017a); Doyon (2018); Møller  _et al._ (2020), and
the large-deviation theory for long-time ballistic transport Doyon and Myers
(2019); Myers  _et al._ (2020); Perfetto and Doyon (2021); Perfetto and
Gambassi (2020).

Euler hydrodynamics is based on a simple extension of equilibrium
thermodynamics Spohn (1991); Doyon (2020): in every mesoscopic region, or
“fluid cell”, the many-body system is assumed to maximise its entropy, under
the local constraints provided by the extensive conserved charges
$\hat{Q}_{i}=\int{\rm d}x\,\hat{q}_{i}(x,t)$. Here $\hat{q}_{i}(x,t)$’s are
the microscopic densities, related to the currents $\hat{j}_{i}(x,t)$ via
continuity equations $\partial_{t}\hat{q}_{i}+\partial_{x}\hat{j}_{i}=0$; we
concentrate on one-dimensional systems for simplicity. Thus, in every fluid
cell, a different Gibbs state, or generalised Gibbs ensemble (in integrable
models) Vidmar and Rigol (2016), is reached:
$\langle\bullet\rangle_{\underline{\beta}}={\rm
Tr}[\mbox{exp}(-\sum_{i}\beta^{i}\hat{Q}_{i})\bullet]/Z$ characterized by
Lagrange parameters $\underline{\beta}$. “Mesoscopic” refers to a size, $L$,
which is much greater than microscopic sizes $\ell_{\rm micro}$ – such as
inter-particle distances and interaction ranges – but much smaller than the
size $\ell$ where finite variations of local observables can be seen,
$\ell_{\rm micro}\ll L\ll\ell$. Perhaps the most natural setup where these
principles apply is when an initial state is prepared in the presence of
external, long-wavelength fields. Maximising entropy with fields
$\beta^{i}(x/\ell)$ coupled to the conserved densities and varying on large
scales $\ell$, the state is

$\langle\bullet\rangle_{\ell}={\rm
Tr}\Big{(}\mbox{exp}\Big{[}-\sum_{i}\int{\rm d}x\,\beta^{i}_{\rm
ini}(x/\ell)\hat{q}_{i}(x)\Big{]}\bullet\Big{)}/Z.$ (1)

The emergent, slow dynamics from the initial state (1) is that induced by the
continuity equations for the mesoscopic densities
$\mathtt{q}_{i}(x,t)=\langle\hat{q}_{i}\rangle_{\underline{\beta}(x,t)}$ and
currents
$\mathtt{j}_{i}(x,t)=\langle\hat{\jmath}_{i}\rangle_{\underline{\beta}(x,t)}$
as measured in the fluid cells,

$\partial_{t}\mathtt{q}_{i}(x,t)+\mathsf{A}_{i}^{\leavevmode\nobreak\
j}(x,t)\partial_{x}\mathtt{q}_{j}(x,t)=0,$ (2)

with $\mathsf{A}_{i}^{\leavevmode\nobreak\
j}=\partial\mathtt{j}_{i}/\partial\mathtt{q}_{j}$ the flux Jacobian. The
mesoscopic densities are related to the generalized inverse temperatures
$\underline{\mathtt{q}}\leftrightarrow\underline{\beta}$ bijectively thanks to
positivity of the static covariance matrix
$\mathsf{C}_{ij}=-\partial\mathtt{q}_{i}/\partial\beta^{j}$.

Euler hydrodynamics asserts that the slow variation in space of the state (1)
induces a corresponding slow variation in time, such that the state keeps the
local equilibrium form; one may propose that at all times (1) correctly
describes the states, with $\beta^{i}_{\rm ini}(x)\to\beta^{i}(x,t)$ (see,
e.g., De Groot and Mazur (1962); Spohn (1991); Doyon (2020)). In (1),
correlations vanish exponentially with the distance, under very broad
conditions including non-zero temperature (finite $\beta^{i}$’s) and the lack
of phase transition Israel (2015); Bratteli and Robinson (1987) (which we
assume here). Thus, spatial correlations should vanish exponentially, even
during time evolution; indeed each fluid cell is at “local equilibrium”, and
entropy maximisation should occur independently in every fluid cell Doyon
(2018); Perfetto and Doyon (2021); Ruggiero  _et al._ (2019); Bouchoule and
Dubail (2022).

In fact, certain long-range, algebraic correlations are known to emerge in
non-equilibrium situations when conservation laws are present. This is well
studied for diffusive systems in non-equilibrium steady states (NESS) Spohn
(1983); Dorfman  _et al._ (1994); de Zárate and Sengers (2004): unbalanced
thermostats at the system’s boundaries lead to nonzero gradients, and
correlations between conserved densities at macroscopic distances decay as
(system size)-1. This is due to the breaking of detailed balance at the
diffusive scale and determined by viscous coefficients, and may be
quantitatively described by fluctuating hydrodynamics and macroscopic
fluctuation theory (see, e.g., Dorfman  _et al._ (1994); de Zárate and Sengers
(2004); Bertini  _et al._ (2007)). But what happens at the Euler scale, where
viscous effects are scaled down to zero size?

In NESS emerging from the partitioning protocol in systems of infinite size
Bernard and Doyon (2016), gradients vanish and correlations are weaker. The
strongest are found in integrable systems, including free particles, where
conserved density correlations decay as (distance)-2 because of
discontinuities in the occupation function of hydrodynamic modes De Nardis and
Panfil (2018); Doyon  _et al._ (2015). But this decay is too quick to
correlate Euler-scale fluid cells (see below).

We note that a similar situation occurs at zero temperature, under the
different physics of quantum fluctuations at Fermi points, and that a theory
for the transport of such weak algebraic correlations on top of moving fluids
is proposed in Ruggiero  _et al._ (2019, 2020) (in GHD). Very far from
equilibrium, stronger long-range correlations may develop: for instance,
global quantum quenches generate finite densities of entangled particles that
may propagate (diffusively or ballistically) and carry nontrivial entanglement
Calabrese and Cardy (2005); Alba and Calabrese (2017, 2018) and correlations
Lux  _et al._ (2014); Del Vecchio  _et al._ (2023). But entangled particle
production is not expected to occur in long-wavelength states, Eq. (1).

Up to now, there has been no prediction, observation or theory for eventual
long-range correlations emerging under ballistic scaling from (1). The
assumption of uncorrelated Euler-scale fluid cells, and that the form (1)
stays valid in time, has remained, and appears to play an important role in
recent studies of the evolution of correlations and fluctuations under
inhomogeneous conditions and nonlinear hydrodynamic response theory Doyon
(2018); Perfetto and Doyon (2021); Ruggiero  _et al._ (2019, 2020); Fava  _et
al._ (2021).

In this manuscript, we show that the assumption of uncorrelated Euler-scale
fluid cells is generically incorrect. We show that correlations of strength
$\propto\ell^{-1}$ develop dynamically from (1), at macroscopic
($\propto\ell$) times and distances, under generic conditions for systems
exhibiting ballistic transport. In particular, if $Q_{i}^{R}(\ell
t),\,Q_{j}^{R^{\prime}}(\ell t)$ are total charges lying on finite but
macroscopically large regions $R,\,R^{\prime}$ that are separated by a
macroscopic distance, $|R|,|R^{\prime}|\propto\ell,\,{\rm
dist}(R,R^{\prime})\propto\ell$, evaluated at macroscopic time $\ell t$, then
their covariance is large, $\langle Q_{i}^{R}(\ell t)Q_{j}^{R^{\prime}}(\ell
t)\rangle^{\rm c}\propto\ell$. This shows strong correlations between
separated cells. The appearance of ballistically scaled long-range
correlations at all macroscopic times is a general phenomenon which, to our
knowledge, has not been discussed before. It holds no matter the nature of the
system, quantum or classical, integrable or not, and is solely controlled by
its Euler hydrodynamics.

This phenomenon is not explained by the theories for diffusive long-range
correlations recalled above, as it does not depend on viscous coefficients or
phenomenological noise, and occurs in ballistic times $t\propto x$. It gives
the dominant correlations on large distances, beyond diffusive broadening and
of larger strength than the $1/x^{2}$ correlations due to occupation
discontinuities. It is not due to quasi-particle excitations, as it is a
universal hydrodynamic effect. By contrast, we show that the phenomenon occurs
in long-wavelength inhomogeneous state (as in (1)), only if the Euler
hydrodynamic theory is interacting, and only if it admits at least two
different fluid modes (with different velocities). Euler-scale long-range
correlations invalidate the assumption that on every time-slice a state such
as (1) is found. This thus calls for a new understanding of the principles of
Euler hydrodynamics, and a re-think of recent studies of hydrodynamic
nonlinear response and the evolution of correlations and fluctuations.

We quantify this phenomenon by proposing that the principle replacing
independent local entropy maximisation of fluid cells is that of relaxation of
fluctuations: local observables relax to fixed, non-fluctuating functions of
conserved densities, which themselves fluctuate. This is developed into a
universal theory, the ballistic macroscopic fluctuation theory (BMFT). The
BMFT is a hydrodynamic large-deviation theory, solely based on the emergent
Euler hydrodynamic data of the model, which characterises all fluctuations and
correlations at the ballistic hydrodynamic scale, including under fluid
motion.

For illustration, we study the paradigmatic hard-rod model of statistical
physics, which is simple enough to be amenable to high-accuracy numerical
simulations, yet truly interacting. We find that the model does indeed develop
long-range correlations, which are quantitatively in excellent agreement with
our theory.

Ballistic long-range correlations.— We show that correlations in the initial
state (1), between macroscopically separated observables $\hat{o}_{1}(\ell
x_{1},\ell t)$ and $\hat{o}_{2}(\ell x_{2},\ell t)$ at macroscopic times,
generically has strength $\ell^{-1}$. That is, the connected correlation
function $\langle\hat{o}_{1}\hat{o}_{2}\rangle^{\rm
c}:=\langle\hat{o}_{1}\hat{o}_{2}\rangle-\langle\hat{o}_{1}\rangle\langle\hat{o}_{2}\rangle$
has a nonzero Euler-scaling limit 111Technically, for Eq. (3) to be
generically true, on its right-hand side observables $o_{i}(\ell x_{i},\ell
t)$ should be replaced by their fluid-cell means by averaging on cells in
space-time of size $L$ around $\ell x_{i},\ell t$. :

$S_{\hat{o}_{1},\hat{o}_{2}}(x_{1},t;x_{2},t):=\lim_{\ell\to\infty}\ell\langle\hat{o}_{1}(\ell
x_{1},\ell t)\hat{o}_{2}(\ell x_{2},\ell t)\rangle^{\rm c}_{\ell}\neq 0,$ (3)

when $x_{1}\neq x_{2}$. If fluid cells were uncorrelated, the result would be
a delta function. We will demonstrate that there is an additional nonzero
contribution if three conditions are satisfied: ( a) the initial state is
inhomogeneous (so there is fluid motion); ( b) the system is interacting (its
hydrodynamic flux Jacobian $\mathsf{A}_{i}^{\leavevmode\nobreak\ j}$ depends
on the state); and ( c) the model has at least two hydrodynamic modes with
different fluid velocities ($\mathsf{A}_{i}^{\leavevmode\nobreak\
j}\not\propto\delta_{i}^{\leavevmode\nobreak\ j}$). The latter two conditions
are generic in many-body systems that admit ballistic transport. This includes
translation invariant Hamiltonian systems such as the anharmonic chain Spohn
(2014), integrable quantum spin chains Bertini  _et al._ (2016), and
integrable gases such as the one-dimensional interacting Bose gas Lieb and
Liniger (1963) and soliton gases El (2021). Taking
$\hat{o}_{1}=\hat{q}_{i},\,\hat{o}_{2}=\hat{q}_{j}$ and integrating over the
regions $R,R^{\prime}$, the result implies that $\langle Q_{i}^{R}(\ell
t)Q_{j}^{R^{\prime}}(\ell t)\rangle^{\rm c}\propto\ell$ as claimed above.

The result (3) means that strong algebraic correlations develop at large times
and distances. If the initial state has short-range correlation, which is the
case for (1), then as time passes, correlations grow, and at macroscopic times
(of the order of the variation scale of the initial state), they are present
on linearly growing macroscopic distances (set by the maximal fluid velocity),
and are of strength inversely proportional to the macroscopic scale. Their
shape on the macroscopic scale and the way they eventually decay in time,
i.e., the precise function $S_{\hat{o}_{1},\hat{o}_{2}}(x_{1},t;x_{2},t)$ –
depends on the model and the precise shape of the initial condition.

The physics of such correlations comes from nonlinear hydrodynamic response.
Correlated fluid modes are emitted from fluid cells in response, at nonlinear
order, to the change of the cells’ states in time. The long-range correlations
are the product of such coherent fluid modes of different velocities
travelling to macroscopic distances and scattering amongst each other.
Nonlinear response is possible only if interaction is present, and long-range
correlations can only occur, after the emission event, if fluid modes
propagate at different velocities; this explains the three conditions above.
This is a new mechanism, inherently nonlinear and distinct from the well-known
mechanism for time correlations due to hydrodynamic linear response Doyon
(2018); Møller  _et al._ (2020); Perfetto and Doyon (2021); Fava  _et al._
(2021), whereby Euler-scale correlations occur when a mode emitted by
$\hat{o}_{1}$ at time $t_{1}$ is probed by $\hat{o}_{2}$ at $t_{2}>t_{1}$. The
emission of coherent fluid modes parallels, but is different from, the
mechanism of entangled particle production found to give long-range
correlations in homogeneous quantum quenches both for diffusive Lux  _et al._
(2014) and ballistic Del Vecchio  _et al._ (2023) transport. Naturally, the
latter is not universal and does not require inhomogeneity or interaction (in
the above hydrodynamic sense), by contrast to the effect we have uncovered.

In order to prove and quantify the result, Eq. (3), we introduce the ballistic
macroscopic fluctuation theory. The BMFT, which applies to the ballistic
scale, is similar in spirit macroscopic fluctuation theory (MFT) Bertini  _et
al._ (2002, 2015); Bodineau and Derrida (2006); Derrida and Gerschenfeld
(2009); Krapivsky  _et al._ (2014, 2015); Mallick  _et al._ (2022); Derrida
(2007), for the diffusive scale. Technically, the long-range correlations that
occur at ballistic and diffusive scales have similar explanations in BMFT and
MFT, respectively: the non-locality of a “free-energy functional” emerging
after time evolution. The physics is nevertheless markedly different; the
former is not an effect of noise, and requires the presence of at least two
fluid modes, in contrast to the latter. For the full development of the BMFT,
see the companion manuscript Doyon  _et al._ (2022).

Example: the hard-rod gas.— Before developing the BMFT, we illustrate the
result on the hard-rod model, with a well-established Euler hydrodynamics
description Boldrighini  _et al._ (1983); Boldrighini and Suhov (1997); Doyon
and Spohn (2017b, a). The rods with unit length $\ell_{\rm micro}=a=1$ have
ordered positions $x_{i}(t)\in{\mathbb{R}}$ that move in time $t$ freely at
velocities $v_{i}$, except at collisions $x_{i+1}(t)-x_{i}(t)=1$ where
velocities are exchanged. We consider the rod density
$\hat{q}_{0}(x)=\sum_{i}\delta(x-x_{i})$ and evaluate
$S_{\hat{q}_{0},\hat{q}_{0}}(x,t;0,t)$ from (3). We take large but finite
values of $\ell$; we find that $\ell>100$ is sufficient for a $<10\%$ relative
accuracy. The initial state $\langle\cdots\rangle_{\ell}$ is a random
distribution of nonoverlapping rods, uniformly with velocities $v_{i}\in\\{\pm
1\\}$ and with the “bump” density profile

$\langle\hat{q}_{0}(x,0)\rangle_{\ell}=\frac{1+3e^{-(x/\ell)^{2}}}{3+3e^{-(x/\ell)^{2}}}\in[1/3,2/3].$
(4)

The particle density $\mathtt{q}_{0}(x,t)$ is shown in the inset of Fig. 1. At
the macroscopic time $t=0.5$ the initial bump splits into two
counterpropagating bumps. We observe that the correlation
$\langle\hat{q}_{0}(\ell x,\ell t)\hat{q}_{0}(0,\ell t)\rangle^{\rm c}_{\ell}$
decay as $1/\ell$ and find a finite result as $\ell\to\infty$ by multiplying
by $\ell$, in accordance with the definition in Eq. (3) of the Euler scaling
limit. Note how both positive and negative correlations occur, and how one
point, $x=0$, is not in the moving bump; these are not simply correlations
between the counter-propagating bumps. The comparison between the result
predicted by BMFT and the simulation is displayed in Fig. 1. Striking are both
the fact that the result is nonzero, thus long-range correlations are present,
and that the BMFT correctly predicts it. Correlations develop everywhere on
the region $x\in(-2,2)$, determined by the maximum fluid velocity, where the
large-scale motion takes place. Outside this interval, the state at time
$t=0.5$ still retains its homogeneous equilibrium structure and long-range
correlations are zero.

Figure 1: Bump release of the hard-rod gas. Rod density mean value
$\mathtt{q}_{0}(x,t)$ (in the inset) and equal-time connected correlation
function $S_{\hat{q}_{0},\hat{q}_{0}}(x,t;0,t)$ from (3) with
$\hat{q}_{0}(x)=\sum_{i}\delta(x-x_{i})$, evaluated at the macroscopic time
$t=0.5$ as a function of the macroscopic space coordinate $x$. The initial
state is the density bump in Eq. (4) (black solid line in the inset). Three
values $\ell=250,500,$ and $1000$ are reported with the associated uncertainty
bars (the data for the scales $\ell=250$ and $1000$ are shifted horizontally
by $\pm 0.04$ to make them clearly visible). The collapse of the numerical
data is evident, numerically confirming the existence of Euler-scale long-
range correlations of order $10^{-2}$. These correlations are captured by the
BMFT (red points) with an excellent agreement. The numerical data are obtained
by simulating the deterministic hard-rod dynamics and by averaging over a
large number of independent realizations of the rods initial positions and
velocities. In the inset, the numerical data (blue dotted curve) are obtained
with $\ell=250$ and match the solution of the Euler equation (2) (red solid
line).

BMFT.— As recalled, Euler-scale hydrodynamics is justified through local
relaxation: local-observable averages are functions of the densities
$\underline{\mathtt{q}}$ only, $\lim_{\ell\to\infty}\langle\hat{o}(\ell x,\ell
t)\rangle_{\ell}=\langle\hat{o}\rangle_{\underline{\beta}(x,t)}:=\mathtt{o}[\underline{\mathtt{q}}(x,t)]$
Spohn (1991); Doyon (2020). Applying this assumption to the current averages
$\mathtt{j_{i}}[\underline{\mathtt{q}}(x,t)]$, the Euler equation (2) follows.
The fundamental principle of the BMFT extends this assumption to fluctuations.
We assume that mesoscopic charge densities can be identified with classical
fluctuating variables $\hat{\underline{q}}(\ell x,\ell t):=\underline{q}(x,t)$
222 For systems with multiple conservation laws, one has a vector
$\underline{q}$ of fluctuating variables with components $q_{i}$. The symbols
$\hat{q}_{i}$ and $\mathtt{q}_{i}$ represent the associated operator and its
(generalized) Gibbs average, respectively.. Local relaxation of fluctuation
then states that mesoscopic fluctuating variables $\hat{o}(\ell x,\ell
t):=o(x,t)$ do not fluctuate independently but are fixed functions of
$\underline{q}(x,t)$. The functional dependence of $o(x,t)$ on
$\underline{q}(x,t)$ is shown to be completely determined by the (generalized)
Gibbs state: $o(x,t)=\mathtt{o}[\underline{q}(x,t)]$ Doyon  _et al._ (2022).
Applying this to currents $j_{i}(x,t)=\mathtt{j}_{i}[\underline{q}(x,t)]$, and
exploiting the conservation laws, the BMFT probability measure for the
fluctuating charge densities is obtained,

${\rm d}\mathbb{P}[\underline{q}(\cdot,\cdot)]={\rm
d}\mu[\underline{q}(\cdot,\cdot)]\,e^{-\ell\mathcal{F}[\underline{q}(\cdot,0)]}\delta[\partial_{t}\underline{q}+\partial_{x}\underline{\mathtt{j}}[\underline{q}]],$
(5)

where ${\rm d}\mu[\underline{q}(\cdot,\cdot)]$ is the flat measure for
functions on the space-time region $\mathbb{S}:=\mathbb{R}\times[0,T]$. This
implies that fluctuations of mesoscopic densities stem from those of the
initial state and are deterministically propagated in space-time according to
the Euler equation. Initial fluctuations of densities $\underline{q}(\cdot,0)$
are determined by $e^{-\ell\mathcal{F}[\underline{q}(\cdot,0)]}$; the large
deviation function for the state (1) is:

$\mathcal{F}[\underline{q}(\cdot,0)]\\!=\\!\int_{\mathbb{R}}{\rm
d}x\,\Big{(}\beta^{i}_{\mathrm{ini}}(x)q_{i}(x,0)\\!-\\!f[\underline{\beta_{\rm
ini}}(x)]\\!-\\!s[\underline{q}(x,0)]\Big{)},$ (6)

with $s[\underline{q}]$ and $f[\underline{\beta}]$ the entropy and free-energy
density, respectively. Equation (6) follows from statistical mechanics Derrida
(2007). It is also used in the MFT (for diffusive systems, where, however, the
measure (5) in space-time is different); see, e.g., Refs. Derrida and
Gerschenfeld (2009); Krapivsky  _et al._ (2014). One may consider a more
general initial state, e.g., itself with long-range correlations Jin  _et al._
(2021). The BMFT average is given by

$\langle\langle\bullet\rangle\rangle_{\ell}=\frac{1}{Z}\int_{(\mathbb{S})}{\rm
d}\mu[\underline{q}(\cdot,\cdot)]\,e^{-\ell\mathcal{F}[\underline{q}(\cdot,0)]}\delta(\partial_{t}\underline{q}+\partial_{x}\underline{\mathtt{j}}[\underline{q}])\bullet,$
(7)

and reproduces the Euler scale of $\langle\bullet\rangle_{\ell}$ at large
$\ell$. The delta function is represented via Laplace transform, leading to
the introduction of an auxiliary field $\underline{H}(x,t)$. In the
large-$\ell$ limit, all results are obtained by saddle-point calculations –
thus we have reduced the problem of Euler-scale correlations to a functional
minimisation problem. From Eq. (5), all BMFT predictions follow.

Note that both the BMFT and the MFT Bertini  _et al._ (2015) are large-
deviation theories based on an action formalism describing space-time
configurations of densities and currents. They, however, pertain to different
hydrodynamic scales, ballistic and diffusive, respectively, and as such they
cannot be derived from one another. In the companion manuscript Doyon  _et
al._ (2022), we also combine these and introduce a multiscale hydrodynamic
fluctuation theory describing both the ballistic and the diffusive scale.

Long-range correlations from the BMFT.— From the above we evaluate (3), for
conserved densities, using the BMFT measure (5) as

$S_{\hat{q}_{i_{1}},\hat{q}_{i_{2}}}(x_{1},t_{1};x_{2},t_{2})=\lim_{\ell\,\to\infty}\ell\langle\langle
q_{i_{1}}(x_{1},t_{1})q_{i_{2}}(x_{2},t_{2})\rangle\rangle^{\rm c}_{\ell}.$
(8)

We introduce a generating function $\exp(\ell\Lambda)$ with
$\Lambda=\lambda_{1}q_{i_{1}}(x_{1},t_{1})+\lambda_{2}q_{i_{2}}(x_{2},t_{2})$.
By saddle point, we have as $\ell\to\infty$,
$\langle\langle\exp(\ell\Lambda)\rangle\rangle_{\ell}\to\exp(-\ell\mathcal{F}_{\Lambda}[\underline{q^{*}}])$
where
$\mathcal{F}_{\Lambda}[\underline{q}(\cdot,\cdot)]=\mathcal{F}[\underline{q}(\cdot,0)]-\Lambda[\underline{q}(\cdot,\cdot)]$.
One calculates $S_{\hat{q}_{i_{1}},\hat{q}_{i_{2}}}(x_{1},t_{1};x_{2},t_{2})$
by taking derivatives with respect to $\lambda$ at $\lambda=0$. Dropping the
superscript $*$, $\underline{q}$ (no longer fluctuating) solves the saddle
point equations:

$\displaystyle H^{i}(x,0)$
$\displaystyle=\beta^{i}_{\mathrm{ini}}(x)-\beta^{i}(x,0),$ (9a)
$\displaystyle H^{i}(x,T)$ $\displaystyle=0,$ (9b)
$\displaystyle\partial_{t}\beta^{i}+\mathsf{A}_{j}^{\,\,i}[\underline{\beta}]\partial_{x}\beta^{j}$
$\displaystyle=0,$ (9c)
$\displaystyle\partial_{t}H^{i}+\mathsf{A}_{j}^{\,\,i}[\underline{\beta}]\partial_{x}H^{j}$
$\displaystyle=-\lambda\delta^{\leavevmode\nobreak\
i}_{i_{1}}\delta(x-x_{1})\delta(t-t_{1}),$ (9d)

with the boundary condition
$\beta^{i}(x\to\pm\infty,t)=\beta^{i}_{\mathrm{ini}}(x\to\pm\infty)$ and
$H^{i}(x\to\pm\infty,t)=0$. Manipulations of these BMFT equations generically
imply long-range correlations, giving, in particular,

$\displaystyle
S_{\hat{q}_{i_{1}},\hat{q}_{i_{2}}}(x_{1},t_{1};x_{2},t_{1})=\mathsf{C}_{i_{1}i_{2}}(x_{1},t_{1})\delta(x_{2}-x_{1})$
$\displaystyle-\partial_{\lambda}\left.\left(U_{\lambda}(t_{1},0)\underline{\beta_{\mathrm{ini}}}\right)^{i}(x_{2})\right|_{\lambda=0}\mathsf{C}_{i\,i_{2}}(x_{2},t_{1}),$
(10)

where $U_{\lambda}(t,t^{\prime})$ is the nonlinear time-evolution operator
associated with Eq. (9c),
$\underline{\beta}(t)=U_{\lambda}(t,t^{\prime})\underline{\beta}(t^{\prime})$.
The first term in (Emergence of hydrodynamic spatial long-range correlations
in nonequilibrium many-body systems) is the linear response contribution Doyon
(2018); Møller  _et al._ (2020); Perfetto and Doyon (2021); Fava  _et al._
(2021), the correlation within the fluid cell, supported at $x_{1}=x_{2}$. The
second term is generically nonzero for $x_{1}\neq x_{2}$ (the operator
$U_{\lambda}$ depends on $x_{1}$) and accounts for long-range correlations. It
vanishes when $t_{1}=0$ [$U_{\lambda}(0,0)=1$], as the initial state (1) has
no Euler-scale correlations; and also when the state is homogeneous
[$U_{\lambda}(t_{1},0)\underline{\beta_{\mathrm{ini}}}=\underline{\beta_{\mathrm{ini}}}$]
and the model is noninteracting [$U_{\lambda}$ does not depend on $\lambda$],
giving conditions (a) and (b) (see after Eq. (3)). The condition (c) follows
from hydrodynamic projections Spohn (1991); Kipnis and Landim (1999); Ayala
_et al._ (2018); Doyon and Spohn (2017a); Doyon (2018); Spohn (1982); Doyon
(2022), see Doyon  _et al._ (2022).

In general, it is challenging to bring the second term on the right-hand side
of Eq. (Emergence of hydrodynamic spatial long-range correlations in
nonequilibrium many-body systems) to a calculable form. Solving the BMFT
equations (9a)-(9d) is generically tricky as shock solutions, where entropy is
not conserved, may appear Bressan (2000); Doyon  _et al._ (2017). In
integrable models, this problem does not arise, as GHD is known to display no
shock solutions. In GHD, Eqs. (9a)-(9d) can be exactly solved using the method
of characteristics Doyon  _et al._ (2018b) leading to an expression for
$S_{\hat{q}_{i_{1}},\hat{q}_{i_{2}}}(x_{1},t_{1};x_{2},t_{1})$ in terms of
solutions to certain integral equations, see Doyon  _et al._ (2022). In the
case of the hard-rod model, this leads to the result reported in Fig. 1.

Conclusions.— We have shown that long-range correlations (3) generically
develop under ballistic scaling over time in many-body systems. Here (and the
companion manuscript Doyon  _et al._ (2022)), we have provided, to our
knowledge, the first observation and description of this universal phenomenon;
This happens under three generic conditions: inhomogeneity in initial
conditions, interaction, and the presence of more than one fluid velocity. It
is different from known long-range effects seen in diffusive NESS and in
global quenches. This result invalidates the fundamental tenet of Euler
hydrodynamics that at each time fluid cells undergo entropy maximization
independent of each other. By introducing the BMFT, the first hydrodynamic
fluctuation theory describing all correlations and fluctuations at the Euler
scale, we have evaluated the ballistic long-range correlations in the hard-rod
model observing excellent agreement with numerical simulations. The long-range
correlations we unveil are expected to have significant impact and open new
research directions in the field of correlations and fluctuations in
inhomogeneous fluids, which is a broad area of current active research. For
example, can we extend quantum GHD Ruggiero  _et al._ (2020) to introduce such
stronger correlations from nonlinear effects at zero or small temperatures? Do
they affect nonlinear response coefficients Fava  _et al._ (2021)? We also
expect that these long-range correlations can be experimentally observed using
quantum gas microscopes Wei  _et al._ (2022).

Acknowledgements—We are grateful to Bruno Bertini, Olalla Castro-Alvaredo,
Jacopo De Nardis, Fabian Essler, Tony Jin, Pierre Le Doussal, Adam Nahum,
Tibor Rakovszky, Paola Ruggiero, and Herbert Spohn for illuminating
discussions. The work of B.D. was supported by the Engineering and Physical
Sciences Research Council (EPSRC) under grants EP/W000458/1 and EP/W010194/1.
G.P. acknowledges support from the Alexander von Humboldt Foundation through a
Humboldt research fellowship for postdoctoral researchers. The work of T.S.
has been supported by JSPS KAKENHI Grant Nos. JP16H06338, JP18H03672,
JP19K03665, JP21H04432, JP22H01143. B.D., T.S. and T.Y. acknowledge
hospitality and support from the Galileo Galilei Institute, and from the
scientific program on “Randomness, Integrability, and Universality”. B.D. and
T.Y. acknowledge hospitality and support from the Isaac Newton Institute, and
from the program “Dispersive hydrodynamics: mathematics, simulation and
experiments, with applications in nonlinear waves”.

## References

  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Eisert  _et al._ (2015) J. Eisert, M. Friesdorf, and C. Gogolin, Nat. Phys. 11, 124 (2015).
  * Gogolin and Eisert (2016) C. Gogolin and J. Eisert, Rep. Prog. Phys. 79, 056001 (2016).
  * D'Alessio  _et al._ (2016) L. D'Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Adv. Phys. 65, 239 (2016).
  * Langen  _et al._ (2016) T. Langen, T. Gasenzer, and J. Schmiedmayer, J. Stat. Mech. Theory Exp. 2016, 064009 (2016).
  * Spohn (1991) H. Spohn,  _Large Scale Dynamics of Interacting Particles_ (Springer Berlin Heidelberg, 1991).
  * Castro-Alvaredo  _et al._ (2016) O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, Phys. Rev. X 6, 041065 (2016).
  * Bertini  _et al._ (2016) B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, Phys. Rev. Lett. 117, 207201 (2016).
  * De Nardis  _et al._ (2018) J. De Nardis, D. Bernard, and B. Doyon, Phys. Rev. Lett. 121, 160603 (2018).
  * Doyon and Yoshimura (2017) B. Doyon and T. Yoshimura, SciPost Phys. 2, 014 (2017).
  * Doyon  _et al._ (2018a) B. Doyon, T. Yoshimura, and J.-S. Caux, Phys. Rev. Lett. 120, 045301 (2018a).
  * Bastianello  _et al._ (2022) A. Bastianello, B. Bertini, B. Doyon, and R. Vasseur, J. Stat. Mech. Theory Exp. 2022, 014001 (2022).
  * Schemmer  _et al._ (2019) M. Schemmer, I. Bouchoule, B. Doyon, and J. Dubail, Phys. Rev. Lett. 122, 090601 (2019).
  * Møller  _et al._ (2021) F. Møller, C. Li, I. Mazets, H.-P. Stimming, T. Zhou, Z. Zhu, X. Chen, and J. Schmiedmayer, Phys. Rev. Lett. 126, 090602 (2021).
  * Malvania  _et al._ (2021) N. Malvania, Y. Zhang, Y. Le, J. Dubail, M. Rigol, and D. S. Weiss, Science 373, 1129 (2021).
  * El (2021) G. A. El, J. Stat. Mech. Theory Exp. 2021, 114001 (2021).
  * Bonnemain  _et al._ (2022) T. Bonnemain, B. Doyon, and G. El, J. Phys. A: Math. Theor. 55, 374004 (2022).
  * Doyon and Spohn (2017a) B. Doyon and H. Spohn, SciPost Phys. 3, 039 (2017a).
  * Doyon (2018) B. Doyon, SciPost Phys. 5, 54 (2018).
  * Møller  _et al._ (2020) F. S. Møller, G. Perfetto, B. Doyon, and J. Schmiedmayer, SciPost Phys. Core 3, 16 (2020).
  * Doyon and Myers (2019) B. Doyon and J. Myers, Ann. Henri Poincaré 21, 255 (2019).
  * Myers  _et al._ (2020) J. Myers, M. J. Bhaseen, R. J. Harris, and B. Doyon, SciPost Phys. 8, 7 (2020).
  * Perfetto and Doyon (2021) G. Perfetto and B. Doyon, SciPost Phys. 10, 116 (2021).
  * Perfetto and Gambassi (2020) G. Perfetto and A. Gambassi, Phys. Rev. E 102, 042128 (2020).
  * Doyon (2020) B. Doyon, SciPost Phys. Lect. Notes , 18 (2020).
  * Vidmar and Rigol (2016) L. Vidmar and M. Rigol, J. Stat. Mech. Theory Exp. 2016, 064007 (2016).
  * De Groot and Mazur (1962) S. De Groot and P. Mazur,  _Non-Equilibrium Thermodynamics_ (North Holland Pu. Co., Amsterdam, 1962).
  * Israel (2015) R. B. Israel,  _Convexity in the Theory of Lattice Gases_ (Princeton University Press, 2015).
  * Bratteli and Robinson (1987) O. Bratteli and D. W. Robinson,  _Operator Algebras and Quantum Statistical Mechanics 2_ (Springer, Berlin, 1987).
  * Ruggiero  _et al._ (2019) P. Ruggiero, Y. Brun, and J. Dubail, SciPost Phys. 6, 051 (2019).
  * Bouchoule and Dubail (2022) I. Bouchoule and J. Dubail, J. Stat. Mech.: Theory Exp. 2022, 014003 (2022).
  * Spohn (1983) H. Spohn, J. Phys. A: Math. 16, 4275 (1983).
  * Dorfman  _et al._ (1994) J. R. Dorfman, T. R. Kirkpatrick, and J. V. Sengers, Annu. Rev. Phys. Chem. 45, 213 (1994).
  * de Zárate and Sengers (2004) J. M. O. de Zárate and J. V. Sengers, J. Stat. Phys. 115, 1341 (2004).
  * Bertini  _et al._ (2007) L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, and L. C., arXiv:0705.2996 (2007).
  * Bernard and Doyon (2016) D. Bernard and B. Doyon, J. Stat. Mech. Theory Exp. 2016, 064005 (2016).
  * De Nardis and Panfil (2018) J. De Nardis and M. Panfil, Phys. Rev. Lett. 120, 217206 (2018).
  * Doyon  _et al._ (2015) B. Doyon, A. Lucas, K. Schalm, and M. Bhaseen, J. Phys. A Math. Theor. 48, 1 (2015).
  * Ruggiero  _et al._ (2020) P. Ruggiero, P. Calabrese, B. Doyon, and J. Dubail, Phys. Rev. Lett. 124, 140603 (2020).
  * Calabrese and Cardy (2005) P. Calabrese and J. Cardy, J. Stat. Mech. Theory Exp. 2005, P04010 (2005).
  * Alba and Calabrese (2017) V. Alba and P. Calabrese, PNAS 114, 7947 (2017).
  * Alba and Calabrese (2018) V. Alba and P. Calabrese, SciPost Phys. 4, 017 (2018).
  * Lux  _et al._ (2014) J. Lux, J. Müller, A. Mitra, and A. Rosch, Phys. Rev. A 89, 053608 (2014).
  * Del Vecchio  _et al._ (2023) G. D. V. Del Vecchio, B. Doyon, and P. Ruggiero, arXiv:2301.02326 (2023).
  * Fava  _et al._ (2021) M. Fava, S. Biswas, S. Gopalakrishnan, R. Vasseur, and S. Parameswaran, PNAS 118, e2106945118 (2021).
  * Note (1) Technically, for Eq. (3) to be generically true, on its right-hand side observables $o_{i}(\ell x_{i},\ell t)$ should be replaced by their fluid-cell means by averaging on cells in space-time of size $L$ around $\ell x_{i},\ell t$.
  * Spohn (2014) H. Spohn, J. Stat. Phys. 154, 1191 (2014).
  * Lieb and Liniger (1963) E. H. Lieb and W. Liniger, Phys. Rev. 130, 1605 (1963).
  * Bertini  _et al._ (2002) L. Bertini, A. D. Sole, D. Gabrielli, G. Jona-Lasinio, and C. Landim, J. Stat. Phys. 107, 635 (2002).
  * Bertini  _et al._ (2015) L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, and C. Landim, Rev. Mod. Phys. 87, 593 (2015).
  * Bodineau and Derrida (2006) T. Bodineau and B. Derrida, J. Stat. Phys. 123, 277 (2006).
  * Derrida and Gerschenfeld (2009) B. Derrida and A. Gerschenfeld, J. Stat. Phys. 137, 978 (2009).
  * Krapivsky  _et al._ (2014) P. L. Krapivsky, K. Mallick, and T. Sadhu, Phys. Rev. Lett. 113, 078101 (2014).
  * Krapivsky  _et al._ (2015) P. L. Krapivsky, K. Mallick, and T. Sadhu, J. Stat. Phys. 160, 885 (2015).
  * Mallick  _et al._ (2022) K. Mallick, H. Moriya, and T. Sasamoto, Phys. Rev. Lett. 129, 040601 (2022).
  * Derrida (2007) B. Derrida, J. Stat. Mech. Theory Exp. 2007, P07023 (2007).
  * Doyon  _et al._ (2022) B. Doyon, G. Perfetto, T. Sasamoto, and T. Yoshimura, arXiv:2206.14167 (2022).
  * Boldrighini  _et al._ (1983) C. Boldrighini, R. L. Dobrushin, and Y. M. Sukhov, J. Stat. Phys. 31, 577 (1983).
  * Boldrighini and Suhov (1997) C. Boldrighini and Y. M. Suhov, Commun. Math. Phys. 189, 577 (1997).
  * Doyon and Spohn (2017b) B. Doyon and H. Spohn, J. Stat. Mech. Theory Exp. 2017, 073210 (2017b).
  * Note (2) For systems with multiple conservation laws, one has a vector $\underline{q}$ of fluctuating variables with components $q_{i}$. The symbols $\hat{q}_{i}$ and $\mathtt{q}_{i}$ represent the associated operator and its (generalized) Gibbs average, respectively.
  * Jin  _et al._ (2021) T. Jin, T. Gautié, A. Krajenbrink, P. Ruggiero, and T. Yoshimura, J. Phys. A: Math. Theor. 54, 404001 (2021).
  * Kipnis and Landim (1999) C. Kipnis and C. Landim,  _Scaling Limits of Interacting Particle Systems_ (Springer Berlin Heidelberg, 1999).
  * Ayala  _et al._ (2018) M. Ayala, G. Carinci, and F. Redig, J. Stat. Phys. 171, 980 (2018).
  * Spohn (1982) H. Spohn, Ann. Phys. 141, 353 (1982).
  * Doyon (2022) B. Doyon, Commun. Math. Phys. 391, 293 (2022).
  * Bressan (2000) A. Bressan,  _Hyperbolic systems of conservation laws: the one-dimensional Cauchy problem_, Vol. 20 (Oxford University Press on Demand, 2000).
  * Doyon  _et al._ (2017) B. Doyon, J. Dubail, R. Konik, and T. Yoshimura, Phys. Rev. Lett. 119, 195301 (2017).
  * Doyon  _et al._ (2018b) B. Doyon, H. Spohn, and T. Yoshimura, Nucl. Phys. B 926, 570 (2018b).
  * Wei  _et al._ (2022) D. Wei, A. Rubio-Abadal, B. Ye, F. Machado, J. Kemp, K. Srakaew, S. Hollerith, J. Rui, S. Gopalakrishnan, N. Y. Yao, I. Bloch, and J. Zeiher, Science 376, 716 (2022).


# Maximizing the Set Cardinality of Users Scheduled for Ultra-dense uRLLC
Networks

Shiwen He, , Jun Yuan, Zhenyu An, , Yunshan Yi, Yongming Huang S. He and J.
Yuan are with the School of Computer Science and Engineering, Central South
University, Changsha 410083, China. S. He is also with the Purple Mountain
Laboratories, Nanjing 210096, China. (email: {shiwen.he.hn,
yuanjun}@csu.edu.cn). Z. An and Y. Yi are with the Purple Mountain
Laboratories, Nanjing 210096, China. (email: {anzhenyu,
yiyunshan}@pmlabs.com.cn). Y. Huang is with the National Mobile Communications
Research Laboratory, School of Information Science and Engineering, Southeast
University, Nanjing 210096, China. He is also with the Purple Mountain
Laboratories, Nanjing 210096, China. (email: huangym@seu.edu.cn).

###### Abstract

Ultra-reliability and low latency communication has long been an important but
challenging task in the fifth and sixth generation wireless communication
systems. Scheduling as many users as possible to serve on the limited time-
frequency resource is one of a crucial topic, subjecting to the maximum
allowable transmission power and the minimum rate requirement of each user. We
address it by proposing a mixed integer programming model, with the goal of
maximizing the set cardinality of users instead of maximizing the system sum
rate or energy efficiency. Mathematical transformations and successive convex
approximation are combined to solve the complex optimization problem.
Numerical results show that the proposed method achieves a considerable
performance compared with exhaustive search method, but with lower
computational complexity.

###### Index Terms:

Ultra-reliability low latency communications, user scheduling, beamforming,
non-convex optimization.

##  I Introduction

Ultra-reliability and low latency communication (uRLLC) requires $99.999\%$
reliability and $1$ ms latency [1] for transmitting a $32$ bytes packet.
Recently, discussion on the property of achievable rate with finite
blocklength transmission is one hot topic for uRLLC [2]. The authors of [3]
studied low bound of blocklength and applied it to power allocation with four
different downlink transmission schemes. The authors of [4] proposed a novel
perturbation-based iterative algorithm to jointly optimize the blocklength
allocation and the unmanned aerial vehicle’s location. The authors of [5]
obtained an analytical condition under which the minimum rate requirement of
any user can be satisfied. The authors of [6] consider three different
problems with the objective of maximizing the users’ minimum rate.

For multiuser multi-antenna communication systems, a popular topic is joint
beamforming (BF) and user scheduling (US), which could be formulated as non-
convex mixed integer programming problem. It was usually studied based on
conventional Shannon capacity [7]. The authors of [8] proposed a low
complexity BF and US scheme for the downlink multiple-input multiple-output
based non-orthogonal multiple access (MIMO-NOMA) system, but the non-convex
problem was handled by independently tracking two subproblems, namely, BF
scheme and greedy min-power US scheme, instead of jointly solving them. The
authors of [9] investigated the joint US and analog beam selection problem for
codebook-based massive MIMO downlink systems with hybrid antenna architecture.
Its optimization objective function was to maximizing the user rate with
allowable transmission power constraint and a $0-1$ diagonal digital precoding
matrix, but without considering the demand of minimum rate of each user. The
authors of [10] considered the joint coordinated BF and US for coordinated
multi-point enabled new radio in unlicensed spectrum networks, with also
aiming at maximizing the throughput instead of considering the rate demand for
each user. An alternative approach for user scheduling is formulated as a
series of decision-making problems, where the user relationship between the
previous and current time is considered. More recently, the domain knowledge
and deep learning [11, 12] are integrated and adopted to solve the problem of
user scheduling, but few studies investigated how to maximize the set
cardinality of scheduled users.

For ultra-dense uRLLC systems, compared to maximizing the system sum rate,
scheduling as many users as possible to be served on the same time-frequency
resource is more important under the maximum allowable transmission power and
the minimum rate requirement of each user. Motivated by these observations, in
the letter, we aim at scheduling as many serving users as possible via
optimizing the transmitting precoding simultaneously subject to the maximum
allowable transmission power and the minimum rate requirement of each user. A
resource allocation model which simultaneously considering US and BF is
proposed. Mathematical transformations and successive convex approximation are
combined to solve the original integer mixed optimization problem. Finally,
simulation results reveal its effectiveness.

##  II System Model and Problem Formulation

Consider a downlink multiple-input single-output uRLLC system with one BS and
$K$ single antenna users. BS is equipped with $N_{\mathrm{t}}$ antenna and
$K\geq N_{\mathrm{t}}$. Let $\mathcal{K}=\left\\{1,2,\cdots,K\right\\}$ and
$\mathcal{S}\subseteq\mathcal{K}$ be the set of all users and served users,
respectively. Let $\mathbf{w}_{k}\in\mathbb{C}^{N_{\mathrm{t}}\times 1}$,
$s_{k}$ and $\mathbf{h}_{k}\in\mathbb{C}^{N_{\mathrm{t}}\times 1}$
respectively denote the beamforming vector used by the BS, the baseband signal
for the $k$-th user, and the channel coefficient between the BS and the $k$-th
user.

The received baseband signal at the $k$-th user is expressed as:

$y_{k}=\sum\limits_{l\in\mathcal{S}}\mathbf{h}_{k}^{H}\mathbf{w}_{l}s_{l}+n_{k}{,}$
(1)

where $n_{k}$ is the zero-mean additive white Gaussian noise with variance
$\sigma_{k}^{2}$ at the $k$-th user. Thus, the signal-to-interference-plus-
noise ratio (SINR) of the $k$-th user is calculated as (2) with
$\overline{\mathbf{h}}_{k}=\frac{\mathbf{h}_{k}}{\sigma_{k}}$.

$\gamma_{k}=\frac{\left|\overline{\mathbf{h}}_{k}^{H}\mathbf{w}_{k}\right|^{2}}{\sum\limits_{l\neq
k,l\in\mathcal{S}}\left|\overline{\mathbf{h}}_{k}^{H}\mathbf{w}_{l}\right|^{2}+1}.$
(2)

The achievable rate obtained by taking the parameters of decoding error
probability $\epsilon$ and finite blocklength $n$ ( or $n$ channels used for
transmitting the user data) under consideration is described as [2, 6]

$R\left(\gamma_{k}\right)=C\left(\gamma_{k}\right)-\vartheta\sqrt{V\left(\gamma_{k}\right)},$
(3)

where $C\left(\gamma\right)=\ln\left(1+\gamma\right)$,
$\vartheta=\frac{Q^{-1}\left(\epsilon\right)}{\sqrt{n}}$ and
$Q^{-1}\left(\cdot\right)$ is the inverse of Gaussian Q-function
$Q\left(x\right)=\frac{1}{\sqrt{2\pi}}\int_{x}^{\infty}\mathrm{exp}\left(-\frac{t^{2}}{2}\right)dt$,
and $V\left(\gamma\right)$ is defined as

$V\left(\gamma\right)=1-\frac{1}{\left(1+\gamma\right)^{2}}.$ (4)

It is not difficult to find that compared to the classical Shannon
capacity111The proposed algorithm in the letter is also suitable for
communication with Shannon capacity., this is a more conservative description
of channel capacity. For the simplicity, we assume that the $n$ channels used
remain unchanged, that is, the channel model used is a slow time-varying
channel.

The set $\mathcal{S}$ of served users should be carefully scheduled in order
to effectively achieve the multiuser diversity gain and assure each user
obtain a minimum achievable rate i.e., $R\left(\gamma_{k}\right)\geq
r_{k}\geq\frac{D}{n}$, where $D$ represents transmitting data size. Note that
$D$ bits data is equally divided to $n$ channels since the slow time-varying
channel is assumed in the letter. To schedule as many users as possible under
the condition of the requirement of minimum user rate and the maximum
allowable transmitting power, the joint BF and US problem is formulated as
follows

$\displaystyle\max_{\mathcal{S}\subseteq\mathcal{K},\left\\{\mathbf{w}_{k}\right\\}}\left|\mathcal{S}\right|,$
(5a) $\displaystyle\mathrm{s.t.}$ $\displaystyle~{}r_{k}\leq
R\left(\gamma_{k}\right),\forall k\in\mathcal{S},$ (5b)
$\displaystyle\sum\limits_{k\in\mathcal{S}}\left\|\mathbf{w}_{k}\right\|_{2}^{2}\leq
P.$ (5c)

In (5a), $\left|\mathcal{S}\right|$ denotes the cardinality of set
$\mathcal{S}$ and $P$ is the maximum allowable transmission power of base
station. Note that problem (5) is non-convex problem and thus is hard to
obtain the solution directly. Furthermore, due to the uncertainty of user
scheduling set $\mathcal{S}$, problem (5) is essentially an uncertainty
problem. Generally speaking, brute exhaustive searching may be an effective
selection for solving the considered problem at the cost of high computational
complexity, especially for a large number of users to be scheduled.

##  III  Design of Optimization Algorithm

To release the uncertainty of user scheduling set $\mathcal{S}$ and obtain a
tractable form of problem (5), we introduce auxiliary variable $\kappa_{k}$ to
be an indicator of user scheduling. Specifically, if $\kappa_{k}=1$, the
$k$-th user is scheduled as a communication user, i.e., $k\in\mathcal{S}$.
Otherwise, $\kappa_{k}=0$ means that the $k$-th user is not scheduled as a
communication user, i.e., $k\notin\mathcal{S}$. Let
$\mathbf{\kappa}=\left[\kappa_{1},\kappa_{2},\cdots,\kappa_{k},\cdots,\kappa_{K}\right]^{T}$,
thus, we can rewrite equivalently problem (5) as follows

$\displaystyle\max_{\left\\{\kappa_{k},\mathbf{w}_{k}\right\\}}\left\|\mathbf{\kappa}\right\|_{0},$
(6a) $\displaystyle\mathrm{s.t.}$ $\displaystyle~{}\kappa_{k}r_{k}\leq
R\left(\gamma_{k}\right),\forall k\in\mathcal{K},$ (6b)
$\displaystyle\kappa_{k}\in\\{0,1\\},\forall k\in\mathcal{K},$ (6c)
$\displaystyle\sum\limits_{k\in\mathcal{K}}\left\|\mathbf{w}_{k}\right\|_{2}^{2}\leq
P.$ (6d)

Problem (6) is an equivalent form of problem (5) and the main goal of
objective functions (6a) is to assure that as many users as possible are
scheduled. Obviously, achievable rate $R\left(\gamma_{k}\right)$ in constraint
(6b) is non-convex nor non-concave, with a square root and inverse form of
$\gamma_{k}$. Besides, constraint (6c) make it be a mixed integer continuous
programming problem and also brings difficulties. Therefore, it is very
difficult to obtain the global optimal solution of problem (6), even the local
optimal solution.

In what follows, we focus on designing an effective iterative optimization
algorithm to solve problem (6) using a series of basic mathematical
operations. According to Theorem 2 and Corollary 1 in reference [5], the
required minimum SINR for satisfying the requirement of minimum user rate is
easily obtained with an analytical solution. Consequently, constraint (6b)
could be equivalently transformed to the following form:

$\displaystyle\max_{\left\\{\kappa_{k},\mathbf{w}_{k}\right\\}}\left\|\mathbf{\kappa}\right\|_{0},$
(7a) $\displaystyle\mathrm{s.t.}$
$\displaystyle~{}\kappa_{k}\widetilde{\gamma}_{k}\leq\gamma_{k},\forall
k\in\mathcal{K},~{}\eqref{URLLC06c},~{}\eqref{URLLC06d},$ (7b)

where $\widetilde{\gamma}_{k}$ is the minimum SINR for achieving the minimum
rate $r_{k}$ and can be obtained via using the conclusion of [5, Corollary 1].
Problem (7) can be transformed into the following form by simultaneously
substituting (2) to (7b) and introducing an auxiliary variable $\varphi_{k}$.

$\displaystyle\max_{\left\\{\kappa_{k},\varphi_{k},\mathbf{w}_{k}\right\\}}\sum\limits_{k\in\mathcal{K}}\kappa_{k},$
(8a) $\displaystyle\mathrm{s.t.}~{}$ $\displaystyle 0\leq\kappa_{k}\leq
1,\forall k\in\mathcal{K},$ (8b)
$\displaystyle\sum\limits_{k\in\mathcal{K}}\left(\kappa_{k}-\kappa_{k}^{2}\right)\leq
0,$ (8c)
$\displaystyle\kappa_{k}\widetilde{\gamma}_{k}-\frac{\left|\overline{\mathbf{h}}_{k}^{H}\mathbf{w}_{k}\right|^{2}}{\varphi_{k}}\leq
0,\forall k\in\mathcal{K},$ (8d) $\displaystyle\sum\limits_{l\neq
k,l\in\mathcal{K}}\left|\overline{\mathbf{h}}_{k}^{H}\mathbf{w}_{l}\right|^{2}+1\leq\varphi_{k},\forall
k\in\mathcal{K},~{}\eqref{URLLC06d},$ (8e)

where constraints (8b) and (8c) are equivalent to (6c). Actually, (8b) limits
$\kappa_{k}$ value between 0 and 1, and the inequality (8b) holds if and only
if every $\kappa_{k}$ is 0 or 1 (any $\kappa_{k}$ with value between 0 and 1
will make the summation be larger than 0). Note that the summation in
denominator of (2) is expanded from set $\mathcal{S}$ to $\mathcal{K}$.
According to [13, Proposition 2], the strong Lagrangian duality holds for the
above problem (8). Using similar mathematical tricks on handling
${(\sum\limits_{k\in{\cal K}}{{\kappa_{k}}})^{2}}$, we resort to solve problem
(8) via solving problem (9) by appropriately choosing $\mu>0$.

$\displaystyle\min_{\left\\{\kappa_{k},\varphi_{k},\mathbf{w}_{k}\right\\}}-\sum\limits_{k\in\mathcal{K}}\kappa_{k}+g\left(\mathbf{\kappa}\right)-h\left(\mathbf{\kappa}\right),$
(9a) $\displaystyle\mathrm{s.t.}~{}$
$\displaystyle\eqref{URLLC08b},~{}\eqref{URLLC08d},~{}\eqref{URLLC08e},~{}\eqref{URLLC06d},$
(9b)

where $g\left(\mathbf{\kappa}\right)$ and $h\left(\mathbf{\kappa}\right)$ are
defined respectively as

$\displaystyle g\left(\mathbf{\kappa}\right)$
$\displaystyle\triangleq\mu\sum\limits_{k\in\mathcal{K}}\kappa_{k}+\mu\left(\sum\limits_{k\in\mathcal{K}}\kappa_{k}\right)^{2},$
(10a) $\displaystyle h\left(\mathbf{\kappa}\right)$
$\displaystyle\triangleq\mu\sum\limits_{k\in\mathcal{K}}\kappa_{k}^{2}+\mu\left(\sum\limits_{k\in\mathcal{K}}\kappa_{k}\right)^{2}.$
(10b)

For fixed $\mu$, (9a) and (8d) are non-convex since the differences of convex
form, which brings difficulties for solving problem (9). Fortunately, note
that function $\kappa_{k}^{2}$ is convex, while function
$\frac{\left|\mathbf{h}_{k}^{H}\mathbf{w}_{k}\right|^{2}}{\varphi_{k}}$ is
quadratic-over-affine function, which is jointly convex w.r.t. the involved
variables.

In the sequel, we resort to solve problem (9) via using successive convex
approximation methods [14]. Using the convexity of function
$h\left(\mathbf{\kappa}\right)$, we have

$h\left(\mathbf{\kappa}\right)\geq
h\left(\kappa^{\left(\tau\right)}\right)+\sum\limits_{k\in\mathcal{K}}h^{\prime}\left(\kappa_{k}^{\left(\tau\right)}\right)\left(\kappa_{k}-\kappa_{k}^{\left(\tau\right)}\right)\triangleq\phi\left(\mathbf{\kappa}\right)$
(11)

where
$h^{\prime}\left(\kappa_{k}\right)=2\mu\left(\kappa_{k}+\left(\sum\limits_{k\in\mathcal{K}}\kappa_{k}\right)\right)$
and the superscript $\tau$ denotes the $\tau$-th iteration of the iterative
algorithm presented shortly. Note that $\phi\left(\mathbf{\kappa}\right)$ is
in fact the first-order Taylor series expansion approximation of function
$h\left(\mathbf{\kappa}\right)$ around the point
$\kappa_{k}^{\left(\tau\right)}$. Similarly, using the convex property of
function
$\frac{\left|\mathbf{h}_{k}^{H}\mathbf{w}_{k}\right|^{2}}{\varphi_{k}}$, we
can obtain their low boundary approximation as follows

$\begin{split}&\frac{\left|\mathbf{h}_{k}^{H}\mathbf{w}_{k}\right|^{2}}{\varphi_{k}}\geq\overline{\varphi}\left(\mathbf{w}_{k},\varphi_{k}\right)\\\
&\triangleq\frac{2\Re\left(\left(\mathbf{w}_{k}^{\left(\tau\right)}\right)^{H}\mathbf{h}_{k}\mathbf{h}_{k}^{H}\mathbf{w}_{k}\right)}{\varphi_{k}^{\left(\tau\right)}}-\left(\frac{\left|\mathbf{h}_{k}^{H}\mathbf{w}_{k}^{\left(\tau\right)}\right|}{\varphi_{k}^{\left(\tau\right)}}\right)^{2}\varphi_{k}.\end{split}$
(12)

From the aforementioned discussions, the convex approximation problem solved
at the $\left(\tau+1\right)$-th iteration of the proposed algorithm is given
by

$\displaystyle\min_{\left\\{\kappa_{k},\varphi_{k},\mathbf{w}_{k}\right\\}}-\sum\limits_{k\in\mathcal{K}}\kappa_{k}+g\left(\mathbf{\kappa}\right)-\phi\left(\mathbf{\kappa}\right),$
(13a) $\displaystyle\mathrm{s.t.}~{}$
$\displaystyle\eqref{URLLC08b},\eqref{URLLC08e},\eqref{URLLC06d},$ (13b)
$\displaystyle\kappa_{k}\widetilde{\gamma}_{k}-\overline{\varphi}\left(\mathbf{w}_{k},\varphi_{k}\right)\leq
0,\forall k\in\mathcal{K}.$ (13c)

Note that the objective function of problem (13) is an upper boundary of the
objective function of problem (9), while the feasible set of problem (13) is a
subset of the feasible set of problem (9). Consequently, the solution obtained
via problem (13) is an upper boundary of problem (9). Here, CVX [16] is
adopted to finally obtain the solution.

Algorithm 1 Solution of problem (13)

1: Let $\tau=0$. Initialize beamforming vector
$\mathbf{w}_{k}^{\left(0\right)}$, $\forall k\in\mathcal{K}$, $\mu=0.05$,
$\delta=10^{-3}$, and $\kappa_{k}^{\left(0\right)}$, $\forall
k\in\mathcal{K}$, such that constraint (7b) is satisfied.

2: Initialize $\varsigma^{\left(0\right)}$ and $\varphi_{k}^{\left(0\right)}$.

3: Let $\tau\leftarrow\tau+1$. Solve problem (13) to obtain
$\kappa_{k}^{\left(\tau\right)}$, $\varphi_{k}^{\left(\tau\right)}$, and
$\mathbf{w}_{k}^{\left(\tau\right)}$, $\forall k\in\mathcal{K}$.

4: Calculate objective value $\varsigma^{(\tau)}$. If
$\left|\frac{\varsigma^{\left(\tau\right)}-\varsigma^{\left(\tau-1\right)}}{\varsigma^{\left(\tau-1\right)}}\right|\leq\delta$,
stop iteration. Otherwise, go to step 3.

An iterative algorithm summarized as Algorithm 1 is outlined for solving
problem (13), where $\varsigma^{\left(\tau\right)}$ denotes the objective
value of optimization problem (13) at the $\tau$-th iteration iteration. To
speed up the convergence of Algorithm 1, we can first filter out the users who
meet constraints (5b) and (5c) by solving single user communication with
maximum ratio transmission and full power transmission.

(a) n = 128, $\epsilon$ = 1e-6.

(b) K = 8, $\epsilon$ = 1e-6.

(c) K = 8, n = 128.

Figure 1: Experiment results under different conditions.

######  Remark 1.

For convex problem (13), its optimal solution of the $\tau$-th iteration is
also a feasible point of the problem in the $(\tau+1)$-th iteration.
Therefore, Algorithm 1 is a non-increasing process. Besides, the power
constraint (6d) limits the lower bound of the problem and guarantees the
convergence[15].

######  Remark 2.

Parameter $\mu$ penalizes the objective function for any $\kappa_{k}$ with its
value between 0 and 1. A large value of $\mu$ will be helpful for convergence
of the algorithm, but may miss some users. In this letter, a smaller value of
$\mu$ would be chosen to find a better solution of (13) for the same initial
solution.

######  Remark 3.

An additional tuning process is adopted to solve the problem when
$\kappa_{k}=0,||\mathbf{w}_{k}||^{2}\geq 0$ even iteration stops, since the
residual power is larger than zeros but insufficient for serving more users.
The tuning process begins by choosing the $k$-th user with this situation.
Then its initial beamforming vector is set to be 0, and the whole optimization
steps 2-4 will be executed again. Here, algorithm 1 with and without tuning
operations are respectively referred as Algorithm 1 - with tuning and
Algorithm 1 - without tuning.

##  IV  Numerical Results

In this section, we compare different algorithms on simulated data. Simulated
channel follows the same model as in the reference [5], where channel
coefficient
${{\bf{h}}_{k}}=\sqrt{1/{{\left({1+({d_{k}}/{d_{0}}}\right)}^{\varrho}})}{{{\bf{\tilde{h}}}}_{k}}$.
Here, elements of ${{{\bf{\tilde{h}}}}_{k}}$ are assumed to be independent and
identically distributed, fading exponent $\varrho$ = 3, reference distance
$d_{0}$ = $30$ m, cell radius is $300$ m. All users have the same noise
variance, i.e., $\sigma_{k}^{2}=\sigma^{2}$, $\forall k\in\mathcal{K}$. SNR is
defined as $=10\log_{10}\left(\frac{P}{\sigma^{2}}\right)$ in dB for easy of
notation. Minimum rate for each user is $\frac{D}{n}$.Monte Carlo simulations
are implemented and the final results are all the mean values. Other main
parameters are listed as follows: $\sigma$ = $1$, antenna number
$N_{\mathrm{t}}$ = $4$, data size $D$ = 256 bits and SNR = $10$ dB. We compare
the proposed algorithm with exhaustive search (ES) method, with feasibility
criterion as in [5, IV-D], and the method using conventional Shannon capacity.
The latter problem could be realized via solving (13c) with replacing
${\tilde{\gamma}_{k}}$ with ${\tilde{\gamma}^{*}_{k}}={2^{\frac{D}{n}}}-1$.
Among them, computational complexity of ES method is
$O(\sum\limits_{\hat{k}=1}^{K}{C_{K}^{\hat{k}}{{(\hat{k}{N_{t}})}^{3}}(K-\hat{k}+1)\hat{k}})$,
where ${C_{K}^{\hat{k}}}$ is combinatorial number. The others share the same
lower computational complexity, namely, $O(({K}{N_{t}})^{3}(4K+1))$.

The number of scheduled users, latency and reliability are three crucial
factors in designing uRLLC system, which are respectively investigated in Fig.
1(a), 1(b) and 1(c). Among them, the reliability is described by decoding
error probability $\epsilon$, while the latency could be reflected by
blocklength $n$, since it is proportional to multiplication of blocklength $n$
and sampling interval $T_{s}$. In the figure, more users would be scheduled
with the increase of number of candidate users, $n$ and $\epsilon$ for the
four algorithms in general. Algorithm 1 with tuning obtains the closest
results compared with ES method, and its performance outperforms Algorithm 1
without tuning at the expense of computational complexity. It is also worth
noting that some scheduled users using conventional Shannon capacity might be
unavailable for transmitting, since the requirement of SINR for users are
relaxed in such circumstance. In the experiments, those users are removed from
scheduled user set before we plot figures. This also explains the phenomenon
that the results using Shannon capacity are much worse than the proposed
methods in Fig. 1(a) and Fig. 1(b).

##  V  Conclusions

In this letter, maximizing the set cardinality of users scheduled is
investigated for ultra-dense uRLLC networks. Instead of maximizing the system
sum rate with Shannon rate, a joint US and BF optimization model, which aims
at scheduling as many users as possible under the maximum allowable
transmission power and the minimum rate requirement of each user, is proposed
and effectively solved. Simulation results show the effectiveness of the
proposed algorithm. Note that the letter provides a novel view for resource
allocation in ultra-dense uRLLC networks, and the proposed algorithm could
also be extended to serving all users by using the Round Robin scheduling or
considering the user fairness problem. Nevertheless, how to schedule all users
with taking the fairness of users under consideration beyond the discussion of
this letter and it would be investigated in the future work.

## References

  * [1] 3GPP, “TR38.913 Study on scenarios and requirements for next generation access technologies (release 16),” _3rd Generation Partnership Project, Technical Specification Group Radio Access Network,_ Jul. 2020.
  * [2] Y. Polyanskiy, H. Poor, and S. Verdú, “Channel coding rate in the finite blocklength regime,” _IEEE Trans. Inf. Theory_ , vol. 56, no. 5, pp. 2307-2359, May 2010.
  * [3] H. Ren, C. Pan, Y. Deng, M. Elkashlan and A. Nallanathan, “Joint power and blocklength optimization for URLLC in a factory automation scenario,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 3, pp. 1786-1801, Mar. 2020.
  * [4] C. Pan, H. Ren, Y. Deng, M. Elkashlan and A. Nallanathan, “Joint blocklength and location optimization for URLLC-enabled UAV relay systems,” _IEEE Commun. Lett._ , vol. 23, no. 3, pp. 498-501, Mar. 2019.
  * [5] S. He, Z. An, J. Zhu, J. Zhang, Y. Huang, and Y. Zhang, “Beamforming design for multiuser uRLLC with finite blocklength transmission,” _to be published in IEEE Trans. Wireless Commun., 2021._ DOI:10.1109/TWC.2021.3090197.
  * [6] A. Nasir, H. Tuan, H. Nguyen, M. Debbah, and H. Poor, “Resource allocation and beamforming design in the short blocklength regime for URLLC,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 2, pp. 1321-1335, Feb. 2021.
  * [7] C. Shannon, “A mathematical theory of communication,” _Bell Syst. Tech. J,_ vol. 27, no. 7, pp. 379-423, Jul. 1948.
  * [8] C. Chen, W. Cai, X. Cheng, L. Yang and Y. Jin, “Low complexity beamforming and user selection schemes for 5G MIMO-NOMA systems,” _IEEE J. Select. Areas Commun._ , vol. 35, no. 12, pp. 2708-2722, Dec. 2017.
  * [9] Z. Jiang, S. Chen, S. Zhou and Z. Niu, “Joint user scheduling and beam selection optimization for beam-based massive MIMO downlinks,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 4, pp. 2190-2204, Apr. 2018.
  * [10] Q. Chen, K. Yang, H. Jiang and M. Qiu, “Joint beamforming coordination and user selection for CoMP enabled NR-U networks,” _IEEE Internet of Things J._ , Early Access, doi: 10.1109/JIOT.2021.3063733.
  * [11] C. She et al., “A tutorial on ultrareliable and low-latency communications in 6G: Integrating domain knowledge into deep learning,” _Proceedings of the IEEE_ , vol. 109, no. 3, pp. 204-246, March 2021.
  * [12] Z. Gu et al., “Knowledge-assisted deep reinforcement learning in 5G scheduler design: From theoretical framework to implementation,” _IEEE J. Sel. in Areas Commun._ , vol. 39, no. 7, pp. 2014-2028, July 2021.
  * [13] E. Che, H. D. Tuan and H. H. Nguyen, “Joint optimization of cooperative beamforming and relay assignment in multi-user wireless relay networks,” _IEEE Trans. Wireless Commun._ , vol. 13, no. 10, pp. 5481-5495, Oct. 2014.
  * [14] S. He, Y. Wu, D. W. K. Ng and Y. Huang, “Joint optimization of analog beam and user scheduling for millimeter wave communications,” _IEEE Commun. Lett._ , vol. 21, no. 12, pp. 2638-2641, Dec. 2017.
  * [15] L. Tran, M. Hanif, A. Tolli, and M. Juntti, “Fast converging algorithm for weighted sum rate maximization in multicell MISO downlink,” _IEEE Signal Process. Lett._ , vol. 19, no. 12, pp. 872-875, Dec. 2012.
  * [16] CVX: MATLAB Software for Disciplined Convex Programming. Jan, 2020. [Online]. Available: http://cvxr.com/cvx/