# Controllability maximization of large-scale systems using projected gradient
method

Kazuhiro Sato and Akiko Takeda K. Sato is with the Department of Mathematical
Informatics, Graduate School of Information Science and Technology, The
University of Tokyo, Tokyo 113-8656, Japan, email:
<EMAIL_ADDRESS>Takeda is with the Department of Creative
Informatics, Graduate School of Information Science and Technology, The
University of Tokyo, Tokyo 113-8656, Japan, and RIKEN Center for Advanced
Intelligence Project, 1-4-1, Nihonbashi, Chuo-ku, Tokyo 103-0027, Japan,
email<EMAIL_ADDRESS>
###### Abstract

In this work, we formulate two controllability maximization problems for
large-scale networked dynamical systems such as brain networks: The first
problem is a sparsity constraint optimization problem with a box constraint.
The second problem is a modified problem of the first problem, in which the
state transition matrix is Metzler. In other words, the second problem is a
realization problem for a positive system. We develop a projected gradient
method for solving the problems, and prove global convergence to a stationary
point with locally linear convergence rate. The projections onto the
constraints of the first and second problems are given explicitly. Numerical
experiments using the proposed method provide non-trivial results. In
particular, the controllability characteristic is observed to change with
increase in the parameter specifying sparsity, and the change rate appears to
be dependent on the network structure.

###### Index Terms:

Controllability, large-scale system, projected gradient method, sparsity

##  I Introduction

Controllability, which refers to the possibility to change the present network
state to a desired state is a fundamental concept in large-scale networked
dynamical systems [1, 2, 3, 4, 5, 6, 7, 8, 9]. Although several real networks
are nonlinear [10, 11, 12, 13], there are no general principles to determine
control inputs for a general nonlinear system to reach a desired state.
Moreover, in the context of neural systems, the assumption that a networked
system is linear is reasonable, because the baseline firing rates of the
neurons pertaining to the case with no stimulation are considerably lower than
those during stimulation and far from the saturation case [14]. Consequently,
several researchers [1, 2, 3, 4, 5, 6, 7, 8, 9] considered the following
linear system that is an approximation around an equilibrium point of a
nonlinear system:

$\displaystyle\dot{x}(t)=Ax(t)+Bu(t),$ (1)

where $x(t)\in{\mathbb{R}}^{n}$ and $u(t)\in{\mathbb{R}}^{m}$ denote the state
and input vectors, respectively; $A\in{\mathbb{R}}^{n\times n}$ is a fixed
constant matrix that reflects a network structure; and
$B\in{\mathbb{R}}^{n\times m}$ can be designed to suit the system
requirements. In other words, we can adjust the influence of input $u$ to the
autonomous dynamical networked system $\dot{x}(t)=Ax(t)$ by introducing a
decision variable $B$.

The authors in [15, 16, 17] considered a matrix $B_{S}$ instead of any $B$,
where $S\subset\\{1,2,\ldots,M\\}$ specifies the vectors from $M(\geq m)$
candidate column vectors of $B$, and addressed the following combinatorial
problem.

$\displaystyle\mathop{\rm
maximize}\limits_{S\subset\\{1,2,\ldots,M\\},\,|S|=m}\quad f(B_{S}),$ (2)

where $f(B_{S})$ denotes an index of controllability. In particular, in [15,
17], problem (2) was solved based on submodular optimization, whereas in [16],
problem (2) was transformed into a linear optimization problem and then solved
using a primal-dual distributed algorithm. The combinatorial problems related
to (2) were considered in [18, 19]. In addition, the authors of [20]
considered a controllability maximization problem in the form of a time-
varying actuator problem. Specifically, $\mathcal{B}V(t)$ was considered as
the matrix $B$, where $\mathcal{B}\in{\mathbb{R}}^{n\times m}$ and
$V(t)\in\\{0,1\\}^{m\times m}$ denote a fixed constant matrix and time-varying
diagonal matrix, respectively, and an optimization problem to determine the
diagonal entries of $V(t)$ was examined. This optimization problem was then
completely resolved, as reported in [21]. However, the methods to determine
the $M$ candidate column vectors of $B$ in [15, 16, 17] and the matrix
$\mathcal{B}$ in [20, 21] for a large-scale networked dynamical system remain
unclear.

To overcome this limitation, in this work, we consider the controllability
maximization problems from a different perspective than those considered in
[15, 16, 17, 20, 21]. In particular, we formulate two optimal decision
problems of $B$ in system (1) without using candidates of column vectors
considered in [15, 16, 17] and a fixed $\mathcal{B}$ considered in [20, 21].
The first problem is a sparsity constraint optimization problem with a box
constraint that specifies the possible values of each element in $B$. The
second problem is a modified problem of the first problem when $A$ is Metzler,
i.e., the off-diagonal elements of $A$ are nonnegative. Specifically, the
possible values of each element in $B$ in the second problem are restricted to
nonnegative values. The second problem is thus a realization problem for a
positive system, which is an important problem in the control community [22,
23].

The contributions of this work can be summarized as follows.

  * •

Although the considered problems are essentially combinatorial problems, we
solve the problems using a continuous optimization approach. In particular, we
propose a simple projected (that is, proximal) gradient method to solve the
considered problems. Moreover, we demonstrate the global convergence of the
proposed algorithm to a stationary point of the objective function.
Furthermore, considering a result recently reported in [24], it is shown that
the convergence rate of the sequences generated by the proposed algorithm is
locally linear.

  * •

The projections onto the constraints of the first and second problems are
given explicitly.

  * •

Numerical experiments using the proposed algorithm provide non-trivial
results. In particular, the controllability characteristic is observed to
change with increase in the parameter specifying sparsity, and the change rate
appears to be dependent on the network structure that determines the structure
of the matrix $A$.

The remaining paper is organized as follows. The problem formulation is
described in Section II. Section III describes the projected gradient methods
used to solve the problems and provides a proof of global convergence to a
stationary point with a convergence rate. Section IV describes the
experimental results, and the conclusions are presented in Section V.

Notation: The set of real numbers is denoted by ${\mathbb{R}}$. Given matrices
$A,B\in{\mathbb{R}}^{n\times m}$, we define $\langle A,B\rangle$ and
$\|A\|_{F}$ as the Euclidean inner product and the Frobenius norm,
respectively; i.e., $\langle A,B\rangle:={\rm tr}(A^{\top}B)$ and
$\|A\|_{F}:=\sqrt{\langle A,A\rangle}$, where the superscript $\top$ denotes
the transpose and ${\rm tr}(A)$ denotes the sum of the diagonal elements of
$A$. $\|A\|_{0}$ is defined as $l^{0}$ norm; i.e., $\|A\|_{0}$ denotes the
number of nonzero elements in $A$. For any matrices
$A,B\in{\mathbb{R}}^{n\times m}$, we write $A\geq B$ $(A\leq B)$ if all the
elements of $A$ are greater (less) than or equal to those of $B$. The symbol
${\bf E}\in{\mathbb{R}}^{n\times m}$ denotes a matrix whose elements are only
$1$. The symbol $I_{n}\in{\mathbb{R}}^{n\times n}$ denotes the identity
matrix.

##  II Problem settings

System (1) is termed controllable [25, 26] if for any desired final state
$x_{T}$ at any desired final time $T$, there exists an input $u$ such that
$x_{T}=\int_{0}^{T}\exp(A(T-t))Bu(t)dt$. That is, there exists an input $u$
satisfying $x(0)=0$ and $x(T)=x_{T}$. However, the input $u$ might be required
to have a high energy. In other words, even if system (1) is controllable, it
may be difficult to control the system state in practice. Thus, it is
important to consider a degree of controllability. In this section, we
formulate two controllability maximization problems using a controllability
index.

First, we consider a general case that $A$ in (1) is a fixed matrix that may
be unstable, and the controllability Gramian

$\displaystyle\mathcal{C}_{T}(B):=\int_{0}^{T}\exp(A(T-t))BB^{\top}\exp(A^{\top}(T-t))dt.$

The controllability Gramian can be related to the minimum-energy control
problem

$\displaystyle{\rm minimize}\quad\int_{0}^{T}\|u(t)\|^{2}dt$
$\displaystyle{\rm
subject\,\,to}\,\,\,\,\eqref{system},\,\,x(0)=0,\,\,x(T)=x_{T},$

where $x_{T}\in{\mathbb{R}}^{n}$ is any final state. In fact, if system (1) is
controllable, the minimum energy, i.e., the optimal objective value of the
above problem, is given by $x_{T}^{\top}\mathcal{C}^{-1}_{T}(B)x_{T}$, as
shown in [10]. Moreover, the Rayleigh-Ritz theorem implies that

$\displaystyle\frac{1}{\lambda_{\rm
max}(\mathcal{C}_{T}(B))}\leq\frac{x_{T}^{\top}\mathcal{C}^{-1}_{T}(B)x_{T}}{x_{T}^{\top}x_{T}}\leq\frac{1}{\lambda_{\rm
min}(\mathcal{C}_{T}(B))}$ (3)

for any $x_{T}\neq 0$, where $\lambda_{\rm min}(\mathcal{C}_{T}(B))$ and
$\lambda_{\rm max}(\mathcal{C}_{T}(B))$ denote the minimum and maximum
eigenvalues of $\mathcal{C}_{T}(B)$, respectively. Because inequality (3)
holds, $\lambda_{\rm min}(\mathcal{C}_{T}(B))$ and ${\rm
tr}(\mathcal{C}_{T}^{-1}(B))$, which denotes the sum of the inverse of all the
eigenvalues of $\mathcal{C}_{T}(B)$, are frequently adopted as the
controllability indices [1, 18, 17, 27]. Moreover, ${{\rm
tr}(\mathcal{C}_{T}^{-1}(B))}>\frac{n}{{\rm tr}(\mathcal{C}_{T}(B))}$. Thus,
to decrease ${\rm tr}(\mathcal{C}_{T}^{-1}(B))$, ${\rm
tr}(\mathcal{C}_{T}(B))$, which is the sum of all eigenvalues of
$\mathcal{C}_{T}(B)$, must be increased. If ${\rm tr}(\mathcal{C}_{T}(B))$ is
sufficiently large, there exists a direction that can be specified by an
eigenvector of $\mathcal{C}_{T}(B)$ such that the state $x(t)$ of system (1)
can move to the direction with a low input energy. Consequently, in the
existing studies [16, 20, 21], ${\rm tr}(\mathcal{C}_{T}(B))$ was adopted as a
controllability index. Note that unlike ${\rm tr}(\mathcal{C}_{T}^{-1}(B))$,
${\rm tr}(\mathcal{C}_{T}(B))$ can be defined even if system (1) is not
controllable.

In this study, we use ${\rm tr}(\mathcal{C}_{T}(B))$ as a controllability
index and consider the following modified problems:

$\displaystyle{\rm minimize}\quad h(B):=-{\rm tr}(\mathcal{C}_{T}(B))$ (4)
$\displaystyle{\rm subject\,\,to}\,\,\,\,\|B\|_{0}\leq s,$

where $s$ is a specified nonnegative value. The constraint $\|B\|_{0}\leq s$
ensures that the nonzero elements of $B$ are less than or equal to $s$. That
is, we can determine sparsity of $B$ by specifying $s\in\\{1,2,\ldots,nm\\}$.

However, in general, $h(B)$ is not bounded below, subject to $\|B\|_{0}\leq
s$. That is, a global optimal solution for (4) does not exist. To demonstrate
this aspect, we consider any $\beta>0$ and $\tilde{B}\in{\mathbb{R}}^{n\times
m}$. Because $h(\beta\tilde{B})=\beta^{2}h(\tilde{B})$ and
$\|\beta\tilde{B}\|_{0}=\|\tilde{B}\|_{0}$, $h(\tilde{B})<0$ (this relation
holds if system (1) with $B=\tilde{B}$ is controllable) and
$\|\tilde{B}\|_{0}\leq s$ imply that
$\lim_{\beta\rightarrow\infty}h(\beta\tilde{B})=-\infty$ and
$\lim_{\beta\rightarrow\infty}\|\beta\tilde{B}\|_{0}\leq s$.

To guarantee the existence of an optimal solution, we consider the following
problem with a box constraint.

Problem 1: Given any matrix $A\in{\mathbb{R}}^{n\times n}$, $T>0$, and
$s\in\\{1,2,\ldots,nm\\}$, find $B\in{\mathbb{R}}^{n\times m}$ that solves

$\displaystyle{\rm minimize}\quad h(B)$ $\displaystyle{\rm
subject\,to}\quad\|B\|_{0}\leq s,\,\,-{\bf E}\leq B\leq{\bf E}.$

The above discussion implies that if $B^{*}$ is an optimal solution to Problem
1, $\beta B^{*}$ is that to a modified problem in which $-{\bf E}\leq
B\leq{\bf E}$ in Problem 1 is replaced with $-\beta{\bf E}\leq B\leq\beta{\bf
E}$, where $\beta>0$.

Next, we consider a positive case in which $A$ in (1) is a fixed Metzler
matrix that may be unstable. Positive systems with Metzler matrix $A$ and non-
negative matrix $B$ are important, as mentioned in Section I. To realize a
positive system when $A$ is Metzler, we consider the following problem.

Problem 2: Given any Metzler matrix $A\in{\mathbb{R}}^{n\times n}$, $T>0$, and
$s\in\\{1,2,\ldots,nm\\}$, find $B\in{\mathbb{R}}^{n\times m}$ that solves

$\displaystyle{\rm minimize}\quad h(B)$ $\displaystyle{\rm
subject\,to}\quad\|B\|_{0}\leq s,\,\,0\leq B\leq{\bf E}.$

That is, we replace the constraint $-{\bf E}\leq B\leq{\bf E}$ in Problem 1
with $0\leq B\leq{\bf E}$.

###### Remark 1

The $l^{0}$ norm constraint is frequently replaced with the $l^{1}$ norm
constraint when an objective function and other constraints are convex,
because the modified problem then becomes convex [28, 29]. However, our
objective function $h$ is non-convex, as shown in Section III. Thus, even if
we replace $l^{0}$ with $l^{1}$, the modified problems are non-convex. Hence,
in this study, we do not replace $l^{0}$ with $l^{1}$.

###### Remark 2

The objective values at the global optimal solutions to Problems 1 and 2
monotonically increase as $s$ increases. That is, when we use the global
optimal solutions, controllability index $-h(B)$ increases as $s$ increases.

##  III Projected gradient methods for Problems 1 and 2

In this section, we develop projected gradient methods for solving Problems 1
and 2. To this end, we consider

$\displaystyle{\rm minimize}\quad\mathcal{I}_{Z_{s}\cap W}(B)+h(B),$ (5)
$\displaystyle{\rm minimize}\quad\mathcal{I}_{Z_{s}\cap W_{[0,1]}}(B)+h(B),$
(6)

where (5) and (6) are unconstrained problems equivalent to Problems 1 and 2,
respectively,

$\displaystyle\mathcal{I}_{\mathcal{S}}(B):=\begin{cases}0\quad\,\,\,(B\in\mathcal{S})\\\
\infty\quad(B\not\in\mathcal{S})\end{cases}$

is the indicator function of any set $\mathcal{S}$, and

$\displaystyle Z_{s}$ $\displaystyle:=\\{B\in{\mathbb{R}}^{n\times
m}\,|\,\|B\|_{0}\leq s\\},$ $\displaystyle W$
$\displaystyle:=\\{B\in{\mathbb{R}}^{n\times m}\,|\,-{\bf E}\leq B\leq{\bf
E}\\},$ $\displaystyle W_{[0,1]}$ $\displaystyle:=\\{B\in{\mathbb{R}}^{n\times
m}\,|\,0\leq B\leq{\bf E}\\}.$

We first note that

$\displaystyle h(B)$ $\displaystyle=-\int_{0}^{T}{\rm
tr}(\exp(A(T-t))BB^{\top}\exp(A^{\top}(T-t)))dt$
$\displaystyle=-\int_{0}^{T}{\rm
tr}(BB^{\top}\exp(A^{\top}(T-t))\exp(A(T-t)))dt$ $\displaystyle={\rm
tr}(BB^{\top}H(A,T)),$ (7)

where

$\displaystyle H(A,T):=-\int_{0}^{T}\exp(A^{\top}(T-t))\exp(A(T-t))dt.$

Using (7), we prove the following theorem.

###### Theorem 1

The function $h(B)$ is strictly concave, and the gradient is

$\displaystyle\nabla h(B)=2H(A,T)B.$ (8)

###### Proof

From (7), the directional derivative of $h$ at $B$ along $B^{\prime}$ is given
by ${\rm D}h(B)[B^{\prime}]=2{\rm tr}\left(B^{\prime\top}H(A,T)B\right)$.
Hence, the gradient of $h(B)$ is given by (8). Also, the Hessian of $h$ at any
$B\in{\mathbb{R}}^{n\times m}$ is given by ${\rm
Hess}\,h(B)[B^{\prime}]=2H(A,T)B^{\prime}$. Because $H(A,T)$ is a symmetric
negative definite matrix, $\langle B^{\prime},{\rm
Hess}\,h(B)[B^{\prime}]\rangle<0$ for any $B^{\prime}\in{\mathbb{R}}^{n\times
m}\backslash\\{0\\}$. Hence, $h(B)$ is strictly concave. $\Box$

The following corollary follows from (8).

###### Corollary 1

The gradient $\nabla h$ is $L(A,T)$-Lipschitz continuous, where

$\displaystyle L(A,T):=2\int_{0}^{T}\|\exp(A(T-t))\|_{F}^{2}dt.$ (9)

That is,

$\displaystyle\|\nabla h(B_{1})-\nabla h(B_{2})\|_{F}\leq
L(A,T)\|B_{1}-B_{2}\|_{F},$ (10)

where $B_{1}$ and $B_{2}$ are any real $n\times m$ matrices.

Algorithm 1 is the proposed algorithm for solving Problems 1 and 2. In
practice, we terminate the iteration if $\|\nabla h(B_{k})\|_{F}$ is
sufficiently small. Note that we must choose a nonzero $B_{0}$ at step 1. This
is because it follows from (8) that $B_{0}=0$ implies $\nabla h(B_{0})=0$.

Algorithm 1 Projected gradient methods for Problem 1 and 2.

1: Set $B_{0}\in\mathbb{R}^{n\times m}\backslash\\{0\\}$ and $t>L(A,T)$, where
$L(A,T)$ is defined as (9).

2: for $k=0,1,2,\ldots$ do

3: $B_{k+1}\in{\rm P}\left(B_{k}-\frac{1}{t}\nabla h(B_{k})\right)$, where
${\rm P}$ is the projection onto $Z_{s}\cap W$ for Problem 1 and $Z_{s}\cap
W_{[0,1]}$ for Problem 2.

4: end for

Using Corollary 1, we can obtain the following theorem regarding global
convergence and convergence rate. To show this and for the following
subsections, we define $\Lambda$ as the index set of ${\mathbb{R}}^{n\times
m}$. That is,
$\Lambda:=\\{(i,j)\,|\,i\in\\{1,2,\ldots,n\\},\,\,j\in\\{1,2,\ldots,m\\}\\}$.

###### Theorem 2

Any sequence $\\{B_{k}\\}$ generated by Algorithm 1 for Problem 1 (Problem 2)
globally converges to a stationary point of $\mathcal{I}_{Z_{s}\cap
W}(B)+h(B)$ ($\mathcal{I}_{Z_{s}\cap W_{[0,1]}}(B)+h(B)$) with locally linear
convergence rate.

###### Proof

We only prove the claim on Problem 1, because the proof on Problem 2 is
similar.

First, we show the global convergence. The objective function of (5), that is,
$\mathcal{I}_{Z_{s}\cap W}(B)+h(B)$, is a proper lower semicontinuous KL
function [30]. Moreover, $\mathcal{I}_{Z_{s}\cap W}(B)+h(B)$ is bounded below.
In fact, $B\in Z_{s}\cap W$ implies $\mathcal{I}_{Z_{s}\cap
W}(B)+h(B)>-\infty$, and $B\not\in Z_{s}\cap W$ yields $\mathcal{I}_{Z_{s}\cap
W}(B)=\infty>-h(B)$. Thus, if $B\not\in Z_{s}\cap W$, then
$\mathcal{I}_{Z_{s}\cap W}(B)+h(B)>-h(B)+h(B)=0$. From Corollary 1, $\nabla h$
is $L(A,T)$-Lipschitz continuous. Thus, Theorem 5.1 and Remark 5.2 in [31]
imply global convergence.

Next, we show locally linear convergence rate. Let $\Gamma_{k}$ be the set of
all sets of $k$ different indices $(i,j)\in\Lambda$. Using $\Gamma_{k}$ and
(7), Problem 1 can be rewritten as

$\displaystyle\min_{B}\quad\min_{\Omega\in\Gamma_{nm-s}}\,\,{\rm
tr}(BB^{\top}H(A,T))$ $\displaystyle{\rm subject\,to}\quad B\in H_{\Omega}\cap
W.$

with $H_{\Omega}:=\\{B\in{\mathbb{R}}^{n\times m}\,|\,B_{ij}=0\,\,{\rm
for}\,\,(i,j)\in\Omega\\}$. Moreover, this problem is equivalent to
$\min_{B}\,\,F(B)$, where
$F(B):=\min_{\Omega\in\Gamma_{nm-s}}\mathcal{I}_{H_{\Omega}\cap W}(B)+{\rm
tr}(BB^{\top}H(A,T))$. Because $F(B)$ is continuous on
$\\{B\in{\mathbb{R}}^{n\times m}\,|\,\partial F(B)\neq\emptyset\\}$, Corollary
5.2 in [24] implies locally linear convergence rate. $\Box$

In the following subsections, we show that the projections onto $Z_{s}\cap W$
and $Z_{s}\cap W_{[0,1]}$ in Algorithm 1 can be easily calculated. To this
end, we define

$\displaystyle\|X\|_{\Gamma}:=\sqrt{\sum_{(i,j)\in\Gamma}X_{ij}^{2}},\quad\langle
X,Y\rangle_{\Gamma}:=\sum_{(i,j)\in\Gamma}X_{ij}Y_{ij}$

for any $X,Y\in{\mathbb{R}}^{n\times m}$ and any $\Gamma\subset\Lambda$.

###  III-A Projection onto $Z_{s}\cap W$

The projections of $B\in{\mathbb{R}}^{n\times m}$ onto $Z_{s}$, $W$, and
$Z_{s}\cap W$ are defined by

$\displaystyle{\rm P}_{Z_{s}}(B)$ $\displaystyle:=\mathop{\rm
argmin}\limits_{X\in{\mathbb{R}}^{n\times
m}}\left\\{\|X-B\|_{F}^{2}\,|\,\|X\|_{0}\leq s\right\\},$ $\displaystyle{\rm
P}_{W}(B)$ $\displaystyle:=\mathop{\rm
argmin}\limits_{X\in{\mathbb{R}}^{n\times m}}\left\\{\|X-B\|_{F}^{2}\,|\,-{\bf
E}\leq X\leq{\bf E}\right\\},$ $\displaystyle=\min(\max(-{\bf E},B),{\bf E})$
$\displaystyle{\rm P}_{Z_{s}\cap W}(B)$ $\displaystyle:=\mathop{\rm
argmin}\limits_{X\in{\mathbb{R}}^{n\times m}}\left\\{\|X-B\|_{F}^{2}\,|\,-{\bf
E}\leq X\leq{\bf E},\right.$ $\displaystyle\quad\,\left.\,\,\|X\|_{0}\leq
s\right\\},$ (11)

respectively. Note that ${\rm P}_{Z_{s}}$ can be a set, in general. In fact,
$({\rm P}_{Z_{s}}(B))_{ij}=B_{ij}$ if $B_{ij}$ is contained in the $s$ first
largest entries in absolute value of $B$, and $({\rm P}_{Z_{s}}(B))_{ij}=0$
otherwise. That is, the set-valuedness of ${\rm P}_{Z_{s}}$ arises from the
fact that the $s$ largest entries may not be uniquely defined. Also, for any
$B\in{\mathbb{R}}^{n\times m}$, let $\Lambda_{s}(B)$ be any index set that
indicates the $s$ first largest entries in $|B_{ij}|$, $(i,j)\in\Lambda$. By
definition, for any $B,X\in{\mathbb{R}}^{n\times m}$,

$\displaystyle\|X\|_{F}^{2}=\|X\|^{2}_{\Lambda_{s}(B)}+\|X\|_{\Lambda\backslash\Lambda_{s}(B)}^{2}.$
(12)

Furthermore, we have the following lemma.

###### Lemma 1

For any $B,X\in{\mathbb{R}}^{n\times m}$,

$\displaystyle\|X-B\|_{\Lambda_{s}(B)}^{2}+\|X\|_{\Lambda\backslash\Lambda_{s}(B)}^{2}=\|X-{\rm
P}_{Z_{s}}(B)\|_{F}^{2}.$ (13)

###### Proof

It follows from (12) that

$\displaystyle\|X-{\rm P}_{Z_{s}}(B)\|_{F}^{2}=$ $\displaystyle\|X-{\rm
P}_{Z_{s}}(B)\|^{2}_{\Lambda_{s}(B)}$ $\displaystyle+\|X-{\rm
P}_{Z_{s}}(B)\|_{\Lambda\backslash\Lambda_{s}(B)}^{2}.$

By definition,

$\displaystyle||X-{\rm P}_{Z_{s}}(B)||^{2}_{\Lambda_{s}(B)}$
$\displaystyle=||X-B||^{2}_{\Lambda_{s}(B)}$ $\displaystyle\|X-{\rm
P}_{Z_{s}}(B)\|_{\Lambda\backslash\Lambda_{s}(B)}^{2}$
$\displaystyle=\|X\|_{\Lambda\backslash\Lambda_{s}(B)}^{2}.$

Thus, (13) holds. $\Box$

Using (12) and Lemma 1, we can show that the projection onto $Z_{s}\cap W$ is
the composition of ${\rm P}_{Z_{s}}$ and ${\rm P}_{W}$.

###### Theorem 3

$\displaystyle{\rm P}_{Z_{s}\cap W}={\rm P}_{W}\circ{\rm P}_{Z_{s}}.$

###### Proof

It follows from (11) and (12) that

$\displaystyle{\rm P}_{Z_{s}\cap W}(B)$ $\displaystyle=\mathop{\rm
argmin}\limits_{X\in{\mathbb{R}}^{n\times
m}}\\{\|X\|^{2}_{\Lambda_{s}(B)}+\|X\|_{\Lambda\backslash\Lambda_{s}(B)}^{2}$
$\displaystyle\quad-2\langle B,X\rangle_{\Lambda_{s}(B)}-2\langle
B,X\rangle_{\Lambda\backslash\Lambda_{s}(B)}+\|B\|_{F}^{2}\,|$
$\displaystyle\quad-{\bf E}\leq X\leq{\bf E},\,\,\|X\|_{0}\leq s\\}$ (14)

Because each $|B_{ij}|$, $(i,j)\in\Lambda_{s}(B)$ is greater than all
$|B_{ij}|$, $(i,j)\in\Lambda\backslash\Lambda_{s}(B)$ and $\|B\|_{F}^{2}$ is
constant, (14) is equivalent to

$\displaystyle{\rm P}_{Z_{s}\cap W}(B)$ $\displaystyle=\mathop{\rm
argmin}\limits_{X\in{\mathbb{R}}^{n\times
m}}\\{\|X\|^{2}_{\Lambda_{s}(B)}-2\langle B,X\rangle_{\Lambda_{s}(B)}\,|$
$\displaystyle\quad\|X\|_{\Lambda\backslash\Lambda_{s}(B)}^{2}=0,-{\bf E}\leq
X\leq{\bf E},\|X\|_{0}\leq s\\}$ $\displaystyle=\mathop{\rm
argmin}\limits_{X\in{\mathbb{R}}^{n\times
m}}\\{\|X\|^{2}_{\Lambda_{s}(B)}-2\langle B,X\rangle_{\Lambda_{s}(B)}\,|$
$\displaystyle\quad\|X\|_{\Lambda\backslash\Lambda_{s}(B)}^{2}=0,-{\bf E}\leq
X\leq{\bf E}\\}$ $\displaystyle=\mathop{\rm
argmin}\limits_{X\in{\mathbb{R}}^{n\times
m}}\\{\|X-B\|^{2}_{\Lambda_{s}(B)}+\|X\|_{\Lambda\backslash\Lambda_{s}(B)}^{2}\,|$
$\displaystyle\quad-{\bf E}\leq X\leq{\bf E}\\}$ (15)

From (13) in Lemma 1, (15) is equivalent to

$\displaystyle{\rm P}_{Z_{s}\cap W}(B)$ $\displaystyle=\mathop{\rm
argmin}\limits_{X\in{\mathbb{R}}^{n\times m}}\left\\{\|X-{\rm
P}_{Z_{s}}(B)\|_{F}^{2}\,|\,-{\bf E}\leq X\leq{\bf E}\right\\}$
$\displaystyle={\rm P}_{W}({\rm P}_{Z_{s}}(B)).$

This completes the proof. $\Box$

In general, ${\rm P}_{Z_{s}\cap W}\neq{\rm P}_{Z_{s}}\circ{\rm P}_{W}$. That
is, the order of the projections ${\rm P}_{Z_{s}}$ and ${\rm P}_{W}$ is not
commutative. In fact, for example, suppose that $B=\begin{pmatrix}3\\\
-4\end{pmatrix}$ and $s=1$. Then, $P_{W}(P_{Z_{1}}(B))=\begin{pmatrix}0\\\
-1\end{pmatrix}$ and $P_{Z_{1}}(P_{W}(B))=\begin{pmatrix}1\\\ 0\end{pmatrix}$.
Thus,
$\|P_{W}(P_{Z_{1}}(B))-B\|_{F}^{2}=18<20=\|P_{Z_{1}}(P_{W}(B))-B\|_{F}^{2}$.

###  III-B Projection onto $Z_{s}\cap W_{[0,1]}$

The projections of $B\in{\mathbb{R}}^{n\times m}$ onto
$\\{B\in{\mathbb{R}}^{n\times m}\,|\,B\geq 0\\}$,
$\\{B\in{\mathbb{R}}^{n\times m}\,|\,B\leq{\bf E}\\}$, and $Z_{s}\cap
W_{[0,1]}$ are defined by

$\displaystyle{\rm P}_{\geq 0}(B)$ $\displaystyle:=\mathop{\rm
argmin}\limits_{X\in{\mathbb{R}}^{n\times m}}\left\\{\|X-B\|_{F}^{2}\,|\,X\geq
0\right\\},$ $\displaystyle=\max(0,B),$ $\displaystyle{\rm P}_{\leq 1}(B)$
$\displaystyle:=\mathop{\rm argmin}\limits_{X\in{\mathbb{R}}^{n\times
m}}\left\\{\|X-B\|_{F}^{2}\,|\,X\leq{\bf E}\right\\},$
$\displaystyle=\min({\bf E},B),$ $\displaystyle{\rm P}_{Z_{s}\cap
W_{[0,1]}}(B)$ $\displaystyle:=\mathop{\rm
argmin}\limits_{X\in{\mathbb{R}}^{n\times m}}\left\\{\|X-B\|_{F}^{2}\,|\,0\leq
X\leq{\bf E},\right.$ $\displaystyle\quad\,\left.\,\,\|X\|_{0}\leq
s\right\\},$ (16)

respectively. For any $B\in{\mathbb{R}}^{n\times m}$, the index sets
$\Gamma_{<0}(B)$, $\Gamma_{[0,1]}(B)$, and $\Gamma_{>1}(B)$ are defined by

$\displaystyle\Gamma_{<0}(B)$ $\displaystyle:=\\{(i,j)\in\Lambda\,|\,B<0\\},$
$\displaystyle\Gamma_{[0,1]}(B)$ $\displaystyle:=\\{(i,j)\in\Lambda\,|\,0\leq
B\leq{\bf E}\\},$ $\displaystyle\Gamma_{>1}(B)$
$\displaystyle:=\\{(i,j)\in\Lambda\,|\,B>{\bf E}\\},$

respectively. By definition, for any $B,X\in{\mathbb{R}}^{n\times m}$,

$\displaystyle\|X\|_{F}^{2}=\|X\|^{2}_{\Gamma_{<0}(B)}+\|X\|_{\Gamma_{[0,1]}(B)}+\|X\|_{\Gamma_{>1}(B)}^{2}.$
(17)

Also, for any $B\in{\mathbb{R}}^{n\times m}$, let $\tilde{\Lambda}_{s}(B)$ be
the index set that indicates the $\min(s,|\tilde{\Lambda}(B)|)$ first largest
entries in $|B_{ij}|$, $(i,j)\in\tilde{\Lambda}(B)$, where

$\displaystyle\tilde{\Lambda}(B):=\Gamma_{[0,1]}(B)\cup\Gamma_{>1}(B).$

The projection onto $Z_{s}\cap W_{[0,1]}$ is the composition of ${\rm P}_{\geq
0}$, ${\rm P}_{Z_{s}}$, and ${\rm P}_{\leq 1}$ as follows.

###### Theorem 4

$\displaystyle{\rm P}_{Z_{s}\cap W_{[0,1]}}={\rm P}_{\leq 1}\circ{\rm
P}_{Z_{s}}\circ{\rm P}_{\geq 0}.$

###### Proof

It follows from (16) and (17) that

$\displaystyle{\rm P}_{Z_{s}\cap W_{[0,1]}}(B)$ $\displaystyle=\mathop{\rm
argmin}\limits_{X\in{\mathbb{R}}^{n\times
m}}\\{\|X-B\|^{2}_{\Gamma_{<0}(B)}+\|X-B\|_{\Gamma_{[0,1]}(B)}^{2}$
$\displaystyle\quad+\|X-B\|_{\Gamma_{>1}(B)}^{2}\,|\,0\leq X\leq{\bf
E},\,\|X\|_{0}\leq s\\}.$

If $(i,j)\in\Gamma_{<0}(B)$, then $B_{ij}<0$. Thus,
$\|X-B\|_{\Gamma_{<0}(B)}^{2}=\|X\|_{\Gamma_{<0}(B)}^{2}-2\langle
B,X\rangle_{\Gamma_{<0}(B)}+\|B\|_{\Gamma_{<0}(B)}^{2}$ subject to $0\leq
X\leq{\bf E}$, $\|X\|_{0}\leq s$ is minimized when $\|X\|_{\Gamma_{<0}(B)}=0$.
Thus,

$\displaystyle{\rm P}_{Z_{s}\cap W_{[0,1]}}(B)$ $\displaystyle=\mathop{\rm
argmin}\limits_{X\in{\mathbb{R}}^{n\times
m}}\\{\|X-B\|^{2}_{\tilde{\Lambda}(B)}\,|\,\|X\|_{\Gamma_{<0}(B)}=0,$
$\displaystyle\quad 0\leq X\leq{\bf E},\,\|X\|_{0}\leq s\\}.$ (18)

Because $B_{ij}\geq 0$ for any $(i,j)\in\tilde{\Lambda}(B)$, (18) implies

$\displaystyle{\rm P}_{Z_{s}\cap W_{[0,1]}}(B)$ $\displaystyle=\mathop{\rm
argmin}\limits_{X\in{\mathbb{R}}^{n\times
m}}\\{\|X-B\|^{2}_{\tilde{\Lambda}(B)}\,|\,\|X\|_{\Gamma_{<0}(B)}=0,$
$\displaystyle\quad X\leq{\bf E},\,\|X\|_{0}\leq s\\}.$ (19)

Moreover, because

$\displaystyle\|X-B\|_{\tilde{\Lambda}(B)}^{2}=$
$\displaystyle\|X\|_{\tilde{\Lambda}_{s}(B)}^{2}+\|X\|_{\tilde{\Lambda}(B)\backslash\tilde{\Lambda}_{s}(B)}^{2}-2\langle
B,X\rangle_{\tilde{\Lambda}_{s}(B)}$ $\displaystyle-2\langle
B,X\rangle_{\tilde{\Lambda}(B)\backslash\tilde{\Lambda}_{s}(B)}+\|B\|_{\tilde{\Lambda}(B)}^{2}$

and each $B_{ij}$ for $(i,j)\in\tilde{\Lambda}_{s}(B)$ is greater than all
$B_{ij}$, $(i,j)\in\tilde{\Lambda}(B)\backslash\tilde{\Lambda}_{s}(B)$, (19)
yields

$\displaystyle{\rm P}_{Z_{s}\cap W_{[0,1]}}(B)$ $\displaystyle=\mathop{\rm
argmin}\limits_{X\in{\mathbb{R}}^{n\times
m}}\\{\|X-B\|^{2}_{\tilde{\Lambda}_{s}(B)}+\|X\|_{\Gamma_{<0}(B)}^{2}$
$\displaystyle\quad+\|X\|_{\tilde{\Lambda}(B)\backslash\tilde{\Lambda}_{s}(B)}^{2}\,|\,X\leq{\bf
E}\\}$ $\displaystyle=\mathop{\rm argmin}\limits_{X\in{\mathbb{R}}^{n\times
m}}\\{\|X-{\rm P}_{Z_{s}}({\rm P}_{\geq
0}(B))\|^{2}_{\tilde{\Lambda}_{s}(B))}$ $\displaystyle\quad+\|X-{\rm
P}_{Z_{s}}({\rm P}_{\geq 0}(B))\|_{\Gamma_{<0}(B)}^{2}$
$\displaystyle\quad+\|X-{\rm P}_{Z_{s}}({\rm P}_{\geq
0}(B))\|_{\tilde{\Lambda}(B)\backslash\tilde{\Lambda}_{s}(B)}^{2}\,|\
X\leq{\bf E}\\}$ $\displaystyle=\mathop{\rm
argmin}\limits_{X\in{\mathbb{R}}^{n\times m}}\\{\|X-{\rm P}_{Z_{s}}({\rm
P}_{\geq 0}(B))\|^{2}_{F}\,|\,X\leq{\bf E}\\}$ $\displaystyle={\rm P}_{\leq
1}({\rm P}_{Z_{s}}({\rm P}_{\geq 0}(B))).$

Here, the first equality follows from a similar discussion to the proof of
Theorem 2, and the second equality follows from

$\displaystyle\|X-{\rm P}_{Z_{s}}({\rm P}_{\geq
0}(B))\|_{\tilde{\Lambda}_{s}(B)}$
$\displaystyle=\|X-B\|_{\tilde{\Lambda}_{s}(B)},$ $\displaystyle||X-{\rm
P}_{Z_{s}}({\rm P}_{\geq 0}(B))||_{\Gamma_{<0}(B)}$
$\displaystyle=||X||_{\Gamma_{<0}(B)},$ $\displaystyle\|X-{\rm P}_{Z_{s}}({\rm
P}_{\geq 0}(B))\|_{\tilde{\Lambda}(B)\backslash\tilde{\Lambda}_{s}(B)}$
$\displaystyle=\|X\|_{\tilde{\Lambda}(B)\backslash\tilde{\Lambda}_{s}(B)}.$

This completes the proof. $\Box$

Note that, similarly to in Theorem 3, we cannot change the order of ${\rm
P}_{\leq 1}$, ${\rm P}_{Z_{s}}$, and ${\rm P}_{\geq 0}$. In fact, for example,
suppose that $B=\begin{pmatrix}3\\\ -4\end{pmatrix}$ and $s=1$. Then, $P_{\leq
1}(P_{Z_{1}}(P_{\geq 0}(B)))=\begin{pmatrix}1\\\ 0\end{pmatrix}$ and $P_{\geq
0}(P_{Z_{1}}(P_{\leq 1}(B)))=\begin{pmatrix}0\\\ 0\end{pmatrix}$. Thus,
$\|P_{\leq 1}(P_{Z_{1}}(P_{\geq 0}(B)))-B\|_{F}^{2}=20<25=\|P_{\geq
0}(P_{Z_{1}}(P_{\leq 1}(B)))-B\|_{F}^{2}$.

##  IV Numerical Experiments

This section describes the results of the numerical experiment performed using
Algorithm 1. In all the cases, we set $m=1$, although Algorithm 1 can also be
used for $m>1$. That is, we only considered a single input case because the
case is already sufficiently difficult. In fact, it has been known that the
problem of finding $B\in{\mathbb{R}}^{n\times 1}$ such that system (1) is
controllable is NP hard [7]. Moreover, we set final time $T=10$ and
$t=1.1L(A,T)$ in Algorithm 1.

For Problem 1, we constructed $A$ in system (1) by using MATLAB command ${\rm
sprandn}$ (that is a sparse normally distributed random matrix generator) and
the Watts–Strogats model with $n$ nodes, $6$ average degree, and $0.05$
rewiring probability [32]. For Problem 2, we constructed a Metzler matrix $A$
in system (1) by using MATLAB command ${\rm sprand}$ (that is a sparse
uniformly distributed random matrix generator on the interval $(0,1)$) and the
Watts–Strogats model that has the same parameters as those used in Problem 1.

Figs. 1 and 2 show the relations between controllability index $-h(B)$ and
sparsity parameter $s$ in Problems 1 and 2, respectively. Here, initial point
$B_{0}$ in Algorithm 1 was the same for all $s$ when $n$ was fixed. As shown
in Fig. 1, $-h(B)$ tended to increase as $s$ increased. However, $-h(B)$ did
not monotonically increase, although $-h(B)$ is higher as $s$ increases when
we use global optimal solutions, as mentioned in Remark 2. This means that
local optimal solutions to Problem 1 could be obtained using Algorithm 1. In
contrast, $-h(B)$ monotonically increased as $s$ increased for Problem 2.

|
---|---
Figure 1: The relations between controllability index $-h(B)$ and sparsity parameter $s$ in Problem 1. The top-left and top-right figures show the simulation results in $n=50$ and $n=100$, respectively. The bottom-left and bottom-right figures show the simulation results in $n=200$ and $n=300$, respectively. The solid and dashed lines indicate the results of the sprandn and Watts–Strogats cases, respectively. |
---|---
Figure 2: The relations between controllability index $-h(B)$ and sparsity
parameter $s$ in Problem 2. The top-left and top-right figures show the
simulation results in $n=50$ and $n=100$, respectively. The bottom-left and
bottom-right figures show the simulation results in $n=200$ and $n=300$,
respectively. The solid and dashed lines indicate the results of the sprand
and Watts–Strogats cases, respectively.

The following non-trivial results were obtained.

  1. 1.

For Problems 1 and 2, the controllability characteristic changed as $s$
increased. That is, the slope of $-h(B)$ varied for small and sufficiently
large values of $s$.

  2. 2.

For Problems 1 and 2, the controllability index $-h(B)$ did not saturate as
$s$ increased.

  3. 3.

In the case of $s=n$, $B$ generated by Algorithm 1 had $n$ non-zero elements
for Problems 1 and 2. That is, the controllability in terms of $-h(B)$ was
maximized when all the states $x_{i}(t)$, $i=1,2,\ldots,n$ were stimulated by
a common single input $u(t)\in{\mathbb{R}}$. Moreover, although the elements
of $B$ had values of only $1$ for Problem 2, as can be shown theoretically,
the corresponding elements in Problem 1 were vectors composed of non-trivial
combinations of $1$ and $-1$.

  4. 4.

For Problem 1 (Problem 2), although controllability in the sprandn (sprand)
cases were higher than that in the Watts–Strogats cases for $n=50$ and
$n=100$; however, the relations were reversed for $n=200$ and $n=300$.

###### Remark 3

We confirmed that ${\rm rank}\,\mathcal{C}_{T}(B)<10$ for Problems 1 and 2
when $n=50$, $100$, $200$, and $300$. That is, system (1) resulting from
Algorithm 1 was not controllable. Thus, $\lambda_{\min}(\mathcal{C}_{T}(B))=0$
and ${{\rm tr}(\mathcal{C}_{T}^{-1}(B))}$ could not be defined. However, we
could increase $-h(B)$. This means that the system controllability can be
increased on a low dimensional subspace, and it is expected that the such a
subspace is determined by the structure of $A$.

##  V Conclusion

We formulated two novel controllability maximization problems and developed a
simple projected gradient method for solving the problems. We proved that a
sequence generated by our method has global convergence with locally linear
convergence rate. Moreover, the projections used in the proposed method were
given explicitly. Numerical experiments demonstrated the effectiveness of our
method, and provided non-trivial results. In particular, it is indicated that
controllability characteristic changes as the parameter specifying sparsity
increases, and the change rate appears to be dependent on a network structure.
The analysis of the change rates for various network structures would be
considered in future work.

## Acknowledgment

This work was supported by Japan Society for the Promotion of Science KAKENHI
under Grant 19H04069.

## References

  * [1] S. Gu, F. Pasqualetti, M. Cieslak, Q. K. Telesford, A. B. Yu, A. E. Kahn, J. D. Medaglia, J. M. Vettel, M. B. Miller, S. T. Grafton, and D. S. Bassett, “Controllability of structural brain networks,” _Nature communications_ , vol. 6, pp. 1–10, 2015.
  * [2] T. M. Karrer, J. Z. Kim, J. Stiso, A. E. Kahn, F. Pasqualetti, U. Habel, and D. Bassett, “A practical guide to methodological considerations in the controllability of structural brain networks,” _Journal of Neural Engineering_ , 2020.
  * [3] J. Z. Kim, J. M. Soffer, A. E. Kahn, J. M. Vettel, F. Pasqualetti, and D. S. Bassett, “Role of graph architecture in controlling dynamical networks with applications to neural systems,” _Nature physics_ , vol. 14, no. 1, pp. 1–8, 2018.
  * [4] A. Li, S. P. Cornelius, Y.-Y. Liu, L. Wang, and A.-L. Barabási, “The fundamental advantages of temporal networks,” _Science_ , vol. 358, no. 6366, pp. 1042–1046, 2017.
  * [5] Y.-Y. Liu, J.-J. Slotine, and A.-L. Barabási, “Controllability of complex networks,” _Nature_ , vol. 473, no. 7346, p. 167, 2011.
  * [6] T. Nepusz and T. Vicsek, “Controlling edge dynamics in complex networks,” _Nature Physics_ , vol. 8, no. 7, pp. 568–573, 2012.
  * [7] A. Olshevsky, “Minimal controllability problems,” _IEEE Transactions on Control of Network Systems_ , vol. 1, no. 3, pp. 249–258, 2014.
  * [8] G. Yan, P. E. Vértes, E. K. Towlson, Y. L. Chew, D. S. Walker, W. R. Schafer, and A.-L. Barabási, “Network control principles predict neuron function in the Caenorhabditis elegans connectome,” _Nature_ , vol. 550, no. 7677, pp. 519–523, 2017.
  * [9] Z. Yuan, C. Zhao, Z. Di, W.-X. Wang, and Y.-C. Lai, “Exact controllability of complex networks,” _Nature communications_ , vol. 4, no. 1, pp. 1–9, 2013\.
  * [10] Y.-Y. Liu and A.-L. Barabási, “Control principles of complex systems,” _Reviews of Modern Physics_ , vol. 88, no. 3, pp. 1–58, 2016.
  * [11] E. Tang and D. S. Bassett, “Colloquium: Control of dynamics in brain networks,” _Reviews of modern physics_ , vol. 90, no. 3, p. 031003, 2018\.
  * [12] A. J. Whalen, S. N. Brennan, T. D. Sauer, and S. J. Schiff, “Observability and controllability of nonlinear networks: The role of symmetry,” _Physical Review X_ , vol. 5, no. 1, p. 011005, 2015.
  * [13] J. G. T. Zañudo, G. Yang, and R. Albert, “Structure-based control of complex networks with nonlinear dynamics,” _Proceedings of the National Academy of Sciences_ , vol. 114, no. 28, pp. 7234–7239, 2017.
  * [14] R. F. Galán, “On how network architecture determines the dominant patterns of spontaneous neural activity,” _PloS one_ , vol. 3, no. 5, 2008.
  * [15] A. Clark, B. Alomair, L. Bushnell, and R. Poovendran, “Submodularity in input node selection for networked linear systems: Efficient algorithms for performance and controllability,” _IEEE Control Systems Magazine_ , vol. 37, no. 6, pp. 52–74, 2017.
  * [16] L. Romao, K. Margellos, and A. Papachristodoulou, “Distributed actuator selection: achieving optimality via a primal-dual algorithm,” _IEEE control systems letters_ , vol. 2, no. 4, pp. 779–784, 2018.
  * [17] T. H. Summers, F. L. Cortesi, and J. Lygeros, “On submodularity and controllability in complex dynamical networks,” _IEEE Transactions on Control of Network Systems_ , vol. 3, no. 1, pp. 91–101, 2016.
  * [18] F. Pasqualetti, S. Zampieri, and F. Bullo, “Controllability metrics, limitations and algorithms for complex networks,” _IEEE Transactions on Control of Network Systems_ , vol. 1, no. 1, pp. 40–52, 2014.
  * [19] V. Tzoumas, M. A. Rahimian, G. J. Pappas, and A. Jadbabaie, “Minimal actuator placement with bounds on control effort,” _IEEE Transactions on Control of Network Systems_ , vol. 3, no. 1, pp. 67–78, 2016.
  * [20] T. Ikeda and K. Kashima, “Sparsity-constrained controllability maximization with application to time-varying control node selection,” _IEEE Control Systems Letters_ , vol. 2, no. 3, pp. 321–326, 2018.
  * [21] A. Olshevsky, “On a relaxation of time-varying actuator placement,” _arXiv:1912.09454_ , 2019.
  * [22] W. M. Haddad, V. Chellaboina, and Q. Hui, _Nonnegative and compartmental dynamical systems_. Princeton University Press, 2010.
  * [23] K. Sato and A. Takeda, “Construction methods of the nearest positive system,” _IEEE Control Systems Letters_ , vol. 4, no. 1, pp. 97–102, 2020.
  * [24] G. Li and T. K. Pong, “Calculus of the exponent of Kurdyka–Łojasiewicz inequality and its applications to linear convergence of first-order methods,” _Foundations of computational mathematics_ , vol. 18, no. 5, pp. 1199–1232, 2018.
  * [25] R. E. Kalman, “On the general theory of control systems,” in _Proceedings First International Conference on Automatic Control, Moscow, USSR_ , 1960, pp. 481–492.
  * [26] ——, “Mathematical description of linear dynamical systems,” _Journal of the Society for Industrial and Applied Mathematics, Series A: Control_ , vol. 1, no. 2, pp. 152–192, 1963.
  * [27] E. Wu-Yan, R. F. Betzel, E. Tang, S. Gu, F. Pasqualetti, and D. S. Bassett, “Benchmarking measures of network controllability on canonical graph models,” _Journal of Nonlinear Science_ , pp. 1–39, 2018.
  * [28] E. J. Candes and T. Tao, “Decoding by linear programming,” _IEEE Transactions on Information Theory_ , vol. 51, no. 12, pp. 4203–4215, 2005.
  * [29] R. Tibshirani, “Regression shrinkage and selection via the lasso,” _Journal of the Royal Statistical Society: Series B (Methodological)_ , vol. 58, no. 1, pp. 267–288, 1996.
  * [30] J. Bolte, S. Sabach, and M. Teboulle, “Proximal alternating linearized minimization or nonconvex and nonsmooth problems,” _Mathematical Programming_ , vol. 146, no. 1-2, pp. 459–494, 2014.
  * [31] H. Attouch, J. Bolte, and B. F. Svaiter, “Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward–backward splitting, and regularized Gauss–Seidel methods,” _Mathematical Programming_ , vol. 137, no. 1-2, pp. 91–129, 2013.
  * [32] D. J. Watts and S. H. Strogatz, “Collective dynamics of ‘small-world’ networks,” _Nature_ , vol. 393, no. 6684, pp. 440–442, 1998.


# Observation of universal behaviour of ultracold quantum critical gases

Hongwei Xiong  State Key Laboratory of Magnetic Resonance and Atomic and
Molecular Physics, Wuhan Institute of Physics and Mathematics, Chinese Academy
of Sciences, Wuhan 430071, P. R. China  Xinzhou Tan  State Key Laboratory of
Magnetic Resonance and Atomic and Molecular Physics, Wuhan Institute of
Physics and Mathematics, Chinese Academy of Sciences, Wuhan 430071, P. R.
China  Graduate School of the Chinese Academy of Sciences, P. R. China  Bing
Wang  State Key Laboratory of Magnetic Resonance and Atomic and Molecular
Physics, Wuhan Institute of Physics and Mathematics, Chinese Academy of
Sciences, Wuhan 430071, P. R. China  Graduate School of the Chinese Academy of
Sciences, P. R. China  Lijuan Cao  State Key Laboratory of Magnetic Resonance
and Atomic and Molecular Physics, Wuhan Institute of Physics and Mathematics,
Chinese Academy of Sciences, Wuhan 430071, P. R. China  Graduate School of the
Chinese Academy of Sciences, P. R. China  Baolong Lű  State Key Laboratory of
Magnetic Resonance and Atomic and Molecular Physics, Wuhan Institute of
Physics and Mathematics, Chinese Academy of Sciences, Wuhan 430071, P. R.
China

###### Abstract

Quantum critical matter has already been studied in many systems, including
cold atomic gases. We report the observation of a universal behaviour of
ultracold quantum critical Bose gases in a one-dimensional optical lattice. In
the quantum critical region above the Berezinskii-Kosterlitz-Thouless
transition, the relative phase fluctuations between neighboring subcondensates
and spatial phase fluctuations of quasi-2D subcondensates coexist. We study
the density probability distribution function when both these two phase
fluctuations are considered. A universal exponential density probability
distribution is demonstrated experimentally, which agrees well with a simple
theoretical model by considering these two phase fluctuations.

The nature of quantum criticality Hertz ; Sachdev ; Cole  driven by quantum
fluctuations is still a great puzzle, despite of the remarkable advances in
heavy-fermion metals and rare-earth-based intermetallic compounds, etc Gegen .
New understanding of quantum criticality is widely believed to be a key to
resolving open questions in metal-insulator transitions Imada , high
temperature superconductivity HS  and novel material design, etc. Cold atoms
in optical lattices provide a unique chance to not only simulate other
strongly correlated systems RMP-Bloch , but also study some models
unaccessible in solid state systems, particularly for the Bose-Hubbard model
Fisher ; Jaksch ; Bloch . Despite of its complexity, a strongly correlated
system in quantum critical regime is expected to exhibit a universal behaviour
described by a certain physical quantity.

Here we report the observation of universal behaviour for ultracold quantum
critical Bose gases in a one-dimensional optical lattice. Density probability
distributions of the released gases are measured for different depths of the
lattice potential. It was found that the density probability follows a simple
exponential law when the Bose gases reach the quantum critical region above
the Berezinskii-Kosterlitz-Thouless (BKT) transition BKTtheory1 ; BKTtheory2 ;
Pol ; BKTexp . This universal behaviour can be well understood in terms of our
theoretical model considering both the relative phase fluctuations of quasi-2D
subcondensates and spatial phase fluctuations of individual subcondensates
above the BKT transition. The method of density probability distribution
should provide a unique tool for identifying certain quantum phases of optical
lattice systems.

Ultracold Bose atoms in 1D, 2D and 3D optical lattices are widely studied by
the Bose-Hubbard model Fisher . At zero temperature, there is a continuous
quantum phase transition from superfluid (SF) to Mott insulator (MI). Because
of the strongly correlated quantum behaviour, the quantum critical point (QCP)
at absolute zero temperature distorts strongly the structure of the phase
diagram at finite temperatures, leading to the emergence of an unconventional
‘V-shaped’ quantum critical region Sachdev ; Ho ; Cap ; Kato ; Chin ; Haz
(Fig. 1). In analogy with a black hole, the crossover to quantum critical gas
involves crossing a ‘material event horizon’, which implies strongly that the
quantum critical gas has a simple and universal behaviour Cole .

Figure 1: Finite temperature phase diagram of ultracold Bose gases in a one-
dimensional optical lattice. $U$ and $J$ represent the strengths of the on-
site repulsion and of the nearest neighbor hopping in the Bose-Hubbard model,
respectively. By increasing the strength of an optical lattice, $U/J$ can be
increased by several orders of magnitude. The existence of the quantum
critical point (QCP) at zero temperature leads to a ‘V-shaped’ quantum
critical (QC) region. Apart from the ordinary SF-QC-MI transition, the BKT
transition for quasi-2D subcondensates can also occur at sufficiently low
temperatures. In the QC region above the BKT transition (Type-I QC region),
there exists a universal behaviour in the density probability distribution.
All the phase boundary lines, except for the dashed line between QC and MI
regions, are calculated with a total atom number of $1.1\times 10^{5}$ and
trap parameters used in the experiment. $T_{har}$ is the critical temperature
for the purely harmonic trap in absence of the optical lattice. The gray line
displays a phase transition path from SF to Type-I QC region. Data points of
(a4)-(a6) along the gray line are obtained from the corresponding absorption
images in Fig. 2

In the commonly used 3D optical lattice systems, the atoms in a lattice site
are strongly confined in all directions, with a mean occupation number per
lattice site of about $1\sim 3$ Bloch . Thus the spatial phase fluctuations
for the atoms in a single lattice site can be omitted, and the relevant
theoretical calculations based on the Bose-Hubbard model Fisher ; Jaksch  and
Wannier function Kohn  can give quantitative descriptions of almost all
experimental phenomenaRMP-Bloch . In contrast, for ultracold Bose atoms in a
1D optical lattice, unique characteristics may arise due to the following two
reasons:

(1) Quasi-2D Bose gas: In an experiment of 1D optical lattice system as ours,
there can be hundreds of atoms in a lattice site. The Bose gas in a lattice
site becomes quasi-two-dimensional (quasi-2D) Pedri ; Burger ; Had , if the
local trapping frequency $\widetilde{\omega}_{z}$ of a lattice site in the
lattice direction satisfies the condition
$\hbar\widetilde{\omega}_{z}>>k_{B}T$ Pethick .

(2) BKT transition: For such a quasi-2D Bose gas, the critical temperature in
the occupied lattice site is given by
$T_{2D}\approx\hbar\omega_{\perp}\left(N_{l}\zeta\left(2\right)\right)^{1/2}/k_{B}$,
with $N_{l}$ being the atomic number in a lattice site. Well below $T_{2D}$,
the whole system becomes a chain of subcondensates. If there is no correlation
between the subcondensates in different lattice sites, a quasi-2D gas
undergoes a BKT transition at $T_{BKT}=T_{2D}/4$ BKTtheory1 ; BKTtheory2 .
Beyond this critical value of $T_{BKT}$, it is favorable to create vortices in
quasi-2D subcondensates, and the unbinding of bound vortices will lead to
strong spatial phase fluctuations within each subcondensate.

Since the BKT transition is crucial to revealing the property of cold atoms in
a 1D optical lattice, its critical temperature is also plotted in the phase
diagram (red line in Fig. 1). The BKT transition line divides the QC region
into two parts (Type I and Type II). Our experiments focus on the transition
from SF region to QC region above the BKT transition.

Figure 2: Density probability distribution for different lattice depths.
(a1)-(a6) in the left column are absorption images showing the density
distribution of the released atomic clouds. The field of view is $0.9\times
0.9$ mm2, and the pixel size is $\Delta^{2}=9.0\times 9.0$ $\mathrm{\mu
m}^{2}$. $T_{c}$ denotes the critical temperature of the atomic gas in the
combined trap. For moderate $s$, the interference fringes are prominent.
Further increasing $s$, however, we see complete disappearance of the
interference fringes (a6). This is due to the crossover from SF to QC region
as $s$ is increased. (b1)-(b6) in the right column are the corresponding
density probability distributions. Open circles are the data points calculated
from the density distributions of the images in the left column. The blue
solid curves give the exponential density probability distribution $P_{0}^{e}$
defined in text. When $s$ becomes large enough, the data points agree well
with the exponential curve (see b5 and b6), which shows clearly a universal
behaviour in the QC region above the BKT transition.

The experiments started with pure 87Rb condensates confined in a magnetic trap
with axial and transverse trapping frequencies of
$\left\\{\omega_{\bot},\omega_{z}\right\\}=2\pi\left\\{83.7,7.6\right\\}$ Hz.
The 1D optical lattice was formed by a retroreflected laser beam of
$\lambda=800\,$nm along the axis ($z$ direction) of the condensate. This laser
beam was ramped up to a given intensity over a time of $50$ ms, yielding a
lattice potential $V_{opt}=sE_{R}\sin^{2}\left(2\pi z/\lambda\right)$, with
$E_{R}$ being the recoil energy of an atom absorbing one lattice photon. After
a holding time of $10$ ms, we suddenly switched off the combined potential and
allowed the cold atomic cloud to expand freely for a time of $30$ ms.

The expanded atomic gases were probed using the conventional absorption
imaging technique. Since the probe beam is applied along $x$ direction, what
an absorption image records is the two-dimensional ($y-z$) column density
profile of the atomic cloud. We denote by $N_{1}^{ph}\left(y,z\right)$ and
$N_{2}^{ph}\left(y,z\right)$ the number of photons detected in the pixel at
position $\left(y,z\right)$ with and without the atomic cloud, respectively.
Then the density distribution is written as
$n_{2D}\left(y,z\right)=\ln\left[N_{2}^{ph}\left(y,z\right)/N_{1}^{ph}\left(y,z\right)\right]\Delta^{2}/\sigma_{e}$,
with $\sigma_{e}$ being the absorption cross section of a single atom and
$\Delta$ the pixel size.

The left column of Fig. 2 ((a1)-(a6)) displays the density distributions for
different lattice depths $s$. The temperatures of the system were inferred
from the condensate fractions as described in the Appendix. From the top three
images ((a1)-(a3)), we see an increasing of the side peaks in the interference
patterns as $s$ is increased. Further increasing $s$, we find the interference
fringes become blurred and disappear eventually (see (a4)-(a6)), similar to
the experimental observation in Ref. Orzel . The gradual disappearance of the
interference fringes can be partially explained by the increasing of the
random relative phase Orzel  during the SF-QC transition. In Figs. 2(a5)-(a6),
except for the density fluctuations along $z$ direction associated with the
random relative phase between different subcondensates, we also see
significant density fluctuations along $y$ direction. The simultaneous
existence of density fluctuations in these two directions is further
demonstrated in Fig. 3 (c) and (d) for a typical case of $s=30.7$. These
experimental results suggest that in the QC region above the BKT transition,
both the random relative phases and spatial phase fluctuations play important
roles in the expanded density distribution.

Figure 3: Density fluctuations of the released atomic clouds. (a)-(d) are one-
dimensional cuts through the corresponding images. For $s=30.7$, the density
distribution is highly fluctuated along both the horizontal ($z$) and vertical
($y$) directions.

To qualitatively analyze the density fluctuations, we consider the following
density probability distribution

$P\left(n\right)=\frac{S\left(n-\delta n/2\right)-S\left(n+\delta
n/2\right)}{\delta n\cdot S_{\mathrm{total}}}\,,$ (1)

where $\delta n$ is the width of a density interval, $S_{\mathrm{total}}$ is
the total area of the region occupied by the atomic gas, while
$S\left(n\right)$ is the area of the region where the density larger than $n$.
The averaged density is $n_{s}=1/S_{\mathrm{total}}$. For convenience, we
define a dimensionless density as $n_{0}=n/n_{s}$. The dimensionless density
probability distribution is then
$P_{0}\left(n_{0}\right)=n_{s}P\left(n\right)$.

The right column of Fig. 2 displays $P_{0}\left(n_{0}\right)$ calculated from
the corresponding density distribution. In each image, the pixel region chosen
for the calculation of $P_{0}\left(n_{0}\right)$ is that occupied by the cold
atoms, as enclosed by the red dashed ellipse. In all our calculations of
$P_{0}\left(n_{0}\right)$, $\delta n_{0}$ is just the horizontal spacing of
discrete data points. Due to the optical noise, $P_{0}\left(n_{0}\right)$ can
be nonzero even for negative $n_{0}$. This can be understood from the atomic
density formula based on the absorption imaging signal. Assuming $n_{0b}$ is
the minimum negative density for nonzero $P_{0}\left(n_{0}\right)$, optical
noise concentrates in the region of $n_{0b}<n_{0}<\left|n_{0b}\right|$ in each
$P_{0}\left(n_{0}\right)$ plot. Thus, $n_{0}>\left|n_{0b}\right|$ gives the
effective region where the optical noise is negligible and
$P_{0}\left(n_{0}\right)$ reflects the true density probability distribution.

The blue solid lines in Fig. 2(b1)-(b6) represent the following exponential
density probability distribution

$P_{0}^{e}=e^{-\left|n_{0}\right|}.$ (2)

It is obvious that, with increased $s$ (and hence $U/J$),
$P_{0}\left(n_{0}\right)$ has a tendency to $P_{0}^{e}$. In Fig. 2(b6), we see
that $P_{0}\left(n_{0}\right)$ agrees well with $P_{0}^{e}$ in the whole
effective region. Despite of the extremely complex many-body state in the QC
region above the BKT transition, our results show a simple universal
behaviour.

We have numerically simulated our experiments to explain the exponential
density probability distribution. When a released gas has experienced a free
expansion over a time of $t$, the 2D density distribution is given by

$n_{2D}\left(y,z,t\right)=\int
dx\left|\sum_{k}\sqrt{N_{k}}e^{i\left(\phi_{k\perp}\left(x,y,t\right)+\phi_{k}^{s}\left(z,t\right)\right)}\varphi_{k\perp}\left(x,y,t\right)\varphi_{kz}\left(z,t\right)\right|^{2}+n_{2D}^{opt}\left(x,y\right).$
(3)

Here $\varphi_{k\perp}\left(x,y,t\right)\varphi_{kz}\left(z,t\right)$ is the
wave packet of the expanded subcondensate initially in the $k$th lattice site,
with $|\varphi_{k\perp}\left(x,y,t\right)|^{2}$ being the Thomas-Fermi density
distribution in transverse ($x$-$y$) directions, and
$|\varphi_{kz}\left(z,t\right)|^{2}$ the Gaussian density distribution along
$z$ direction. Two classes of phase fluctuations enter the expression.
$\phi_{k\perp}\left(x,y,t\right)$ represents the spatial phase fluctuations in
transverse directions, while $\phi_{k}^{s}\left(z,t\right)$ comprises the
relative phase fluctuations between different subcondensates. $N_{k}$ denotes
the atomic number of the $k$th subcondensate. The term $n_{2D}^{opt}$ is added
to account for the optical noise. In numerical calculations, it was simulated
with the realistic optical noise distribution in the region having no atoms.

Figure 4: Theoretical simulation. Left column: simulated density distributions
based on Eq. (3). (a1), (a2) and (a3) are three typical cases in which an
atomic gas initially stays in SF, Type-I QC and Type-II QC regions,
respectively. In (a1), $s=11.2$, $\delta\phi_{k\perp}=0.2\pi$. In (a2),
$s=30.7$, both $\phi_{k}^{s}$ and $\phi_{k\perp}$ are completely random. In
(a3), $s=30.7$, $\phi_{k\perp}=0$, but $\phi_{k}^{s}$ is completely random.
Each of the color map has a size of $0.9\times 0.9$ mm2. Right column: Black
lines are the density probability distribution calculated from the
corresponding density distributions in the left column. Red lines are obtained
under the same conditions except that the optical noise is not included. The
exponential density probability distribution (blue lines) is also plotted for
comparison. The data in (b2) clearly shows an exponential density probability
distribution in the QC region above the BKT transition.

Figs. 4(a1) and (b1) give the simulated density distribution and density
probability distribution for an atomic gas initially in the combined trap with
$s=11.2$. Spatial phase fluctuations ($\phi_{k\perp}$) have been assumed to be
zero. Figs. 4(a2) and (b2) display similar calculations for $s=30.7$, but with
completely random phases $\phi_{k\perp}$ and $\phi_{k}^{s}$. The theoretical
results in Figs. 4(b2) match the exponential density probability distribution,
in agreement with our experimental results. The coexistence of two random
phases of $\phi_{k\perp}$ and $\phi_{k}^{s}$ makes the system similar to the
situation of the interference of randomly scattered waves Sheng , where
exponential distribution is also found. More general analytical derivations
for Eq. (3) and the exponential density probability distribution are given in
the Appendix. Figs. 4(a3) and (b3) give the simulation of $s=30.7$ for a cold
atomic gas in the QC region below the BKT transition. Accordingly,
$\phi_{k\perp}$ is assumed to be zero, while there exists a completely random
relative phase in $\phi_{k}^{s}$. No exponential density probability
distribution is found in this simulation. In addition, we do not notice
regular interference fringes either. Note that in Ref. Had , the QC region far
below the BKT transition was studied experimentally. However, because there
are only about $30$ subcondensates, high-contrast interference fringes were
still observed. Since our calculation is based on the experimental parameters
in this work, the total number of subcondensates is much higher (up to $250$).

To further explain the experimental phenomena, we stress that there are two
completely different situations in determining the BKT transition temperature
for cold atoms in 1D optical lattices.

(1) BKT transition in QC and MI regions: In these two regions, the quasi-2D
Bose gases in different lattice sites can be regarded as independent
subcondensates as far as the BKT transition is concerned. By minimizing the
free energy $F=E-TS$ ($E$ and $S$ are the energy and entropy of the Bose gas
in a lattice site when vortices are considered), one can get the BKT
transition temperature $T_{BKT}=T_{2D}/4$ Pethick . In the MI region above the
BKT transition, Eq. (3) can be used to simulate the density distribution by
including a completely random relative phase in $\phi_{k}^{s}$ Had , if the
averaged particle number per lattice site is much larger than $1$. Thus, in
the MI region above the BKT transition, one still expects the exponential
density probability distribution for a large number of more than one hundred
independent subcondensates.

(2) BKT transition in SF region: In the SF region, although the Bose gases in
the lattice sites become quasi-2D, they are highly correlated. For a series of
$M$ highly correlated quasi-2D Bose gases, the free energy becomes $F=ME-TS$.
If there are vortices in the quasi-2D subcondensates, vortices in different
lattice sites are highly correlated in spatial locations. Thus, $S$ can be
approximated as the entropy of a single subcondensate. Then the BKT transition
temperature becomes $T_{BKT}^{M}=MT_{BKT}$. Since $M>100$, we see that
$T_{BKT}^{M}>>T_{2D}$. This means that the temperatures of quasi-2D
condensates are always much lower than the BKT transition temperature. For the
situation of $s=11.2$ in our experiment, based on the estimation of the
particle number fluctuations $\delta N_{k}$ Jav  and the relation of
$\delta\phi_{k}^{s}\sim 2\pi/\delta N_{k}$ Orzel , we deduced a phase
fluctuation of $\delta\phi_{k}^{s}\approx 0.2\pi<<2\pi$, which shows strong
correlation between different lattice sites. As shown in Fig. 3(b), there are
no noticeable density fluctuations along $y$ direction, which confirms that
the initially confined gas is well below the BKT transition even though $T$ is
already beyond $T_{BKT}$.

We should also emphasize that, in the QC region above the BKT transition, the
system cannot be regarded as a classical thermal cloud although the
temperature is higher than the critical temperature $T_{c}$. A classical
thermal cloud features density fluctuations as $\delta^{2}n\approx n$. In
contrast, the exponential density probability distribution corresponds to much
larger density fluctuations of $\delta^{2}n\approx n^{2}$ XiongWu . This
exponential density probability distribution physically originates from the
superposition of quantum sates with a completely random phase, while the
universal behaviour lies in that the distribution is always exponential for
completely random superposition of quantum state XiongWu . The present work
supports the long-standing belief that the quantum states in the QC region are
strongly correlated, and there should be a simple universal behaviour, once
appropriate physical quantity is found Sachdev ; Cole ; Gegen .

###### Acknowledgements.

H. W. thank the cooperation with Biao Wu about the dynamical universal
behaviour of quantum chaos, which stimulates the present work. This work was
supported by National Key Basic Research and Development Program of China
under Grant No. 2006CB921406 and NSFC under Grant No. 10634060.

## References

  * (1) J. Hertz, Phys. Rev. B 14, 1165 (1976).
  * (2) S. Sachdev, Quantum Phase Transitions (Cambridge Univ. Press, Cambridge, 1999).
  * (3) P. Coleman, A. J. Schofield, Nature 433, 226 (2005).
  * (4) P. Gegenwart, Q. Si, F. Steglich, Nature Phys. 4, 186 (2008).
  * (5) M. Imada, A. Fujimori, Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).
  * (6) S. Sachdev, Rev. Mod. Phys. 75, 913 (2003).
  * (7) I. Bloch, J. Dalibard, W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * (8) M. P. A. Fisher, P. B. Weichman, G. Grinstein, D. S. Fisher, Phys. Rev. B 40, 546 (1989).
  * (9) D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).
  * (10) M. Greiner, O. Mandel, T. Esslinger, T. W. Ha̋nsch, I. Bloch, Nature 415, 39 (2002).
  * (11) V. L. Berezinskii, Sov. Phys. JETP 34, 610 (1972).
  * (12) J. M. Kosterlitz, D. J. Thouless, J. Phys. C 6, 1181 (1973).
  * (13) A. Polkovnikov, E. Altman, E. Demler, Proc. Natl Acad. Sci. USA 103, 6125 (2006).
  * (14) Z. Hadzibabic, P. Krüger, M. Cheneau, B. Battelier, J. Dalibard, Nature 441, 1118 (2006).
  * (15) R. B. Diener, Q. Zhou, H. Zhai, T. L. Ho, Phys. Rev. Lett. 98, 180404 (2007); Q. Zhou and T. L. Ho, Phys. Rev. Lett. 105, 245702 (2010).
  * (16) B. Capogrosso-Sansone, N. V. Prokof ’ev, B. V. Svistunov, Phys. Rev. B 75, 134302 (2007).
  * (17) Y. Kato, Q. Zhou, N. Kawashima, N. Trivedi, Nature Phys. 4, 617 (2008).
  * (18) N. Gemelke, X. B. Zhang, C. L. Hung, C. Chin, Nature 460, 995 (2009).
  * (19) K. R. A. Hazzard, E. J. Mueller, arXiv: 1006.0969 (2010).
  * (20) W. Kohn, Phys. Rev. 115, 809 (1959).
  * (21) P. Pedri et al., Phys. Rev. Lett. 87, 220401 (2001).
  * (22) S. Burger et al., Europhys. Lett. 57, 1 (2002).
  * (23) Z. Hadzibabic et al., Phys. Rev. Lett. 93, 180403 (2004).
  * (24) C. J. Pethick, H. Smith, Bose-Einstein condensation in dilute gases. (Cambridge University Press, Cambridge, 2008).
  * (25) C. Orzel, A. K. Tuchman, M. L. Fenselau, M. Yasuda, M. A. Kasevich, Science 291, 2386 (2001).
  * (26) P. Sheng, Introduction to Wave Scattering, Localization, and Mesoscopic Phenomena. Chap. 10 (Academic Press. Inc., California, 1995).
  * (27) J. Javanainen, M. Y. Ivanov, Phys. Rev. A 60, 2351 (1999).
  * (28) H. Xiong, B. Wu, arXiv: 1007.2771 (2010).

Appendix

Bose-Hubbard model. The combined potential of a one-dimensional (1D) optical
lattice and a harmonic trap is given by

$V=\frac{1}{2}m\omega_{\perp}^{2}\left(x^{2}+y^{2}\right)+\frac{1}{2}m\omega_{z}^{2}z^{2}+sE_{R}\sin^{2}\left(\frac{2\pi
z}{\lambda}\right).$ (4a) Here $E_{R}$ is the recoil energy of an atom
absorbing one lattice photon with wavelength $\lambda$. For an atom in a
lattice site, it experiences an effective harmonic potential along $z$
direction with angular frequency $\widetilde{\omega}_{z}\simeq
2\sqrt{s}E_{R}/\hbar$.

Bose-condensed gases in 1D, 2D and 3D optical lattices are widely studied by
the following Bose-Hubbard model FisherS

$\widehat{H}=-\sum\limits_{\left\langle
ij\right\rangle}J\widehat{a}_{i}^{{\dagger}}\widehat{a}_{j}+\frac{U}{2}\sum\limits_{i}\widehat{n}_{i}\left(\widehat{n}_{i}-1\right)+\sum\limits_{i}\left(\varepsilon_{i}-\mu\right)\widehat{n}_{i}.$
(4b)

Here $\widehat{n}_{i}=\widehat{a}_{i}^{{\dagger}}\widehat{a}_{i}$ is the
number operator at the $i$th site, with $\widehat{a}_{i}^{{\dagger}}$
($\widehat{a}_{i})$ being the boson creation (annihilation) operator. $U$ and
$J$ represent the strengths of the on-site repulsion and of the nearest
neighbor hopping, respectively. $\varepsilon_{i}$ describes an energy offset
due to the harmonic trap, and $\mu$ is the chemical potential. The phase
diagram given by Fig. 1 in the text is calculated with the above combined
potential and Bose-Hubbard model.

Universal exponential density probability distribution. The many-body quantum
state in the QC region can be written as in a general way

$\left|\Psi_{qcr}\right\rangle=\sum_{\left\\{\Sigma
N_{k}=N\right\\}}C\left(N_{1},\cdot\cdot\cdot,N_{2k_{M}+1}\right)\prod_{k}\frac{1}{\sqrt{N_{k}!}}\left(\widehat{a}_{k}^{{\dagger}}\right)^{N_{k}}\left|0\right\rangle.$
(4c)

Here $N_{k}$ denotes the atomic number in the $k$th lattice site. The total
atomic number is then $N=\Sigma_{k=-k_{M}}^{k=k_{M}}N_{k}$, with $2k_{M}+1$
being the total number of occupied lattice sites. For strongly correlated
quantum gases in the QC region, the function
$C\left(N_{1},\cdot\cdot\cdot,N_{2k_{M}+1}\right)$ is extremely complex, and
no analytical expression exists. After switching off the combined potential,
the density distribution is then

$n\left(x,y,z,t\right)=\left\langle\Psi_{qcr}\left(t\right)\right|\widehat{\Psi}^{{\dagger}}\widehat{\Psi}\left|\Psi_{qcr}\left(t\right)\right\rangle.$
(4d)

Assuming that the initial wave function in the $k$th lattice site is
$\varphi_{k}\left(x,y,z,t=0\right)=\varphi_{k\perp}\left(x,y,t=0\right)\varphi_{kz}\left(z,t=0\right)$,
for long-time evolution so that the expanded subcondensates have sufficient
overlapping with each other, we have

$n\left(x,y,z,t\right)\approx\left|\sum_{k}\sqrt{N_{k}}e^{i\phi_{k}^{s}\left(z,t\right)}\varphi_{k\perp}\left(x,y,t\right)\right|^{2}\left|\varphi_{0z}\left(z,t\right)\right|^{2}.$
(4e)

This formula holds in the QC region, because there are non-negligible particle
number fluctuations in each lattice site. Here
$\phi_{k}^{s}\left(z,t\right)=\phi_{k}^{r}+\phi_{k}^{0}$, with $\phi_{k}^{r}$
being a completely random relative phase in the QC region.
$\phi_{k}^{0}=m\left(z-k\lambda/2\right)^{2}/2\hbar t$ can be obtained
directly from the free expansion along $z$ direction. Without the random phase
$\phi_{k}^{r}$, the phase $\phi_{k}^{0}$ will lead to clear interference
fringes along $z$ direction with a period of $4\pi\hbar t/m\lambda$ PedriS .
It is understood that there would be no interference fringes along $z$
direction if $\phi_{k}^{s}$ is a completely random relative phase, as shown in
Figs. 2(a5)-(a6) in the text.

In the QC region above the BKT transition, there are two sorts of random
phases: (i) the completely random phase $\phi_{k}^{s}\left(z,t\right)$ in the
QC region; (ii) the spatially relevant random phase
$\phi_{k\perp}\left(x,y,t\right)$ above the BKT transition. In this region,
the density distribution can be further written as

$\displaystyle n\left(x,y,z,t\right)$ $\displaystyle\approx$
$\displaystyle\left|\sum_{k}\sqrt{N_{k}}e^{i\phi_{k\perp}\left(x,y,t\right)}e^{i\phi_{k}^{s}\left(z,t\right)}\right|^{2}$
(4f)
$\displaystyle\left|\varphi_{0\perp}\left(x,y,t\right)\varphi_{0z}\left(z,t\right)\right|^{2}.$

Here $\phi_{k\perp}\left(x,y,t\right)$ physically originates from the spatial
phase fluctuations of quasi-2D Bose gas above the BKT transition.

The density distribution recorded by a CCD is a 2D density distribution
$n_{2D}\left(y,z,t\right)=\int dxn\left(x,y,z,t\right)$. By assuming that
$\left|\varphi_{0\perp}\left(x,y,t\right)\varphi_{0z}\left(z,t\right)\right|^{2}$
is locally constant (i.e. the local density approximation), after integrating
out the $x$ variable, we have

$\displaystyle n_{2D}\left(y,z,t\right)$ $\displaystyle\approx$
$\displaystyle\left|\frac{1}{\sqrt{2k_{M}+1}}\sum_{k}\alpha_{k}\left(y,z,t\right)e^{i\beta_{k}\left(y,z,t\right)}\right|^{2}$
(4g) $\displaystyle n_{2D}^{0}\left(y,z,t\right),$

with $n_{2D}^{0}\left(y,z,t\right)=\int
dx\left|\sqrt{N}\varphi_{0\perp}\left(x,y,t\right)\varphi_{0z}\left(z,t\right)\right|^{2}$.
$\alpha_{k}\left(y,z,t\right)$ and $\beta_{k}\left(y,z,t\right)$ are
statistically independent random real function, because they originates from
two different random phases $\phi_{k}^{s}\left(z,t\right)$ and
$\phi_{k\perp}\left(x,y,t\right)$. $\beta_{k}\left(y,z,t\right)$ distributes
uniformly between $-\pi$ and $\pi$. The normalization condition requests
$\left\langle\alpha_{k}^{2}\right\rangle=1$.

Based on the central limit theorem Feller , we get the following exponential
density probability distribution

$P\left(n_{2D}\left(x,y,t\right)\right)=\frac{e^{-n_{2D}\left(x,y,t\right)/n_{2D}^{0}\left(x,y,t\right)}}{n_{2D}^{0}\left(x,y,t\right)}.$
(4h)

For uniform $n_{2D}^{0}$, it is straightforward to get Eq. (2) in the text.
For nonuniform $n_{2D}^{0}$, with the local density approximation, Eq. (2) in
the text is a very good approximation obtained from Eq. (4h).

In Figs. 2(b5)-(b6) in the text, we see that the density probability
distribution agrees well with the exponential density probability
distribution. However, because the system temperature is slightly smaller than
the critical temperature for Fig. 2(b5) ($s=22.3$), we expect Fig. 2(b6)
agrees better with the exponential density probability distribution. Here we
give the semilog plot of the density probability distribution for $s=22.3$ and
$s=30.7$, respectively. We do find that the situation of $s=30.7$ agrees
better with the exponential density probability distribution.

Figure 5: Semilog plot of the density probability distribution for $s=22.3$
and $s=30.7$. The blue line gives the exponential density probability
distribution.

System temperature. In our experiment, the total atomic number can be
accurately measured. The atomic numbers of Figs. 2(a4)-(a6) in the text are
$1.10\times 10^{5}$, $1.19\times 10^{5}$ and $1.13\times 10^{5}$,
respectively. To drive the system into the deeper region of superfluidity, we
have to further lower the temperature of an atomic sample before loading it to
the lattice. This is implemented by kicking more atoms out of the magnetic
trap during the evaporative cooling stage. Specifically, for Figs. 2(a2) and
(a3) in the text, the atomic numbers are reduced by about a factor of two
compared with Figs. 2(a4)-(a6) in the text.

However, the measurement of the system temperature is hampered by the lack of
accurate thermometry in optical lattices. Actually, the ratio $T/T_{c}$ is a
more important quantity than $T$ in characterizing the transition from
superfluid to quantum critical regions. Therefore, the inaccuracy of system
temperature does not affect the physics of the present work, as long as
$T/T_{c}$ can be somehow determined with a satisfying accuracy.

As predicted in Ref. Duan  and demonstrated in a recent experiment Phillips ,
the appearance of superfluid in an inhomogeneous lattice is associated with a
bimodal structure in the density profile. The condensate fraction $N_{0}/N$
can be obtained by fits to the bimodal distribution. As we have known,
$N_{0}/N$ shows a temperature dependence having a characteristic shape as

$N_{0}/N=1-(T/T_{c})^{\alpha}.$ (4i)

Apparently, $T/T_{c}$ can be directly determined if the parameter $\alpha$ is
known. Blakie _et al._ Tc-Blakie  had performed numerical calculations of the
critical temperature as well as the condensate fraction, for a combined trap
similar to that in our work. Their results show that $\alpha$ is close to $3$
at a moderate lattice depth, but it may drop slightly with increasing lattice
depth. Nevertheless, at a high depth level, the condensate fraction is usually
very small, and, as a result, the value of $T/T_{c}$ derived from Eq. (4i)
becomes insensitive to the potential errors of $\alpha$. Therefore, we just
set $\alpha=3$ in the calculation of $T/T_{c}$, which seems to be a reasonable
assumption.

Figure 6: A bimodal fit to the interference pattern in Fig. 2(a3) in the text.
The total atom number is measured as $N\simeq 5.9\times 10^{4}$, and the
atomic cloud was released from the combined trap with a lattice strength of
$s=11.2$. Black line, the measured atomic density which is obtained by
integrating over each column of pixels. Red line, a fitting curve including
six Gaussian peaks, where the three narrow peaks at the upper part are
attributed to superfluid. The condensate fraction is $N_{0}/N\approx 0.58$,
corresponding to $T/T_{c}\approx 0.75$.

For each absorption image, we first obtain the atomic density along the $z$
dimension by integrating pixels in each column. The narrow peaks riding over
broad ones are regarded as the condensate parts. We then fit both the broad
and the narrow peaks using Gaussian profiles, so as to get the condensate
fraction. Apparently, this method works only when $T<T_{c}$. In Fig. 2(a6) in
the text, the interference peaks vanish completely, indicating that $T>T_{c}$.
In this case, we are unable to extract the quantity $T/T_{c}$ any more with
this method. The universal exponential density probability distribution found
in this work also implies that the accurate measurement of system temperature
in the QC region above BKT transition is difficult.

## References

  * (1) M. P. A. Fisher, P. B. Weichman, G. Grinstein, D. S. Fisher, Phys. Rev. B 40, 546 (1989).
  * (2) P. Pedri et al., Phys. Rev. Lett. 87, 220401 (2001).
  * (3) G. D. Lin, W. Zhang, L. M. Duan, Phys. Rev. A 77, 043626 (2008).
  * (4) I. B. Spielman, W. D., Phillips, J. V. Porto, Phys. Rev. Lett. 100, 120402 (2008).
  * (5) P. B. Blakie, W. X. Wang, Phys. Rev. A 76, 053620 (2007).
  * (6) W. Feller, An Introduction to Probability Theory and its Applications Vol. 1 (John Wiley & Sons, New York, 1957).