# A New Perspective of Accelerated Gradient Methods: The Controlled Invariant
Manifold Approach

Revati Gunjal, Sushama Wagh, Syed Shadab Nayyer, Alex Stankovic, and Navdeep
M. Singh Revati Gunjal, Sushama Wagh, Syed Shadab Nayyer, and Navdeep M. Singh
are with Control and Decision Research Centre (CDRC), EED, Veermata Jijabai
Technological Institute, Mumbai 400019, India. Sushama Wagh is a visiting
scientist and Alex Stankovic is a distinguished scientist at SLAC National
Accelerator Laboratory, Menlo Park, CA, USA<EMAIL_ADDRESS>
###### Abstract

Gradient Descent (GD) is a ubiquitous algorithm for finding the optimal
solution to an optimization problem. For reduced computational complexity, the
optimal solution $\mathrm{x^{*}}$ of the optimization problem must be attained
in a minimum number of iterations. For this objective, the paper proposes a
genesis of an accelerated gradient algorithm through the controlled dynamical
system perspective. The objective of optimally reaching the optimal solution
$\mathrm{x^{*}}$ where $\mathrm{\nabla f(x^{*})=0}$ with a given initial
condition $\mathrm{x(0)}$ is achieved through control.

###### Index Terms:

Controlled Dynamical system, Gradient Descent, Manifold stabilization,
Nesterov’s Accelerated Gradient method, Passivity, and Immersion Approach

##  I Introduction

Optimization plays a central role in various areas such as machine learning,
data analysis, game theory, and a wide range of scientific and engineering
applications. The gradient descent (GD) method is a ubiquitous algorithm for
finding the optimal solution to an optimization problem through an iterative
scheme. In the literature, the GD algorithm is analyzed as a first-order
dynamical system, which helps one to get a better insight into the pros and
cons of the algorithm. The gradient-based strategies for optimization are
prevalent due to benefits like low iteration cost and easier deployment on
parallel and distributed processing architecture. However, these methods are
often inadequate to accommodate the growing complexity and scalability owing
to large-scale data analysis that demands faster convergence rates. These
issues motivated the requirement for accelerated gradient methods that achieve
convergence to the optimal solution with the minimum number of iterations. B.
T. Polyak in 1964 and Y. Nesterov in 1983 developed a class of accelerated
gradient algorithms with improved convergence rates and lower complexity as
compared to the GD method. The notion of acceleration has also been
incorporated with many optimization schemes such as stochastic optimization
[1] [2], non-convex optimization [3][4], composite optimization [5], and conic
programming [6], to name a few. These ideas have also been integrated into
non-Euclidean optimization and higher-order optimization algorithms.

The literature discusses various continuous-time dynamical system formulations
for accelerated gradient techniques, such as the Heavy Ball (HB) method [7],
Nesterov’s Accelerated Gradient (NAG) [8] method, and their extensions. The
above approaches analyze the GD algorithm through a second-order ordinary
differential equation (ODE). The momentum-based optimization algorithms,
namely, HB and its extension Heavy Ball with Friction (HBF), are formulated
based on the analogies of mechanical systems. The first-order GD system mimics
the movement of a drop of water as it slides on the profile defined by the
function $\mathrm{f(x)}$. The response of first-order GD is slow, and due to
the inappropriate selection of step size, the algorithm might get stuck at the
local minima. To tackle these limitations, B. T. Polyak introduced the HB
algorithm [7] that accelerates the GD scheme by adding the momentum term.
However, the response of the HB scheme exhibits an oscillatory behavior. To
improve the response of the HB method, a damping term is added, establishing
another strategy known as the HBF approach [9]. The ODE of HBF includes a
“geometric damping term” that involves the Hessian of the objective function.

Analysis of gradient-based techniques often incorporates the limiting ODE
perspective. However, these ODEs cannot differentiate between the two
fundamental algorithms, the HB method and Nesetrov’s Accelerated Gradient
method for Strongly Convex function (NAG-SC). In view of this, Bin Shi et. al.
[10] attempted an alternative limiting process to obtain the high-resolution
ODEs for these methods. As a result of the high-resolution ODEs a geometric
correction term appears. This geometric correction term includes Hessian of
the function that provides information about the local geometry and the
curvature of the descent. The knowledge of local geometry allows the selection
of larger step sizes and hence faster convergence without risking the
algorithm’s stability or loss of any information.

As evident from the above-mentioned methods, the accelerated response can be
achieved with fewer iterations by interpreting GD in terms of a second-order
ODE. However, these methods throw no light on the origins of these algorithms
and they still remain a conceptual mystery. Our work attempts to provide
insight into the origins of accelerated gradient algorithms and the approach
so taken may be helpful in the design of algorithms rather than just providing
an ODE-based interpretation of the existing accelerated gradient-based
methods. Heuristically, one could think of the algorithms as simply the
numerical integration of a continuous time dynamical system. And there is a
rich history of interpreting optimization algorithms as continuous time
dynamical systems [11] [12][13]. Further, L. Lassard et. al. [14] incorporated
the control theoretic perspective in the accelerated gradient methods for
optimization through the notion of Integral Quadratic Constraints (IQC). This
was utilized for the stability analysis of the optimization algorithms to
establish the upper bounds on the convergence rates of the gradient-based
algorithms. Furthermore, a continuous time variational, continuous-time
approach is proposed in [15] for understanding NAG method and its
generalizations. Motivated by the above approaches, this paper attempts to
present the controlled dynamical system perspective to optimization algorithms
and also explains the origin of the NAG in the continuous time domain.

The optimization problem is visualized as a manifold stabilization problem.
Manifold stabilization is a methodology developed in linear/nonlinear control
theory that aims to design a control input or feedback scheme that guarantees
the convergence (attractivity) of the system dynamics to a certain invariant
manifold. The manifold selection depends on the desired objectives and if one
considers the equilibrium point as the invariant manifold then the above
problem reduces to the equilibrium stabilization problem.

The manifold stabilization in this work is achieved through the Passivity and
Immersion (P&I) approach [16] [17]. This approach integrates the concepts of
immersion intended for immersing the system dynamics into the defined lower-
order target dynamics and the notion of attractivity of the manifold derived
from the passivity theory. P&I is proposed for control and stabilization, and
the concept is extended to manifold stabilization. Recently the approach has
been applied to parameter estimation [18] by visualizing the estimation
problem as a manifold stabilization problem, and the accurate estimation of
the unknown parameters is accomplished via a control perspective. P&I is an
efficient and structured formulation that is useful to solve any problem by
visualizing it as a manifold stabilization problem. An algorithm can be
devised systematically through P&I to obtain a controlled dynamical system
that fulfills a specific objective. This paper adapts the above approach to
the optimization framework.

The controlled invariant manifold is susceptible to numerical errors, which
can be viewed as perturbations, and as a result of which the invariant
manifold may not persist when subjected to such perturbations. A novel
approach is taken in this paper thereby the off-the-manifold dynamics
transversal to the manifold is added to ensure stronger contraction of the
ambient space and hence the persistence of the invariant manifold under
numerical perturbations. This aids in designing a numerically insensitive
algorithm for accelerated gradient descent.

The paper is structured as follows: The detailed procedure for the proposed
controlled invariant manifold approach for the accelerated gradient method is
elaborated in Section II. Section III provides insights into the applicability
of the proposed theory. Section IV provides a brief discussion about the
potential of the proposed philosophy. The paper is concluded along with the
future work in Section V. Mathematical formulations of existing methods and a
detailed procedure of the P&I approach are described in the Appendix.

##  II Controlled Invariant Manifold Approach for Accelerated Gradient Method

###  II-A Problem Definition

Consider an unconstrained optimization problem

$\begin{split}&\mathrm{min}\;\mathrm{f}(\mathrm{x})\\\
&\mathrm{s.t.}\;\mathrm{x}\in\mathbb{R}\end{split}$ (1)

###### Assumption 1

The objective function $\mathrm{f}$ is smooth $\mu$-strongly convex with
$L$-Lipschitz gradients.

###### Definition 1

(Strong Convexity) A twice differentiable function
$\mathrm{f:\mathbb{R}^{n}\rightarrow\mathbb{R}}$ is $\mu$-strongly convex [19]
with $\mu>0$ being the modulus of convexity, if its Hessian matrix
$\nabla^{2}\mathrm{f(x)}$ satisfies the matrix inequality

$\nabla^{2}\mathrm{f(x)\geq\mu I\hskip 56.9055pt\forall x\in\mathbb{R}^{n}}$
(2)

###### Definition 2

(L-Lipschitz Continuity) A differentiable function $\mathrm{f(x)}$ is smooth
iff it has a Lipschitz continuous gradient, i.e. iff $\exists L<\infty$ such
that

$\left\|\nabla\mathrm{f}(\mathrm{x}_{1})-\nabla\mathrm{f}(\mathrm{x}_{2})\right\|_{2}\leq
L\left\|\mathrm{x}_{1}-\mathrm{x}_{2}\right\|_{2}\>\>\>\>\forall\mathrm{x\in\mathbb{R}^{n}}$
(3)

The gradient descent (GD) for the above optimization problem can be written
as:

$\accentset{\mbox{\large.}}{\mathrm{x}}=-\nabla\mathrm{f(x)},$ (4)

The paper is concerned with a controlled dynamical system perspective to find
the optimal point $\mathrm{x}^{*}$ of $\mathrm{f(x)}$ via a Controlled
Invariant Manifold approach.

###  II-B Optimization as a Manifold Stabilisation Problem

The concept of manifold stabilization has been incorporated into the
optimization problem. In the context of control theory, a manifold is defined
as a lower-dimensional sub-manifold immersed in the entire state space of the
system. Instead of evolving in the whole state space the trajectories of the
system dynamics evolve on the lower-dimensional manifolds. As per the desired
behavior of the system or the set of constraints to be followed by the system,
the manifold can be selected as an equilibrium point or a particular region of
the state space of the system.

####  II-B1 Concept of Immersion

Since the focus is on the manifold stabilization problem and its application
in providing the controlled dynamical system perspective to the gradient
method via the P&I approach, the basic idea is to obtain an implicit manifold
by projecting the NLS system (system under consideration) onto a system with
pre-defined characteristics. For example, in the equilibrium stabilization
problem, the ultimate goal is to achieve the desired equilibrium in adaptive
control and stabilization or to obtain a zero equilibrium in the state
estimator design for the state estimation error. Since the P&I approach
utilizes the notion of immersion from Immersion and Invariance (I&I), the
concept of immersing the system dynamics into the defined lower-order target
dynamics is explained as follows.

Consider a NLS system

$\mathrm{\dot{x}=f(x)+g(x)u}$ (5)

where $\mathrm{x\in\mathbb{R}^{n}}$ and $\mathrm{u\in\mathbb{R}^{m}}$, the aim
is to find control law $\mathrm{u=\nu(x)}$ such that stabilizes the closed
loop system at an equilibrium point $\mathrm{x^{*}\in\mathbb{R}^{n}}$.

  * •

Target System: The above system can be projected onto a target dynamical
system

$\accentset{\mbox{\large.}}{\eta}=\beta(\eta)$ (6)

where $\eta\in\mathbb{R}^{\mathrm{p}}$ and $\mathrm{p<n}$, having a GAS
equilibrium at $\eta^{*}\in\mathbb{R}^{\mathrm{p}}$. Under a smooth mapping
$\pi:\mathbb{R}^{\mathrm{p}}\rightarrow\mathbb{R}^{\mathrm{n}}$,

$\mathrm{x^{*}}=\pi(\eta^{*})$ (7)

  * •

Manifold: The immersion of the system under consideration into lower-order
target dynamics (6) via mapping $\mathrm{x}=\pi(\eta)$ helps in the
construction of the implicit manifold $mathcal{M}$,

$\mathcal{M}=\left\\{\mathrm{x\in\mathbb{R}^{n}|x=\pi(\eta),\eta\in\mathbb{R}^{p}}\right\\}$
(8)

  * •

Invariance of the Manifold: The manifold
$\mathcal{M=\left\\{\mathrm{x\in\mathbb{R}^{n}|\pi(x)=0}\right\\}}$, with
$\pi(\mathrm{x})$ being smooth, is said to be invariant, if
$\mathrm{\pi(x(0))=0\Rightarrow\pi(x(t))=0,\;\;\forall t\geq 0}$. i.e. all
trajectories $\mathrm{x(t)}$ that start on the manifold remain on the manifold
and asymptotically converge to the point $\pi(0)$. (Refer to Definition A.10
of [20]) For proof of manifold invariance refer to Lemma 1.

  * •

Attractivity of the Manifold: Once the invariant manifold is obtained the
convergence of the off-the-manifold dynamics to the invariant manifold i.e.
attractivity of the manifold is ensured by linking the passivity theory with
immersion, which is accomplished by the P&I approach discussed in the
Appendix.

For further understanding of the immersion process and the choice of target
dynamics according to the desired objective, a geometric visualization is
depicted in Fig.1.

Figure 1: Schematic representation of the notion of immersion and its
projection onto the target dynamics

Fig. 1 illustrates two scenarios; equilibrium stabilization and orbital
stabilization (with the goal is to induce or generate a periodic orbit in
certain applications such as oscillation mechanisms in biology, electric
motors, power converters, and walking robots).

In Fig. 1(a), desired characteristics of the equilibrium stabilization for the
third-order system (system under consideration) are captured by second-order
target dynamics

$\accentset{\mbox{\large.}}{\eta}_{1}=-\beta_{1}\eta_{1}\hskip
28.45274pt\accentset{\mbox{\large.}}{\eta}_{2}=-\beta_{2}\eta_{2}.$ (9)

Similarly, Fig. 1(a) captures the desired characteristics of the periodic
orbit for the third-order system by second-order target dynamics.

$\accentset{\mbox{\large.}}{\eta}_{1}=\eta_{2}\hskip
28.45274pt\accentset{\mbox{\large.}}{\eta}_{2}=-\mathrm{sin}(\eta_{1})$ (10)

###### Remark 1

The objective of this article is to optimize the cost function. Thus, the
manifold stabilization problem is viewed as an equilibrium stabilization
problem, and focus is given to selecting the target dynamics that have a GES
equilibrium.

####  II-B2 Accelerated Gradient Method through P&I

Consider a second-order controlled dynamical system

$\accentset{\mbox{\large.}}{\mathrm{x}}_{1}=\mathrm{x_{2}}\hskip
28.45274pt\accentset{\mbox{\large.}}{\mathrm{x}}_{2}=\mathrm{u}$ (11)

The control $\mathrm{u}\in\mathbb{R}$ is designed such that the trajectories
of the above dynamical system will converge to the minima $\mathrm{x^{*}}$ of
the objective function $\mathrm{f(x)}$.

$(S_{1})$ Construction of the implicit manifold:

In the paper, the manifold is selected according to the desired objective,
i.e., to achieve the optimal solution $\mathrm{x^{*}}$ of (1) or the
equilibrium point of the dynamical system in (11). The optimization problem in
(1) is geometrically visualized as the manifold stabilization problem to
design an algorithm by adding the virtual control term $\mathrm{u}$. This
input $\mathrm{u}$ is evaluated through the procedure of passivity and
immersion (P&I) approach [21] (The step-wise procedure is overviewed in
Appendix) such that the trajectories converge exponentially to the equilibrium
with a suitable choice of an implicit manifold.

The target dynamics for finding the optimal solution $\mathrm{x^{*}}$ is
defined as $\mathrm{x_{2}=-\beta\nabla f(x_{1})}$. The target dynamics is
defined on the implicit manifold

$\psi(\mathrm{x_{1},x_{2}})=\mathrm{x_{2}+\beta\nabla f(x_{1})}=0$ (12)

###### Lemma 1

The implicit manifold $\psi$ in (12) is invariant.

Proof: According to the relationship $\mathrm{x_{2}=-\beta\nabla f(x_{1})}$,
the manifold is obtained as
$\psi(\mathrm{x_{1},x_{2}})=\mathrm{x_{2}+\beta\nabla f(x_{1})}=0$. The normal
to the velocity vector field $\begin{bmatrix}\dot{\mathrm{x}}_{1}\\\
\dot{\mathrm{x}}_{2}\end{bmatrix}$ is given as

$\nabla\psi(\mathrm{x_{1},x_{2}})=\begin{bmatrix}\beta\nabla^{2}\mathrm{f(x_{1})}&\mathrm{I}\end{bmatrix}$
(13)

The product of (13) and the velocity vector field yields
$\beta\nabla^{2}\mathrm{f(x_{1})}\dot{\mathrm{x}}_{1}+\dot{\mathrm{x}}_{2}\Rightarrow\beta\nabla^{2}\mathrm{f(x_{1})}\dot{\mathrm{x}}_{1}-\beta\nabla^{2}\mathrm{f(x_{1})}\dot{\mathrm{x}}_{1}=0$.
This implies that the velocity vector field is always tangent to the manifold
$\psi$. Hence, the implicit manifold $\psi$ is invariant.

$(S_{2})$ Tangent space structure for control systems:

The psuedo Riemann metric
$\mathrm{R(x)}=\mathrm{\nabla\psi(x_{1},x_{2})^{T}\nabla\psi(x_{1},x_{2})}$ is
defined.

$\begin{split}\mathrm{R(x)}=\begin{bmatrix}\mathrm{\beta\nabla^{2}f(x_{1})}\\\
\mathrm{I}\end{bmatrix}\begin{bmatrix}\mathrm{\beta\nabla^{2}f(x_{1})}&\mathrm{I}\end{bmatrix}=\\\
\begin{bmatrix}(\mathrm{\beta\nabla^{2}f(x_{1})})^{T}(\mathrm{\beta\nabla^{2}f(x_{1})})&\mathrm{\beta\nabla^{2}f(x_{1})}\\\
\mathrm{\beta\nabla^{2}f(x_{1})}&\mathrm{I}\end{bmatrix}\end{split}$ (14)

Comparing (14) with (44), $\mathrm{m_{21}=\beta\nabla^{2}f(x_{1})}$ and
$\mathrm{m_{22}=1}$. Hence this gives $\mathrm{\nabla
q(x)=\beta\nabla^{2}f(x_{1})}$ and $\mathrm{q(x)=\beta\nabla f(x_{1})}$. Since
$\accentset{\mbox{\large.}}{\mathrm{x}}_{2}$ is along the vertical direction,
the passive output is selected as
$\accentset{\mbox{\large.}}{\mathrm{x}}_{2}+\mathrm{m}_{22}^{-1}\mathrm{m}_{21}\accentset{\mbox{\large.}}{\mathrm{x}}_{1}=\accentset{\mbox{\large.}}{\mathrm{x}}_{2}+\mathrm{\beta\nabla^{2}f(x_{1})}\accentset{\mbox{\large.}}{\mathrm{x}}_{1}$
which is in the same direction or parallel to
$\accentset{\mbox{\large.}}{\mathrm{x}}_{2}$.

$(S_{3})$ Passive output: According to the preliminary definition of passive
output, $\mathrm{y=y_{1}+y_{2}}$ is opted such that,
$\mathrm{y_{1}}=\mathrm{\int_{0}^{t}(\nabla_{x_{1}}\varphi(x)\accentset{\mbox{\large.}}{\mathrm{x}}_{1})dt}=\mathrm{\beta\nabla
f(x_{1})}$, and
$\mathrm{y_{2}}=\mathrm{\int_{0}^{t}(\nabla_{x_{2}}\varphi(x)\accentset{\mbox{\large.}}{\mathrm{x}}_{2})dt}=\mathrm{x_{2}}$.

$(S_{4})$ Storage function: The storage function (or candidate Lyapunov
function) is defined with the passive output $\mathrm{y}$ as:

$\mathrm{S(x_{1},x_{2})=\frac{1}{2}(\mathrm{x_{2}+\beta\nabla f(x_{1})})^{2}}$
(15)

The convergence of the trajectories of the off-the-manifold dynamics to the
implicit manifold at an exponential rate $\alpha$ is ensured with the
condition :

$\mathrm{\accentset{\mbox{\large.}}{S}\leq-\alpha_{1}S}$ (16)

Utilizing the condition (16) along with the storage function (15) yeilds

$(\mathrm{x_{2}+\beta\nabla
f(x_{1})})(\accentset{\mbox{\large.}}{\mathrm{x}}_{2}+\mathrm{\beta\nabla^{2}f(x_{1})}\accentset{\mbox{\large.}}{\mathrm{x}}_{1})\leq-\frac{\alpha_{1}}{2}(\mathrm{x_{2}+\beta\nabla
f(x_{1})})^{2}$ (17)

$\Rightarrow\accentset{\mbox{\large.}}{\mathrm{x}}_{2}+\mathrm{\beta\nabla^{2}f(x_{1})}\accentset{\mbox{\large.}}{\mathrm{x}}_{1}+\frac{\alpha_{1}}{2}(\mathrm{x_{2}+\beta\nabla
f(x_{1})})=0$ (18)

From (11), the equation (18) is modified in terms of the final control law as:

$\mathrm{u}=-\mathrm{\beta\nabla^{2}f(x_{1})}\mathrm{x}_{2}-\frac{\alpha_{1}}{2}(\mathrm{x_{2}+\beta\nabla
f(x_{1})})$ (19)

###### Remark 2

For simplicity and understanding, the above formulation is illustrated for
$\mathrm{x}\in\mathbb{R}$. This approach can be extended to
$\mathrm{x}\in\mathbb{R}^{\mathrm{n}}$ with
$\mathrm{u}\in\mathbb{R}^{\mathrm{n}}$.

####  II-B3 ODE of Accelerated Gradient Method through P&I:

Hence with the final control law, the controlled dynamical system modifies to
the following form:

$\begin{split}\accentset{\mbox{\large.}}{\mathrm{x}}_{1}&=\mathrm{x_{2}}\\\
\accentset{\mbox{\large.}}{\mathrm{x}}_{2}&=-\mathrm{\beta\nabla^{2}f(x_{1})}\mathrm{x}_{2}-\frac{\alpha_{1}}{2}(\mathrm{x_{2}+\beta\nabla
f(x_{1})})\end{split}$ (20)

The above system can be written as a second-order ordinary differential
equation (ODE):

$\accentset{\mbox{\large..}}{\mathrm{x}}+\beta\nabla^{2}\mathrm{f(x)}\dot{\mathrm{x}}+\alpha\mathrm{\dot{x}}+\alpha\beta\nabla\mathrm{f(x)}=0$
(21)

where $\alpha=\frac{\alpha_{1}}{2}$.

###### Remark 3

The ODEs of accelerated gradient methods in literature reflect a geometric
damping term involving the Hessian of the objective function. In P&I the
Hessian term $(\beta\nabla^{2}\mathrm{f(x)}\dot{\mathrm{x}})$ appears from the
implicit manifold. The trajectories of the off-the-manifold dynamics take into
account the geometry of the implicit manifold and approach the invariant
manifold exponentially. The geometry of the manifold plays a crucial role in
the acceleration phenomenon. Hence it justifies the term “clever geometric
damping” defined by Alvarez et. al. [22]. If the manifold is changed, the
geometry will be altered, hence the Hessian damping term.

The above formulation is based on the invariant manifold theory which begins
with the assumption that the invariant manifold exists. A further theory is
developed from this point onwards by building upon this basic invariant
manifold. From a numerical analysis viewpoint, it is necessary to analyze
whether or not an invariant manifold persists [23] under perturbation. And if
so, if it maintains, loses, or gains differentiability is important.

For computational purposes, all the algorithms are implemented in the
discrete-time system. These calculations are often subjected to numerical
errors due to approximations. In order to enhance the accuracy in such
circumstances the invariant manifold should persist under the perturbations
(i.e., numerical errors).

###  II-C Persistence of Invariant Manifold

If the considered invariant manifold controlled with P&I possesses the
property of normal hyperbolicity, then the normally hyperbolic invariant
manifold (NHIM) persists under perturbations. The notion of normal
hyperbolicity and the conditions to analyze its persistence property is
introduced by N. Fenichel [24] and M. Hirsch, et.al. [25]. For a dynamical
system, an invariant manifold is considered as NHIM [26] if there is a
splitting of state space into three invariant sub-bundles: a bundle tangent to
the manifold, a stable bundle (dynamics on this is contracting), and an
unstable bundle (dynamics on this is expanding). If there exists a dominance
of dynamics on the stable and unstable bundle over the dynamics on the tangent
bundle, then the manifold is NHIM.

###### Remark 4

The normal hyperbolicity implies that the dynamics transverse to the invariant
manifold is dominant i.e., the growth rate of the vectors transverse to the
manifold dominates the growth rate of vectors tangent to the manifold. This
conceptual explanation is understood from the definition of normal
hyperbolicity and $r$-normally hyperbolicity (Refer [25] )

####  II-C1 Concept of Normally Controlled Hyperbolicity

Since the concerned invariant manifold is a controlled invariant manifold
hence the notion of $r$-normally controlled hyperbolicity introduced in [23]
should be considered.

Given a controlled invariant manifold $\mathbb{M}$, stabilised by the control
law $\mathrm{u}$ designed for the controlled dynamical system
$\mathrm{\accentset{\mbox{\large.}}{x}=X_{C}(x,u)}$, where
$\mathrm{x\in\Omega}$ and $\mathrm{u}\in\mathbb{U}$. The corresponding flow
for the controlled invariant manifold is
$f:\mathbb{M}\times\mathbb{U}\rightarrow\mathbb{M}$. Assume
$\sigma_{\mathrm{C}}$ chosen to associate to a state $\mathrm{x\in\Omega}$, a
stabilising control law $\mathrm{u}\in\mathbb{U}$.

###### Definition 3

(r-normally controlled-hyperbolicity)
In the above-mentioned controlled dynamical system, the manifold $\mathbb{M}$
is an r-normally controlled-hyperbolic with respect to the controlled flow

$\rho:=f\circ\sigma_{\mathrm{C}}$

iff the tangent bundle of $\Omega$ restricted to $\mathbb{M}$ splits into two
continuous sub-bundles:

$\mathbb{T}_{\mathbb{M}}\Omega=\mathbb{T}\mathbb{M}\oplus\mathbb{N}^{\mathrm{st}}$

invariant by the tangent of $\rho$ i.e. $\mathbb{T}{\rho}$, such that for
$p\in\mathbb{M}$, $0\leq\kappa\leq\mathrm{r}$:
$\mathbb{T}\rho$ contracts the stable normal bundle $\mathbb{N}^{\mathrm{st}}$
more sharply than $\mathbb{T}\mathbb{M}$, i.e.,

$\left\|(\mathbb{T}\rho\setminus\mathbb{N}_{p}^{\mathrm{un}}(\mathbb{M}))^{-1}\right\|\cdot\left\|\mathbb{T}\rho\setminus\mathbb{N}^{\mathrm{st}}(\mathbb{M})\right\|^{\kappa}\leq
c(1/\lambda)$

for some $1<\lambda<\infty,0<c<\infty$. Referring to the contraction and
expansion conditions, the concept of normal hyperbolicity can be informally
defined as the normal behavior of $\mathbb{T}\rho$ dominates the tangential
behavior of $\mathbb{T}\rho$. The dominance of the normal behavior over the
tangential behavior allows the controlled invariant manifold to persist
subjected to the perturbations.

####  II-C2  ODE of Accelerated Gradient Method through P&I and Additional
off-the-manifold Dynamics

To ensure the persistence of the invariant manifold under the perturbations, a
contracting off-the-manifold dynamics transverse to the invariant manifold
having the same equilibrium point as (21) is added. Hence the equation (21)
modifies to the following form:

$\ddot{\mathrm{x}}+\beta\nabla^{2}\mathrm{f(x)}\dot{\mathrm{x}}+\alpha\mathrm{\dot{x}}+\alpha\beta\nabla\mathrm{f(x)}+\left[\alpha\dot{\mathrm{x}}+\nabla\mathrm{f(x)}\right]=0$
(22)

The dynamics $\mathrm{\accentset{\mbox{\large.}}{x}=-\frac{1}{\alpha}\nabla
f(x)}$ transverse to the invariant manifold aids to the effect of contraction
of the stable normal bundle $\mathbb{N}^{\mathrm{st}}$. In effect, the
dynamical behavior normal to the invariant manifold is more dominant than the
form of the dynamics defined on the invariant manifold. The additional off-
the-manifold dynamics adhere to the notion of normal hyperbolicity of the
controlled invariant manifold. i.e. the controlled invariance properties are
preserved under perturbations.

###### Remark 5

In the considered problem, the optimal solution $\mathrm{x^{*}}$ is hyperbolic
(Refer [27]). Under normal hyperbolicity, even if the system perturbs, the new
point of a perturbed system preserves the properties of the interest.

Due to the inclusion of the dynamics transverse to the invariant manifold, the
space around the manifold becomes more contracting. The strong contraction
ensures the persistence of the invariant manifold and results in faster
convergence of the trajectories of the off-the-manifold dynamics towards the
invariant manifold. The persistence of the invariant manifold under the
perturbation ensures that the proposed algorithm is insensitive to numerical
errors.

###### Remark 6

The stability of the system remains unaffected by the addition of the dynamics
transverse to the manifold since the response of the additional dynamics is
exponential and its effect diminishes as the trajectories of the off-the-
manifold dynamics approach the equilibrium. Since the effect of dynamics
transverse to the manifold disappears at the equilibrium point, there is no
shift of the equilibrium point.

###  II-D Correlation with the Nesterov’s accelerated gradient method

The ODE derived with manifold stabilization procedure

$\accentset{\mbox{\large..}}{\mathrm{x}}+2\alpha\mathrm{\accentset{\mbox{\large.}}{x}}+\beta\nabla^{2}\mathrm{f(x)}\accentset{\mbox{\large.}}{\mathrm{x}}+(1+\alpha\beta)\nabla\mathrm{f(x)}=0$
(23)

is similar to Nesterov’s Accelerated Gradient for strongly convex functions
(NAG-SC) mentioned in the high-resolution differential equation derived by Bin
Shi et. al. [10]).

$\accentset{\mbox{\large..}}{\mathrm{x}}+2\sqrt{\mu}\mathrm{\accentset{\mbox{\large.}}{x}}+\sqrt{s}\nabla^{2}\mathrm{f(x)}\accentset{\mbox{\large.}}{\mathrm{x}}+(1+\sqrt{\mu
s})\nabla\mathrm{f(x)}=0$ (24)

The Hessian of the function in the above equation works as the geometric
damping and aids in the acceleration of the GD algorithm. The algorithms
designed in the continuous time system are implemented with the help of
numerical integration in the discrete-time system for computation purposes. To
reduce the computational burden on the processor, the computation complexity
should be reduced. To achieve reduced complexity, the selection of the step
size is a matter of concern. The selection of the larger step size allows the
reduction in the number of iterations. But the larger step size poses a
serious problem to the stability of the algorithm. The Hessian of the function
gives information about the curvature of the descent resulting in faster
convergence by appropriate scaling of the step size as per the level of the
curvature.

The Hessian term reflects in the P&I approach by following a systematic
methodology without any approximations. The Hessian damping term in the ODE
derived through the P&I is governed by the geometry of the manifold and allows
the convergence to the solution with a larger step size without disturbing the
stability of the algorithm. Larger step sizes ensure the convergence to the
optimal solution with a lesser number of iterations leading to accelerated
algorithms.

###  II-E Significance of $\alpha$ and $\beta$

The parameters $\alpha$ and $\beta$ in the proposed approach can be selected
by comparing the ODE (23) to the parameters $\mu$ and $s$ in the NAG-sc. The
parameters $\mu$ (the modulus of convexity of the strongly convex function)
and $s$ (step size) appear from the discretization. These parameters play an
important role in the stability of the algorithm and are selected by following
several guidelines.

The parameter $\alpha$ present in the off-the-manifold dynamics is bounded by
the modulus of convexity of the function. $\alpha$ is related to the manifold
attractivity (the trajectories of the off-the-manifold dynamics approach the
implicit manifold exponentially with the rate $\alpha$) and $\mu$ represents
the lower bound to the growth of the function. The trajectories of the off-
the-manifold dynamics should approach the invariant manifold faster than this
bound. Hence $\alpha$ can be related to $\mu$. The other parameter $\beta$
present in the dynamics defined on the invariant manifold is bounded by the
Lipschitz constant of the function. $\beta$ in the definition of the implicit
manifold is similar to the step size $\mathrm{s}$ in the GD iterative
algorithm. The optimal step size is bounded by the Lipschitz constant of the
function (Refer to equation (1.2.12) in [28]). Referring to the NAG-SC in
equation (24), the optimal values of the parameters are selected as
$\alpha=\sqrt{\mu}$ and $\beta=\sqrt{s}$.

###  II-F Stability Analysis

The ODE (22) for the accelerated algorithm can be looked upon as a second-
order linear dynamical system as follows:

$\begin{bmatrix}\accentset{\mbox{\large.}}{\mathrm{x}}_{1}\\\
\accentset{\mbox{\large.}}{\mathrm{x}}_{2}\end{bmatrix}=\begin{bmatrix}0&1\\\
0&0\end{bmatrix}\begin{bmatrix}\mathrm{x_{1}}\\\
\mathrm{x_{2}}\end{bmatrix}+\begin{bmatrix}0&0\\\
1&1\end{bmatrix}\begin{bmatrix}\mathrm{u_{1}}\\\ \mathrm{u_{2}}\end{bmatrix}$
(25)

Since the dynamical system is linear, by superposition property, both the
control inputs in the above system $\mathrm{u_{1}}$ and $\mathrm{u_{2}}$ can
be evaluated independently. Control input $\mathrm{u_{1}}$ is evaluated
through the P&I procedure, whereas the other control input $\mathrm{u_{2}}$ is
added for the persistence of the invariant manifold under perturbations. Hence
the additional off-the-manifold dynamics transversal to the invariant manifold
will formulate the control input $\mathrm{u_{2}}=\mathrm{-\alpha x_{2}-\nabla
f(x_{1})}$.

The Lyapnouv function for the system dynamics with $\mathrm{u_{1}}$ evaluated
through the P&I strategy is mentioned earlier in equation (15). The
convergence of the trajectories of the off-the-manifold dynamics to the
invariant manifold with the control input $\mathrm{u_{1}}$ is ensured with the
condition (16).

The stability of the system with control input $\mathrm{u_{2}}$ can be
evaluated with $\mathrm{u_{1}}=0$. The second order dynamical system with
control input $\mathrm{u_{2}}$ will be of a form:

$\begin{split}\accentset{\mbox{\large.}}{\mathrm{x}}_{1}&=\mathrm{x_{2}}\\\
\accentset{\mbox{\large.}}{\mathrm{x}}_{2}&=-\mathrm{\alpha x_{2}-\nabla
f(x_{1})}\end{split}$ (26)

For the sake of simplicity and understanding, the stability analysis is
illustrated for $\mathrm{x_{1},x_{2}\in\mathbb{R}}$. A similar analysis can be
suitably extended to $\mathrm{x_{1},x_{2}\in\mathbb{R}^{n}}$.

###### Lemma 2

The second-order dynamics in (26) is stable.

Proof: The Lyapunov function for the second order system mentioned in (26)
with the additional off-the-manifold dynamics transversal to the invariant
manifold is given as:

$\mathrm{V(x_{1},x_{2})}=\mathrm{f(x_{1})-f(x_{1}^{*})+\frac{1}{2}{x}_{2}^{2}}$
(27)

For proving stability, a derivative of the above Lyapunov function is
evaluated as:

$\mathrm{\accentset{\mbox{\large.}}{V}(x_{1},x_{2})}=\mathrm{\nabla
f(x_{1})\accentset{\mbox{\large.}}{x}_{1}+x_{2}\accentset{\mbox{\large.}}{x}_{2}}$
(28)
$\Rightarrow\mathrm{\accentset{\mbox{\large.}}{V}(x_{1},x_{2})}=\mathrm{\nabla
f(x_{1})x_{2}-\alpha x_{2}^{2}-\nabla f(x_{1})x_{2}}$ (29)
$\Rightarrow\mathrm{\accentset{\mbox{\large.}}{V}(x_{1},x_{2})}=-\alpha\mathrm{x_{2}^{2}}$
(30)

Since $\mathrm{\dot{V}(x_{1},x_{2})}\leq 0$ the second-order system dynamics
(26) is stable.

###### Remark 7

The asymptotic stability of the above system can be proved using LaSalle’s
theorem.

###### Remark 8

For the exponential stability of the system in (26), the Lyapunov function
from [10] can be referred. From the Theorem 2 of [10] the Lyapunov function
can be selected as

$\mathrm{V(x_{1},x_{2})}=\mathrm{f(x_{1})-f(x_{1}^{*})+\frac{1}{2}{x}_{2}^{2}+\frac{1}{2}(x_{2}+\alpha(x_{1}-x_{1}^{*}))^{2}}$

where $\alpha=\alpha_{1}/2$.

By following the same steps as mentioned in [10], the exponential stability
can be proved with the help of the above Lyapunov function. The exponential
convergence rates for both the dynamics with $\mathrm{u_{1}}$ and
$\mathrm{u_{2}}$ are the same i.e. $\alpha_{1}$.

##  III Numerical Example and Results

To gain deeper insights into the proposed approach consider an optimization
problem with a quadratic cost function.

$\begin{split}&\mathrm{min}\;\mathrm{f}(\mathrm{x})\\\
&\mathrm{s.t.}\;\mathrm{x}\in\mathrm{X}\end{split}$ (31)

where $\mathrm{X}=\left\\{\mathrm{x}\in\mathbb{R}\right\\}$ and
$\mathrm{f(x)=\frac{1}{2}x^{T}x}$.

The gradient descent for the above quadratic objective function will be:

$\accentset{\mbox{\large.}}{\mathrm{x}}=-\mathrm{\nabla f}=-\mathrm{x}$ (32)

For solving the given problem with the control-theoretic approach consider a
second-order system

$\displaystyle\begin{split}\mathrm{\accentset{\mbox{\large.}}{x}}_{1}=\mathrm{x_{2}}\hskip
28.45274pt\mathrm{\accentset{\mbox{\large.}}{x}}_{2}=\mathrm{u}\end{split}$
(33)

As per the equation (12), the implicit manifold is defined as

$\psi(\mathrm{x_{1},x_{2}})=\mathrm{x_{2}+\beta x_{1}}=0$ (34)

Following the remaining steps of P&I the control law is derived.

$\mathrm{u}=\mathrm{-\beta x_{2}-\alpha x_{2}-\alpha\beta x_{1}}$ (35)

The corresponding second-order ODE is

$\mathrm{\accentset{\mbox{\large..}}{x}}+\beta\mathrm{\accentset{\mbox{\large.}}{x}}+\alpha\mathrm{\accentset{\mbox{\large.}}{x}}+\alpha\beta\mathrm{x}=0$
(36)

###  III-A Response before the addition of dynamics transversal to the
invariant manifold

Figure 2: Response of Accelerated Gradient Descent derived through P&I with
the y-axis representing the solution of the quadratic cost function
$\mathrm{f(x)}$ and x-axis representing the time for $\alpha=10$ with various
values of $\beta$ ** Figure 3: Response of Accelerated Gradient Descent
derived through P&I with the y-axis representing the solution of the quadratic
cost function $\mathrm{f(x)}$ and x-axis representing the time for $\alpha=1$
with various values of $\beta$ **

The ODE (36) is simulated with the various values of the parameters $\alpha$
and $\beta$ as demonstrated in Fig. 2 and Fig. 3. It is observed that the
convergence to the optimal solution $\mathrm{x^{*}}$ becomes faster with the
increasing value of the parameter $\beta$. From Fig. 2, it is visible that the
peak overshoots are decreased as the value of $\alpha$ is increased as
compared to Fig. 3.

###  III-B Response after the addition of dynamics transversal to the
invariant manifold

For ensuring the persistence of the invariant manifold under the perturbations
such as numerical errors, another dynamics transverse to the invariant
manifold is added. Hence the ODE (36) modifies as:

$\mathrm{\accentset{\mbox{\large..}}{x}}+\beta\mathrm{\accentset{\mbox{\large.}}{x}}+2\alpha\mathrm{\accentset{\mbox{\large.}}{x}}+(1+\alpha\beta)\mathrm{x}=0$
(37)

Figure 4: Effect of the addition of another dynamics transverse to the
invariant manifold to the Accelerated Gradient Descent derived through P&I
with the y-axis representing the solution of the quadratic cost function
$\mathrm{f(x)}$ and x-axis representing the time for $\alpha=10$ with various
values of $\beta$ ** Figure 5: Effect of the addition of another dynamics
transverse to the invariant manifold to the Accelerated Gradient Descent
derived through P&I with the y-axis representing the solution of the quadratic
cost function $\mathrm{f(x)}$ and x-axis representing the time for $\alpha=1$
with various values of $\beta$ **

The above ODE is simulated for the different values of $\alpha$ and $\beta$
similar to the earlier case and the results are illustrated in Fig. 4 and Fig.
5. For the given quadratic cost function the optimal values of the parameters
are selected as $\alpha=1$ and $\beta=0.9$.

The effect of the additional off-the-manifold dynamics transverse to the
invariant manifold with the optimal values of parameters $\alpha$ and $\beta$
is demonstrated in Fig. 5. It is observed that due to the dominating effect of
the dynamics transverse to the invariant manifold over the dynamics defined on
the invariant manifold, the convergence to the optimal solution
$\mathrm{x^{*}}$ of the quadratic cost function is achieved with the minimum
number of iterations.

###### Remark 9

** Fig. 2 to Fig. 5 are revised in the latest version of the manuscript. The
ODE response is illustrated in continuous time domain and the axes labels are
revised accordingly.

##  IV Discussions

In this paper, a controlled dynamical system perspective for deriving the
accelerated gradient method is proposed. The perspective is motivated by an
understanding that the algorithms are numerical integration of a dynamical
system. Hence for achieving a particular objective, any problem could be
looked upon as a dynamical system. And a control law is designed for the
dynamical system to achieve the desired objective. The controlled dynamical
system perspective is implemented with the incorporation of the notion of
manifold stabilization. A control scheme for the dynamical system is designed
to ensure the convergence (attractivity) of the system dynamics to a
particular invariant manifold. Hence the considered problem is also referred
to as the controlled invariant manifold approach. Since the objective of this
paper is to optimally reach the optimal solution $\mathrm{x^{*}}$ where
$\mathrm{\nabla f(x)=0}$, the manifold is considered as an equilibrium point.
Hence the considered problem is reduced to an equilibrium stabilization
problem.

The equilibrium stabilization is achieved by blending the concept of immersion
and passivity through the P&I approach. The invariant manifold is constructed
through the process of immersion by immersing the system dynamics into the
specified lower-order target dynamics. The convergence of the trajectories of
the system dynamics to the invariant manifold (attractivity of the invariant
manifold) is derived through the passivity theory. The ODE derived through the
systematic procedure of P&I contains a term with Hessian of the function which
provides the knowledge of the local geometry of the curvature and enables the
convergence to the optimal solution in fewer iterations with the selection of
larger step sizes.

Since the algorithms are implemented through the various techniques of
numerical integration in the discrete-time system, some numerical errors occur
due to the approximations involved in the process. The invariant manifold is
often susceptible to these numerical errors (perturbations) and the invariant
manifold might not persist under such perturbations. The proposed theory is
developed on the basic invariant manifold and this invariant manifold should
persist when subjected to numerical errors. Hence, to ensure the persistence
of the controlled invariant manifold, an additional contracting dynamics
transversal to the manifold is added which constitutes another major
contribution of this paper. The additional off-the-manifold dynamics adhere to
the notion of normal hyperbolicity of the controlled invariant manifold. Under
normal hyperbolicity, even if the system perturbs, the new point of the
perturbed system preserves the properties of the interest. The additional
dynamics ensure the stronger contraction of the ambient space and hence the
persistence of the invariant manifold under numerical perturbations. As a
result, a numerically insensitive algorithm for the accelerated gradient
method is developed.

It is observed that the ODE obtained through the proposed approach with the
additional off-the-manifold dynamics resembles the high-resolution ODE for the
NAG-SC method. The basic formulation derived for the NAG-SC could be modified
to other settings for accelerated gradient methods. For various objective
functions, the Hessian of the function term is negligible and could be
omitted. This leads to the equation which is similar to the high-resolution
ODE for the HB algorithm. Since the proposed formulation is developed on the
persistent invariant manifold, the algorithm remains insensitive to the
additional terms to the dynamics. Hence with the addition of another damping
term, the proposed formulation could be extended to the accelerated Triple
Momentum (TM) method [29]. These extensions are summarised with the
corresponding ODEs in the following figure.

Figure 6: Extensions of ODE obtained through proposed approach

The term $[\alpha\accentset{\mbox{\large.}}{\mathrm{x}}+\nabla\mathrm{f(x)}]$
appears in NAG-SC, HB and TM methods. This certifies that the additional
dynamics is necessary for the accelerated gradient methods. By combining the
additional dynamics with the remaining ODEs the above equations are rewritten
as:

  * •

Proposed ODE:

$\accentset{\mbox{\large..}}{\mathrm{x}}+2\alpha\mathrm{\accentset{\mbox{\large.}}{x}}+\beta\nabla^{2}\mathrm{f(x)}\accentset{\mbox{\large.}}{\mathrm{x}}+(1+\alpha\beta)\nabla\mathrm{f(x)}=0$

  * •

Nesterov’s Accelerated Gradient for Strongly Convex function (NAG-SC):

$\accentset{\mbox{\large..}}{\mathrm{x}}+2\sqrt{\mu}\mathrm{\accentset{\mbox{\large.}}{x}}+\sqrt{s}\nabla^{2}\mathrm{f(x)}\accentset{\mbox{\large.}}{\mathrm{x}}+(1+\sqrt{\mu
s})\nabla\mathrm{f(x)}=0,$

where the parameters of the proposed ODE are selected as $\alpha=\sqrt{\mu}$
and $\beta=\sqrt{s}$.

  * •

Heavy ball method (HB):

$\accentset{\mbox{\large..}}{\mathrm{x}}+2\sqrt{\mu}\mathrm{\accentset{\mbox{\large.}}{x}}+(1+\sqrt{\mu
s})\nabla\mathrm{f(x)}=0$

where the parameters of the proposed ODE are selected as $\alpha=\sqrt{\mu}$
and $\beta=\sqrt{s}$ and $\nabla^{2}\mathrm{f(x)}=0$.

  * •

Accelerated Tripple Momentum method (TM):

$\accentset{\mbox{\large..}}{\mathrm{x}}+2\sqrt{\mu}\mathrm{\accentset{\mbox{\large.}}{x}}+\gamma(1+\sqrt{\mu
s})\sqrt{s}\nabla^{2}\mathrm{f(x)}\accentset{\mbox{\large.}}{\mathrm{x}}+(1+\sqrt{\mu
s})\nabla\mathrm{f(x)}=0$

where the parameters of the proposed ODE are selected as $\alpha=\sqrt{\mu}$
and $\beta=\sqrt{s}$ and another damping coefficient $\gamma$ is added.

###### Remark 10

The stability certificate of these methods can be referred from the proofs
given in [10] [29].

##  V Conclusion and Future Work

The controlled invariant manifold approach aims to promote a novel philosophy
to provide a better understanding of the working of the algorithm. The
proposition forwarded by this work is problems in optimization could be solved
by visualizing them as a manifold stabilization problem. And the algorithm
could be generated systematically through the P&I approach to obtain a
controlled dynamical system that achieves a desired objective. The authors are
intending to extend the proposed approach to the geodesic convexity settings.
In future work, the above setting proposed for the unconstrained optimization
case will be extended to the primal-dual optimization problems.

## Acknowledgment

The authors of the paper would like to acknowledge the support provided by
Control and Decision Research Center (CDRC), VJTI, Mumbai, India, and (Project
20627-X2001, FY22-Stankovic Fusion Methods) SLAC National Accelerator
Laboratory, Menlo Park, CA, USA.

## Appendix

###  V-A Existing Approaches for Accelerated Gradient Method

Consider an unconstrained optimization problem of the form

$\mathrm{min\;f(x)}\;\;\;\;\;\;\mathrm{x\in\mathbb{R}^{n}}$

The objective function $\mathrm{f(x)}$ can be a convex function or a strongly
convex function. The optimal value of this function is denoted as
$\mathrm{f^{*}=\underset{x}{min}\;f(x)}$ and the optimal solution of the
above-mentioned problem is represented as
$\mathrm{x^{*}=arg\>\underset{x}{min}\;f(x)}$. The straightforward approach
for finding $\mathrm{x^{*}}$ is to solve the system of equation
$\mathrm{\nabla f(x^{*})=0}$ which is usually not easier to solve from an
analytical viewpoint. Hence, a preferable strategy is an iterative scheme that
involves computing a minimizing sequence of points
$\mathrm{x^{(0)},x^{(1)},\cdots}$ such that $\mathrm{f(x^{(k)})\rightarrow
f(x^{*})}$ as $\mathrm{k\rightarrow\infty}$.

Gradient descent (GD) or steepest descent is a ubiquitous iterative algorithm
for solving optimization problems. It is represented by the first-order
dynamics in the time system as follows:

$\begin{split}\mathrm{\accentset{\mbox{\large.}}{x}(t)+\nabla f(x(t))=0}\\\
\mathrm{x(0)=x_{0}}\end{split}$ (38)

Since GD is a first-order system, its trajectory is solely governed by the
initial point $\mathrm{x_{0}}$. The GD system mimics the movement of a drop of
water as it slides on the profile defined by function $\mathrm{f(x)}$. The GD
algorithm is slower and it often gets stuck in local minima due to
inappropriate selection of step size.

The limitations of GD are overcome by the Heavy Ball (HB) method [7] proposed
by B. T. Polyak which is the accelerated version of the GD represented as:

$\mathrm{\accentset{\mbox{\large..}}{x}(t)+\nabla f(x(t))=0}$ (39)

HB method provides an accelerated response with oscillations. To suppress
these oscillations H. Attouch et. al. [9] added a damping term representing
friction and the modified method is referred to as Heavy Ball with Friction
(HBF) system. It is a non-linear oscillator with damping, given as:

$\mathrm{\accentset{\mbox{\large..}}{x}(t)+\lambda\accentset{\mbox{\large.}}{x}(t)+\nabla
f(x(t))=0}$ (40)

###  V-B The P&I Approach

The immersion of the system under consideration into lower-order target
dynamics helps in the construction of the implicit manifold. The convergence
and manifold attractivity are ensured by linking the passivity theory with
immersion. The P&I approach utilizes the idea of an invariant target manifold
giving rise to a non-degenerate two form, through which the definition of
certain passive outputs and storage functions leads to a generation of control
law for stabilizing the system. Moreover, the involvement of a particular
structure of tangent bundle associated with a suitable pseudo-Riemannian
structure imposed on the control system makes the P&I more methodical and
systematic.

The notion of the convergence of the off-the-manifold dynamics to the implicit
manifold and manifold attractivity can be visualized in Fig. 7.

Figure 7: Geometric representation of the aspects of invariant manifold theory

The P&I approach’s concept and procedures applicable to this research are
presented below. Given a single input system

$\displaystyle\begin{split}\accentset{\mbox{\large.}}{\mathrm{x}}=\mathrm{\mathbb{F}(\mathrm{x},\lambda)}\hskip
42.67912pt\accentset{\mbox{\large.}}{\lambda}=\mathrm{u}\end{split}$ (41)

with $(\mathrm{x},\lambda)\in(\mathbb{R}^{n-1},\mathbb{R})$ and without any
particular structure. The proposed P&I approach ensures NLS stabilization by
executing the four steps $(S_{1}-S_{4})$ outlined below.
$(S_{1})$ Construction of the implicit manifold:

The target dynamics $\accentset{\mbox{\large.}}{\eta}=\beta(\eta)$ with
$\mathrm{x}=\pi(\eta)=\eta x$ is defined such that the subsystem
$\accentset{\mbox{\large.}}{\mathrm{x}}=\mathrm{\mathbb{F}(\mathrm{x},\varphi(\mathrm{\mathrm{x}}))}$
for $\mathrm{C}^{\infty}$ mapping
$\varphi(\mathrm{x}):\mathbb{R}^{\mathrm{n}}\rightarrow\mathbb{R}$ has a
GES/GAS equilibrium at the origin by considering the relationship
$\mathrm{\lambda}=\varphi(\mathrm{\mathrm{x}})$. This defines the implicit
manifold
$\mathbb{M}=\left\\{(\mathrm{x},\lambda)\in\mathbb{R}^{\mathrm{n-1}}\times\mathbb{R}|\Psi(\mathrm{x},\lambda):=\lambda-\varphi(\mathrm{x})=0\right\\}$,
and $\pi(\eta)=\mathrm{col}(\eta,\varphi(\mathrm{\eta}))$.

$(S_{2})$ Tangent space structure for control systems: Consider an
n-dimensional manifold $\mathbb{M}$ with tangent bundle
$\mathbb{T}_{\mathbb{M}}$, such that all $\mathrm{p}\in{\mathbb{M}}$,
$\mathbb{T}_{\mathrm{p}}{\mathbb{M}}$ has the following structure

$\mathbb{T}_{\mathrm{p}}{\mathbb{M}}=\mathbb{H}_{\mathrm{p}}\oplus\mathbb{V}_{\mathrm{p}}:\hskip
8.5359pt\mathbb{H}_{\mathrm{p}}\cap\mathbb{V}_{\mathrm{p}}=0$ (42)

where $\mathbb{H}_{\mathrm{p}}$ is the horizontal space and
$\mathbb{V}_{\mathrm{p}}$ is the vertical space.

$\mathrm{Then}\hskip
25.6073pt\mathbb{T}_{\mathrm{p}}{\mathbb{M}}=\mathbb{H}_{\mathrm{p}}\oplus\mathbb{V}_{\mathrm{p}}=(\accentset{\mbox{\large.}}{\mathrm{x}},0)\oplus(0,\accentset{\mbox{\large.}}{\lambda})=(\accentset{\mbox{\large.}}{\mathrm{x}},\accentset{\mbox{\large.}}{\lambda})$
(43)

is written for given system (41) at any point $\mathrm{p}\in\mathbb{M}$. With
the implicit manifold $\Psi(\mathrm{x},\lambda)$ obtained in $S_{1}$, the
normal vector direction is given by $\triangledown\Psi(\mathrm{x},\lambda)$.
Thus, a PR metric $\mathrm{R}$ on space $\mathbb{T}_{\mathrm{p}}{\mathbb{M}}$
can be defined as

$\displaystyle\begin{split}\mathrm{R}&=\triangledown\Psi(\mathrm{x},\lambda)^{\mathrm{T}}\triangledown\Psi(\mathrm{x},\lambda)\\\
&=\begin{bmatrix}\left(\frac{\partial\varphi}{\partial\mathrm{x}}\right)^{\mathrm{T}}\left(\frac{\partial\varphi}{\partial\mathrm{x}}\right)&\left(-\frac{\partial\varphi}{\partial\mathrm{x}}\right)^{\mathrm{T}}\\\
-\left(\frac{\partial\varphi}{\partial\mathrm{x}}\right)&\mathrm{I}\end{bmatrix}=\begin{bmatrix}\mathrm{\mathrm{m}}_{11}&\mathrm{\mathrm{m}}_{12}\\\
\mathrm{\mathrm{m}}_{21}&\mathrm{\mathrm{m}}_{22}\end{bmatrix}.\end{split}$
(44)

which is intuitively a natural choice.

###### Remark 11

To obtain the passive output, the metric $\chi$ is replaced with semi-
Riemannian metric $\mathrm{R}$ as a natural choice.

For $(\mathbb{M},\mathrm{R})$, the splitting is visualized as follows:

$(\accentset{\mbox{\large.}}{\mathrm{x}},\accentset{\mbox{\large.}}{\lambda})=\left(\accentset{\mbox{\large.}}{\mathrm{x}},-\mathrm{m}_{22}^{-1}\mathrm{m}_{21}\accentset{\mbox{\large.}}{\mathrm{x}}\right)\oplus\left(0,\accentset{\mbox{\large.}}{\lambda}+\mathrm{m}_{22}^{-1}\mathrm{m}_{21}\accentset{\mbox{\large.}}{\mathrm{x}}\right)=\widetilde{\mathbb{H}}_{\mathrm{p}}\oplus\widetilde{\mathbb{V}}_{\mathrm{p}}$
(45)

As $\accentset{\mbox{\large.}}{\lambda}$ is along the vertical direction, the
passive output is chosen as a component
$\accentset{\mbox{\large.}}{\lambda}+\mathrm{m}_{22}^{-1}\mathrm{m}_{21}\accentset{\mbox{\large.}}{\mathrm{x}}$
which is in the same direction or parallel to
$\accentset{\mbox{\large.}}{\lambda}$ [21]. Roughly speaking, the idea is to
bring the component
$\accentset{\mbox{\large.}}{\lambda}+\mathrm{m}_{22}^{-1}\mathrm{m}_{21}$ of
$\widetilde{\mathbb{H}}_{\mathrm{p}}$ to the component
$\mathrm{m}_{22}^{-1}\mathrm{m}_{21}\accentset{\mbox{\large.}}{\mathrm{x}}$ of
$\widetilde{\mathbb{V}}_{\mathrm{p}}$.

Figure 8: Geometrical interpretation: vertical vector
$\mathbb{V}_{\mathrm{p}}$ is along the fiber direction and
$\mathbb{H}_{\mathrm{p}}\oplus\mathbb{V}_{\mathrm{p}}=\mathbb{T}_{\mathrm{p}}{\mathbb{M}}$.

The geometrical interpretation for the splitting tangent vector is shown in
Fig. 8.

$(S_{3})$ Passive output: The component of $\mathrm{u}$ tangent vector along
$\accentset{\mbox{\large.}}{\lambda}$ is used to define the passive output

$\mathrm{y}=\mathrm{y}_{1}+\mathrm{y}_{2}=\int_{0}^{t}\accentset{\mbox{\large.}}{\lambda}\hskip
2.27626pt\mathrm{dt}+{\int_{0}^{t}(\mathrm{m}_{22}^{-1}\mathrm{m}_{21}\accentset{\mbox{\large.}}{\mathrm{x}})\hskip
2.27626pt\mathrm{dt}}$ (46)

with the help of the passivity theory. Here,
$\mathrm{y}_{1}={\int_{0}^{t}\accentset{\mbox{\large.}}{\lambda}\hskip
2.27626pt\mathrm{dt}}$ and
$\mathrm{y}_{2}={\int_{0}^{t}(\mathrm{m}_{22}^{-1}\mathrm{m}_{21}\accentset{\mbox{\large.}}{\mathrm{x}})\hskip
2.27626pt\mathrm{dt}}$ are defined. If $\mathrm{m}_{22}^{-1}\mathrm{m}_{21}$
is a constant then
$\mathrm{y}=({\lambda}+\mathrm{m}_{22}^{-1}\mathrm{m}_{21}{\mathrm{x}})$ is
defined. If $\mathrm{m}_{22}^{-1}\mathrm{m}_{21}$ is a function of
$\mathrm{x}$ then, it can be written as the gradient of any function
$\mathrm{q(x)}$ i.e.,
$\mathrm{m}_{22}^{-1}\mathrm{m}_{21}\mathrm{(x)}=\triangledown\mathrm{q(x)}$.
Then

$\mathrm{y}=\int_{0}^{t}(\accentset{\mbox{\large.}}{\lambda}+\triangledown\mathrm{q(x)}\accentset{\mbox{\large.}}{\mathrm{x}})\mathrm{dt}=(\lambda+\mathrm{q(x)})$
(47)

$(S_{4})$ Storage function: With $\mathrm{y}$, the candidate Lyapunov function
$\mathbb{S}(\mathrm{x},\lambda)$ (i.e., storage function) is defined as

$\mathbb{S}(\mathrm{x},\lambda)=\frac{1}{2}{\mathrm{y}}^{2}=\frac{1}{2}(\lambda+\mathrm{q(x)})^{2}.$
(48)

One can use the condition
$\accentset{\mbox{\large.}}{\mathbb{S}}\leq-\alpha\mathbb{S}$ along with the
storage function (48) and passive output (47) to get

$(\accentset{\mbox{\large.}}{\lambda}+\frac{\partial\mathrm{q(x)}}{\partial\mathrm{x}}\accentset{\mbox{\large.}}{\mathrm{x}})+\frac{\alpha}{2}(\lambda+\mathrm{q(x)})\leq
0$ (49)

The equation (49) is modified by substituting (41) in terms of the final
control law as

$\mathrm{u}=-\frac{\alpha}{2}\lambda-\frac{\alpha}{2}\mathrm{q(x)}-\frac{\partial\mathrm{q(x)}}{\partial\mathrm{x}}\mathrm{f(\mathrm{x},\lambda)}.$
(50)

The above-defined control law ensures the GAS equilibrium point of the system
to zero/origin.

###### Remark 12

The mathematical details & related theorems-proofs of the P&I approach can be
found in [21].

## References

  * [1] C. Hu, W. Pan, and J. Kwok, “Accelerated gradient methods for stochastic optimization and online learning,” _Advances in Neural Information Processing Systems_ , vol. 22, 2009.
  * [2] G. Lan, “An optimal method for stochastic composite optimization,” _Mathematical Programming_ , vol. 133, no. 1-2, pp. 365–397, 2012.
  * [3] S. Ghadimi and G. Lan, “Accelerated gradient methods for nonconvex nonlinear and stochastic programming,” _Mathematical Programming_ , vol. 156, no. 1-2, pp. 59–99, 2016.
  * [4] H. Li and Z. Lin, “Accelerated proximal gradient methods for nonconvex programming,” _Advances in neural information processing systems_ , vol. 28, 2015.
  * [5] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” _SIAM journal on imaging sciences_ , vol. 2, no. 1, pp. 183–202, 2009.
  * [6] G. Lan, Z. Lu, and R. D. Monteiro, “Primal-dual first-order methods with iteration-complexity for cone programming,” _Mathematical Programming_ , vol. 126, no. 1, pp. 1–29, 2011.
  * [7] B. T. Polyak, “Some methods of speeding up the convergence of iteration methods,” _Ussr computational mathematics and mathematical physics_ , vol. 4, no. 5, pp. 1–17, 1964.
  * [8] Y. E. Nesterov, “A method of solving a convex programming problem with convergence rate $o\left(\frac{1}{k^{2}}\right)$,” in _Doklady Akademii Nauk_ , vol. 269, no. 3. Russian Academy of Sciences, 1983, pp. 543–547.
  * [9] H. Attouch, X. Goudou, and P. Redont, “The heavy ball with friction method, i. the continuous dynamical system: global exploration of the local minima of a real-valued function by asymptotic analysis of a dissipative dynamical system,” _Communications in Contemporary Mathematics_ , vol. 2, no. 01, pp. 1–34, 2000.
  * [10] B. Shi, S. S. Du, M. I. Jordan, and W. J. Su, “Understanding the acceleration phenomenon via high-resolution differential equations,” _Mathematical Programming_ , pp. 1–70, 2021.
  * [11] U. Helmke and J. B. Moore, _Optimization and dynamical systems_. Springer Science & Business Media, 2012.
  * [12] R. W. Brockett, “Dynamical systems that sort lists, diagonalize matrices, and solve linear programming problems,” _Linear Algebra and its applications_ , vol. 146, pp. 79–91, 1991.
  * [13] J. Soto and J.-J. E. Slotine, “Nonlinear contraction tools for constrained optimization,” _arXiv preprint arXiv:1206.1838_ , 2012.
  * [14] L. Lessard, B. Recht, and A. Packard, “Analysis and design of optimization algorithms via integral quadratic constraints,” _SIAM Journal on Optimization_ , vol. 26, no. 1, pp. 57–95, 2016.
  * [15] A. Wibisono, A. C. Wilson, and M. I. Jordan, “A variational perspective on accelerated methods in optimization,” _proceedings of the National Academy of Sciences_ , vol. 113, no. 47, pp. E7351–E7358, 2016.
  * [16] S. Shadab Nayyer, S. R. Wagh, and N. M. Singh, “Passivity and immersion (p&i) approach for constructive stabilization and control of nonlinear systems,” _IEEE Control Systems Letters_ , vol. 7, pp. 817–822, 2023.
  * [17] S. S. Nayyer, J. Hozefa, R. Gunjal, S. Bhil, S. Wagh, and N. Singh, “Passivity and immersion (p&i) approach with gaussian process for stabilization and control of nonlinear systems,” _IEEE Access_ , vol. 10, pp. 132 621–132 634, 2022.
  * [18] S. S. Nayyer, G. Revati, S. Wagh, and N. Singh, “Passivity and immersion based-modified gradient estimator: A control perspective in parameter estimation,” _arXiv preprint arXiv:2211.10674_ , 2022.
  * [19] P. M. Wensing and J.-J. Slotine, “Beyond convexity—contraction and global convergence of gradient descent,” _Plos one_ , vol. 15, no. 8, p. e0236661, 2020.
  * [20] A. Astolfi, D. Karagiannis, and R. Ortega, _Nonlinear and adaptive control with applications_. Springer, 2008, vol. 187.
  * [21] S. S. Nayyer, S. R. Wagh, and N. M. Singh, “Towards a constructive framework for stabilization and control of nonlinear systems: Passivity and immersion (p$\backslash$&i) approach,” _arXiv preprint arXiv:2208.10539_ , 2022.
  * [22] F. Alvarez, H. Attouch, J. Bolte, and P. Redont, “A second-order gradient-like dissipative dynamical system with hessian-driven damping.: Application to optimization and mechanics,” _Journal de mathématiques pures et appliquées_ , vol. 81, no. 8, pp. 747–779, 2002.
  * [23] V. Talasila, J. Clemente-Gallardo, and A. Astolfi, “Normally hyperbolic controlled-invariant manifolds,” in _17th International Symposium on Mathematical Theory of Networks and Systems_ , 2006.
  * [24] N. Fenichel and J. Moser, “Persistence and smoothness of invariant manifolds for flows,” _Indiana University Mathematics Journal_ , vol. 21, no. 3, pp. 193–226, 1971.
  * [25] M. Hirsch, C. Pugh, and M. Shub, “Invariant manifolds (lecture notes in mathematics, 583),” 1977.
  * [26] M. J. Capiński and H. Kubica, “Persistence of normally hyperbolic invariant manifolds in the absence of rate conditions,” _Nonlinearity_ , vol. 33, no. 9, p. 4967, 2020.
  * [27] L. Wang, F. Forni, R. Ortega, and H. Su, “Immersion and invariance stabilization of nonlinear systems: A horizontal contraction approach,” in _2015 54th IEEE Conference on Decision and Control (CDC)_. IEEE, 2015, pp. 3093–3097.
  * [28] Y. Nesterov, _Introductory lectures on convex optimization: A basic course_. Springer Science & Business Media, 2003, vol. 87.
  * [29] B. Sun, J. George, and S. Kia, “High-resolution modeling of the fastest first-order optimization method for strongly convex functions,” in _2020 59th IEEE Conference on Decision and Control (CDC)_. IEEE, 2020, pp. 4237–4242.


# Quantum droplet of a two-component Bose gas in an optical lattice near the
Mott insulator transition

Yoshihiro Machida1, Ippei Danshita1, Daisuke Yamamoto2, Kenichi Kasamatsu1
1Department of Physics, Kindai University, Higashi-Osaka, Osaka 577-8502,
Japan
2Department of Physics, Nihon University, Sakurajosui, Setagaya, Tokyo
156-8550, Japan

###### Abstract

We theoretically study dynamical formation of a quantum droplet in a two-
component Bose-Hubbard system with an external trap potential. Specifically,
the superfluid in the central region surrounded by the Mott insulator with
double filling forms a quantum droplet, which is self-bound thanks to the
discontinuous nature of the quantum phase transition between the two phases.
We show how to induce the characteristic behavior of the droplet through the
control of the trap potential by using the time-dependent Gutzwiller
simulations in a two-dimensional system. The static and dynamical properties
of the droplet can be described qualitatively by the effective Ginzburg-Landau
field theory with cubic-quintic nonlinearities, where the attractive cubic
nonlinearlity emerges although all the bare interparticle interactions are
repulsive.

In the past few years, there has been tremendous interest in studying the
properties of quantum droplets in ultracold gases [1, 2]. The droplet state is
sustained by the energy competition between the well-tuned mean-field
interactions and the intrinsic quantum fluctuations [3]. Experimentally, a lot
of efforts have been underway to reveal the properties of the quantum droplets
[4, 5, 6, 7, 5, 8]. While the previous experimental studies on quantum
droplets have addressed only systems in continuum, recent theoretical studies
have found droplet phases in a lattice system, which is specifically two-
component Bose gases in an optical lattice (OL) with the repulsive
intracomponent and attractive intercomponent interactions [9, 10]. Ultracold
atomic gases in an OL are nowadays standard platforms to implement versatile
quantum simulators of quantum many-body systems [11, 12]. The controllability
of the system parameters from a weakly-interacting regime to a strongly-
interacting one enables us to explore various states of matter.

In this work, we propose an alternative mechanism to realize a quantum droplet
in a two-component Bose-Hubbard (BH) system with both intra- and inter-
component interactions are repulsive. The Bose-Bose mixture in an OL is
experimentally accessible [13, 14]. The previous theoretical studies have
considered the ground state phase diagrams of the two-component BH model,
predicting the existence of several phases and phase transitions; especially,
there is the discontinuous (first-order) phase transition between the Mott
insulator (MI) and the superfluid (SF) phase near the tip of the Mott-lobe in
the spatial $d$-dimension with $d\geq 2$ [15, 16, 17, 18, 19, 20]. Also,
nonequilibrium dynamics subject to the two-component BH model has been studied
to reveal the impact of the inter-atomic interactions on dynamical behavior
[21, 22, 23, 24, 25].

Our dynamical simulations based on the Gutzwiller approximation show that the
two-component Bose gas in a combined OL and trap potential can accommodate a
localized structure of the SF component, provided by the discontinuous nature
of the SF-MI transition. It is well known that when the interaction dominates
over the kinetic energy, the density distribution of BH systems with a
parabolic trap forms an inhomogeneous wedding-cake structure [26]. We find
that when the transition between the central SF and the surrounding MI is of
the first-order, the central SF does not expand even after the trap around it
is locally turned off, behaving as a self-bound droplet. We explain the
mechanism of the self-binding from a viewpoint of the effective Ginzburg-
Landau (GL) theory, where the SF order parameter around the transition point
can be described qualitatively by the effective GL equation with the
attractive-cubic and repulsive-quintic nonlinearity. This type of nonlinear
equation can support the solution of a self-bound droplet [27, 28, 29]. Since
the effective attraction is induced by the fourth-order perturbation process
of the hoppings from the atomic limit, the droplet formation originates purely
from quantum nature. Using the effective theory, we study stationary density
profiles and collective excitations of the droplet, which are key features in
the experimental observation.

We consider the BH Hamiltonian for the two-component Bose gases in a two-
dimensional square OL and a trap potential, which is given by

$\displaystyle\hat{H}=$
$\displaystyle\sum_{\alpha=1,2}\biggl{[}-J_{\alpha}\sum_{\braket{i,j}}(\hat{a}^{\dagger}_{\alpha,j}\hat{a}_{\alpha,i}+\hat{a}^{\dagger}_{\alpha,i}\hat{a}_{\alpha,j})-\sum_{j}\mu_{\alpha,j}\hat{n}_{\alpha,j}$
$\displaystyle+\frac{U_{\alpha}}{2}\sum_{j}\hat{n}_{\alpha,j}(\hat{n}_{\alpha,j}-1)\biggl{]}+U_{12}\sum_{j}\hat{n}_{1,j}\hat{n}_{2,j}$
(1)

with the component index $\alpha=1,2$ and the site index $j=(j_{x},j_{y})$.
Here, $\braket{i,j}$ represents the nearest neighbor sites, $J_{\alpha}$ the
hopping coefficient, $\mu_{\alpha,j}$ the local chemical potential,
$U_{\alpha}$ and $U_{12}$ the intracomponent on-site interaction and the
intercomponent one, respectively. The operators $\hat{a}_{\alpha,j}$
$(\hat{a}_{\alpha,j}^{\dagger})$ and
$\hat{n}_{\alpha,j}\equiv\hat{a}^{\dagger}_{\alpha,j}\hat{a}_{\alpha,j}$ are
the annihilation (creation) operator and the number operator at site $j$,
respectively.

We assume below the symmetric choice of the parameters, namely
$J_{1}=J_{2}\equiv J$, $\mu_{j,1}=\mu_{j,2}\equiv\mu_{j}$, and
$U_{1}=U_{2}\equiv U>0$. For example, when we choose the two hyperfine states
of ${}^{87}{\rm Rb}$ atom, $|F=2,m_{F}=-1\rangle$ and $|1,1\rangle$, as the
two components, in a standard optical lattice the equal hoppings hold and the
equal intracomponent interactions approximately do ($U_{1}/U_{2}\simeq 0.95$
[30]. By assuming further that the trap potential is equal for the two
components, the equal local chemical potentials approximately imply the equal
numbers of particles. The parameter $U_{12}$ can be controlled with use of a
magnetic Feshbach resonance [31].

We calculate the ground state and time-dependent dynamics of the system
obeying Eq. (7) through the Gutzwiller approximation for the many-body wave
function [26, 17, 21, 24]. The Gutzwiller variational wave function is given
by

$\ket{\Psi_{\text{G}}(t)}=\prod_{j}\sum_{n_{1},n_{2}=0}f_{n_{1},n_{2}}^{(j)}(t)\ket{n_{1},n_{2}}_{j}$
(2)

where $\ket{n_{1},n_{2}}_{j}$ represents the Fock state associated with the
particle numbers for the two components at site $j$. The probability amplitude
$f_{n_{1},n_{2}}^{(j)}$ satisfies the normalization condition
$\sum_{n_{1},n_{2}}|f_{n_{1},n_{2}}^{(j)}|^{2}=1$. Under the variational
principle, the minimization of
$\bra{\Psi_{\text{G}}}\hat{H}-i\hbar\frac{d}{dt}\ket{\Psi_{\text{G}}}$ with
respect to $f_{n_{1},n_{2}}^{(j)*}$ gives the time-dependent equations for
$f_{n_{1},n_{2}}^{(j)}$ as

$\displaystyle
i\hbar\frac{df_{n_{1},n_{2}}^{(j)}}{dt}=\sum_{\alpha=1,2}\biggl{[}\frac{U}{2}n_{\alpha,j}(n_{\alpha,j}-1)-\mu_{j}n_{\alpha,j}\biggl{]}f_{n_{1},n_{2}}^{(j)}$
$\displaystyle-J\sum_{i\in
N(j)}\biggl{(}\Psi_{1,i}\sqrt{n_{1,j}}f_{n_{1}-1,n_{2}}^{(j)}+\Psi_{1,i}^{*}\sqrt{n_{1,j}+1}f_{n_{1}+1,n_{2}}^{(j)}\biggl{)}$
$\displaystyle-J\sum_{i\in
N(j)}\biggl{(}\Psi_{2,i}\sqrt{n_{2,j}}f_{n_{1},n_{2}-1}^{(j)}+\Psi_{2,i}^{*}\sqrt{n_{2,j}+1}f_{n_{1},n_{2}+1}^{(j)}\biggl{)}$
$\displaystyle+U_{12}n_{1,j}n_{2,j}f_{n_{1},n_{2}}^{(j)}.$ (3)

Here, $N(j)$ represents the nearest neighbor sites of site $j$. The SF order
parameter for each component is written by
$\Psi_{1,j}=\sum_{n_{1},n_{2}}f_{n_{1}-1,n_{2}}^{(j)*}\sqrt{n_{1}}f_{n_{1},n_{2}}^{(j)}$
and
$\Psi_{2,j}=\sum_{n_{1},n_{2}}f_{n_{1},n_{2}-1}^{(j)*}\sqrt{n_{2}}f_{n_{1},n_{2}}^{(j)}$.
In the following, real-time dynamics are simulated by solving Eq. (3) with the
Crank-Nicholson scheme, while the initial ground state is obtained by the
imaginary time propagation of Eq. (3).

Figure 1: The panel (a) shows the ground state phase diagram in the
$zJ/U$-$\mu/U$ plane obtained by the Gutzwiller analysis of the BH model Eq.
(7) [17, 18]. The solid and dashed curves between SF and MI represent the
first and second order transitions, respectively. The vertical dashed lines
show the radial dependence of the equilibrium phase in the bottom panels of
(d), determined by the local chemical potential. The panel (b) shows the
equilibrium condensate density $\rho=\sum_{\alpha}|\Psi_{\alpha}|^{2}$ for
several values of $zJ/U$ as a function of the chemical potential $\mu$
measured from the transition point $\mu_{c}$. For $\mu<\mu_{c}$ the SF state
is metastable. In (c), we show how to change the trap potential for inducing
characteristic behavior of a droplet formed by the SF in the central region of
the trap. The panels (d) show the time evolution of the cross section of the
SF density along the $x$-axis starting from the initial states of the bottom
panels, in which the profiles of the total particle density $n_{1}+n_{2}$ (red
solid curve) and the SF density $\rho$ (blue dashed curve) are shown. The
parameters are $U_{12}/U=0.9$, $\tilde{k}=0.001$, and
$(zJ/U,\mu/U,R)=(0.14,1.56,21)$ for (d-1) and $(0.05,1.85,35)$ for (d-2).
Here, $\mu$ is chosen as $\mu=\mu_{c}+0.05U$ The vertical dashed-dotted lines
show the radius $j_{x}=\pm R$ of the flat potential area. The panel (e) shows
the time evolution of the condensate density at the center for several values
of $zJ/U$. The horizontal thin lines show the gap $\Delta$ of the condensate
density at the transition points for the corresponding $zJ/U$, obtained in
(b).

The mean-field phase diagram of the ground state in the homogeneous system is
obtained in Refs. [17, 18], showing rich phase structures depending on the
values of $U_{12}/U$. In this work, we confine ourselves to $U_{12}/U=0.9$;
the phase diagram in the $(zJ/U)$-$(\mu/U)$ plane is shown in Fig. 1(a), where
$z=4$ is the coordinate number. The mean-field result is in qualitative
agreement with the quantum Monte-Carlo result [19]. There appear the two
tricritical points (TCPs) for $0.68<U_{12}/U<1$ upon the boundary between the
SF and the MI with double filling factors. In Fig. 1(b), the condensate
density $\rho=\sum_{\alpha}|\Psi_{\alpha}|^{2}$ obtained from Eq. (3) is
plotted for several values of $zJ/U$ as a function of the chemical potential
measured from the transition point, which is determined at the crossing point
of the energy $\langle\Psi_{\text{G}}|\hat{H}|\Psi_{\text{G}}\rangle$ of the
MI and that of the SF. At the transition point, the condensate density at the
ground state shows a jump, denoted as $\Delta$, a clear signature of the
first-order transition.

Next, we introduce the trap potential $V_{j}$. In the local density
approximation, the chemical potential consists of the global value $\mu$ and
the parabolic trap $V_{j}$ as $\mu_{j}=\mu-V_{j}=\mu-k(a_{0}j)^{2}/2$ with the
spring constant $k$, the lattice constant $a_{0}$, and
$j^{2}=j_{x}^{2}+j_{y}^{2}$. Then, the equilibrium number distribution
constitutes a wedding-cake structure depending on the global value of $\mu$.
We choose $\mu/U$ such that the SF is positioned at the central region of the
trap and the MI with $n=2$ surrounds the SF. The Thomas-Fermi radius of the
density profile is given by $a_{0}|j|=\sqrt{2\mu/k}$. The bottom panel of Fig.
1(d-1) shows the equilibrium profiles of the particle number $n$ and the SF
density $\rho$ when the phase boundary between the central SF and the
surrounding MI with $n=2$ associates with the first-order transition. In
contrast, that of Fig. 1(d-2) shows the case when the phase boundary
corresponds to the second-order transition.

In the former case, the central SF can possess a characteristics of a droplet
in the sense that the localized structure is dynamically kept even after the
trap potential is turned off. To demonstrate this, we slowly change the
potential $V_{j}$ under the protocol shown in Fig. 1(c) so as to flatten the
area $j^{2}a_{0}^{2}\leq R^{2}$. Here, the radius $R$ of the flat area is
determined at the position sufficiently inside the $n=2$ MI domain. The time
sequence of the trap potential for $j^{2}a_{0}^{2}\leq R^{2}$ is set as
$V_{j}=\delta\epsilon(t)(V_{j}^{(0)}-\epsilon_{0})+\epsilon_{0}$, where
$V_{j}^{(0)}=k(a_{0}j)^{2}/2$ and we change the time-dependent parameter as
$\delta\epsilon(t)=1\to 0$ linearly within the interval $t=[0,T=1000\hbar/U]$.
In the flat area, for $t>T$ the chemical potential takes a constant
$\mu-\epsilon_{0}$, which is set to be the value inside the $n=2$ MI region.
This local manipulation of the trap is necessary since the SF droplet is
sustained only in the presence of the surrounding MI domain. Such local
control of the external potential can be achieved with recent experimental
techniques, e.g., using the degital micromirror device [32].

We calculate the real time dynamics of the SF component using the time-
dependent Gutzwiller approximation. The top panel of Fig. 1(d-1) represents
the time evolution of $\rho$ starting from the initial state of the bottom
panel ($zJ/U=0.14$), where we choose $R=21a_{0}$. Even though the confining
potential is flat, the central SF component does not expand. We also show the
time evolution of $\rho$ at the central site in Fig. 1(e). After
$Ut/\hbar=1000$, the central density makes an oscillation around a certain
steady value. This steady value is almost coincident with the jump $\Delta$ of
the condensate density at the first-order transition point. We confirmed that
this localized structure appears when the central SF is surrounded by the
$n=2$ MI through the first-order transition line. When the simulation starts
from the bottom panel of Fig. 1(d-2) ($zJ/U=0.05$), on the other hands, the SF
density expands and reaches to the edge of the flat region; the density makes
an interference pattern due to the reflection from the potential edge.

The appearance of the localized structure can be interpreted as the nature of
the phase transition. For the first-order transition, the condensate density
cannot continuously change down to zero at the SF-MI boundary. Thus, the
monotonic expansion of the central density is prohibited so that the mean
density cannot decrease below the gap $\Delta$. To understand these properties
more clearly, we employ the effective sixth-order GL action expanded by the SF
order parameter $\Psi_{\alpha}$ [19, 20]. Assuming the symmetric form
$\Psi_{1}=\Psi_{2}=\Psi$ for simplicity, we get the cubic-quintic GL equation

$i\hbar K\frac{\partial\Psi}{\partial
t}=\biggl{[}-\frac{\hbar^{2}}{2m}\nabla^{2}-r_{0}+u_{+}|\Psi|^{2}+w_{+}|\Psi|^{4}\biggl{]}\Psi.$
(4)

Since we confine ourselves to the low-energy dynamics of the quantum droplet,
we take only the first-order time derivative in the GL action. The effective
mass is given by $m=\hbar^{2}/(2Ja_{0}^{2})$. Although the GL parameters have
the spatial dependence since $\mu_{j}$ depends on the position in a trap, we
approximate them by the values upon the transition line between the $n=2$ MI
and the SF. The GL parameters $(K,r_{0},u_{+},w_{+})$ along the first-order
transition line is shown in Fig. 2(a). The negative $r_{0}$ ensures that the
free energy around the MI phase takes a minimum.

Figure 2: (a) The GL parameters along the first-order transition line between
the SF and the $n=2$ MI as a function of $\mu/U$ for $U_{12}/U=0.9$. The
solid, short dashed, and long dashed curves represent $u_{+}/(a_{0}^{d}U)$,
$w_{+}/(a_{0}^{2d}U)$, and $K$, respectively (Here, $d=2$). The densely dotted
curve represents $r_{0}/U$, whose vertical axis is the right side. (b) One of
the fourth-order processes that contribute to the negativity of $u_{+}$.

According to the GL theory, the emergence of the first-order transition
requires $u_{+}<0$, which is satisfied in a regime where the central SF forms
the droplet, as seen in Fig. 2(a). The negativity of $u_{+}$ can be understood
by the fourth-order perturbation processes [33], in which the hopping terms
are regarded as the perturbation. One of the perturbation processes is shown
in Fig. 2(b). When $U\sim U_{12}$, the energy of $|1,1\rangle$ is nearly equal
to that of the intermediate state $|2,0\rangle$ or $|0,2\rangle$. Then, these
processes give rise to an enhanced negative contribution of the coupling
constant in front of $|\Psi_{1}|^{2}|\Psi_{2}|^{2}$, which results in the
attractive cubic-nonlinearlity in Eq. (4). Indeed, the first-order transition
emerges only when $U\sim U_{12}$ (more precisely when $0.68<U/U_{12}<1$)
according to the Gutzwiller analysis [17, 18].

The negative $u_{+}$ also provides a clear evidence of the existence of a
self-bound droplet. First, let us see the stationary solution of Eq. (4) by
assuming the isotropy of the solution $\Psi(\bm{r})=\psi(r)$. The solutions
can be obtained through the imaginary time propagation of Eq. (4) with the
fixed norm $\int d^{2}r|\psi|^{2}\equiv N_{\rho}$, where $r_{0}$ reads the
chemical potential of the condensate. For $u_{+}<0$ we can get localized
solutions even in a free space [27, 28, 29]. With increasing $N_{\rho}$, the
solution exhibits a “flat-top” shape, which is typical of the droplet
structure [2, 3]. When the droplet size is much larger than the length scale
$l_{w}\equiv\hbar/\sqrt{mw_{+}n_{0}^{2}}$, where $n_{0}$ is the density at the
flat region of the droplet, the derivative term $(1/r)\partial_{r}$ in the
Laplacian can be neglected. Then, the droplet solution with the radius
$R_{\text{d}}$, at which the density becomes a half of that of the center, can
be written by [20]

$\psi^{2}=\frac{n_{0}}{2}\left[1-\tanh\left(\sqrt{\frac{2}{3}}\frac{r-R_{\text{d}}}{l_{w}}\right)\right].$
(5)

Here, $n_{0}$ is given by the jump $\Delta_{\text{GL}}=-3u_{+}/(4w_{+})$ of
the SF density at the first-order transition point. The numerical solution of
Eq. (4) can be well approximated by Eq. (5), the central density being
slightly exceeded from the above estimation $2n_{0}$ [33].

Figure 3: Comparison between the Gutzwiller calculation and the GL analysis.
In (a), the red, blue, green curves represent the stationary profiles of the
condensate density $\rho=2|\psi|^{2}$ by solving Eq. (4) for the GL parameters
corresponding to $(zJ/U,\mu/U)=$ (0.14,1.51), (0.12,1.61), and (0.11.1.64),
respectively. The profiles obtained from the dynamical Gutzwiller simulation
are also shown by circles, triangles, and squares for the same parameters. The
panel (b) shows the eigenfrequencies of the collective excitations with the
quantum numbers $(n,l)=(1,0)$, obtained by solving the BdG equation for
$N_{\rho}(2U/zJ)\equiv N^{\prime}=100$ (long-dashed), 200 (short-dashed), 500
(dotted) and 1000 (solid). The circles represent the results by the Gutzwiller
simulations, where the frequency is extracted from the power spectrum of the
oscillation in Fig. 1(e) and the attached number represents $N^{\prime}$.

In Fig. 3(a), we plot the typical stationary solutions of Eq. (4) as well as
the profile of the condensate density obtained from the Gutwiller simulations
in Fig. 1, averaged for $3000\leq Ut/\hbar\leq 5000$. Here, the norm of $\psi$
is determined to reproduce the condensate number $N_{\rho}$ within the droplet
of $|\Psi_{\text{G}}\rangle$ [33]. The profiles of the Gutzwiller calculations
are well described by the GL results as the parameters approach to the TCP.
This is naturally understood that the GL expansion is quantitatively validated
for a smaller values of a jump of the SF order parameter at the transition
point. Nevertheless, the overall qualitative feature of the quantum droplet is
well captured by the GL solutions. For example, the central density of the
Gutzwiller calculation is determined by the density jump $\Delta$ shown in
Fig. 1(b), whose property has been clearly seen in Fig. 1(e).

The small oscillation seen in Fig. 1(e) is related with the collective modes
of the quantum droplet [34, 35, 36] induced by the temporal change of the
potential. We calculate the frequencies of the collective excitations by the
Bogoliubov-de Genne (BdG) analysis, in which we expand the order parameter as
$\Psi(\bm{r},t)=[\Psi(\bm{r})+\delta\Psi(\bm{r},t)]e^{-ir_{0}t/\hbar}$ with
$\delta\Psi(\bm{r},t)=\sum_{n}[u_{n}(r)e^{il\theta-i\omega
t}-v_{n}^{*}(r)e^{-il\theta+i\omega t}]$ and solve the eigenvalue equations.
Figure 3 (b) shows the excitation frequencies with $(n,l)=(1,0)$, i.e.,
breathing mode, along the first-order transition line as a function of $zJ/U$.
The frequency is decreased from the tip of the Mott lobe ($zJ/U\simeq 0.155$)
to the TCP ($zJ/U\simeq 0.096$), being insensitive to the norm for
$N_{\rho}(2U/zJ)>200$. We also plot the oscillation frequencies, obtained from
the main peak of the power spectrum, of the central density in the Gutzwiller
simulations [Fig. 1(e)]. The plots are reasonably coincident with the
frequencies of the BdG analysis, since excitations with $l\neq 0$ cannot be
induced in our potential protocol. We note that since a number of modes with
$n\geq 2$ are lying above the $n=1$ mode with small energy gaps especially
near the TCP, multiple modes are excited simultaneously. Although the GL
approximation should become better near the TCP, it is difficult to make a
quantitative comparison of the BdG analysis since a pure breathing oscillation
with a small amplitude cannot occur there.

In summary, we have shown that a suitable manipulation of an external trap
leads to the formation of a quantum droplet of the SF phase surrounded by the
$n=2$ MI in a two-component BH system. The underlying mechanism of the droplet
formation is the effective attraction between the intercomponent SF order
parameters, which originates from the fact that the fourth-order process of
the perturbative expansion is enhanced for $U_{12}\sim U$, even though the on-
site interactions are all repulsive. We showed through the time-dependent
Gutzwiller analysis that the localized structure of the SF phase can be kept
even after the external trap is locally flattened. The static and dynamical
properties of the droplet can be described qualitatively by the effective GL
theory with the cubic-quintic nonlinearity.

This work was financially supported by JSPS KAKENHI Grants Nos. JP18K03472
(K.K.), JP18H05228 (I.D.), JP21H01014 (I.D.), JP18K03525 (D.Y.) and JP21H05185
(D.Y.), by JST CREST Grant No. JPMJCR1673 (I.D. and D.Y.), by JST FOREST Grant
No. JPMJFR202T (I.D.), and by MEXT Q-LEAP Grant No. JPMXS0118069021 (I.D.).

## References

  * Ferrier-Barbut [2019] I. Ferrier-Barbut, Ultradilute quantum droplets, Physics Today 72, 46 (2019).
  * Luo  _et al._ [2021] Z.-H. Luo, W. Pang, B. Liu, Y.-Y. Li, and B. A. Malomed, A new form of liquid matter: Quantum droplets, Frontiers of Physics 16, 1 (2021).
  * Petrov [2015] D. S. Petrov, Quantum mechanical stabilization of a collapsing bose-bose mixture, Phys. Rev. Lett. 115, 155302 (2015).
  * Ferrier-Barbut  _et al._ [2016] I. Ferrier-Barbut, H. Kadau, M. Schmitt, M. Wenzel, and T. Pfau, Observation of quantum droplets in a strongly dipolar bose gas, Phys. Rev. Lett. 116, 215301 (2016).
  * Chomaz  _et al._ [2016] L. Chomaz, S. Baier, D. Petter, M. J. Mark, F. Wächtler, L. Santos, and F. Ferlaino, Quantum-fluctuation-driven crossover from a dilute bose-einstein condensate to a macrodroplet in a dipolar quantum fluid, Phys. Rev. X 6, 041039 (2016).
  * Semeghini  _et al._ [2018] G. Semeghini, G. Ferioli, L. Masi, C. Mazzinghi, L. Wolswijk, F. Minardi, M. Modugno, G. Modugno, M. Inguscio, and M. Fattori, Self-bound quantum droplets of atomic mixtures in free space, Phys. Rev. Lett. 120, 235301 (2018).
  * Cabrera  _et al._ [2018] C. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarruell, Quantum liquid droplets in a mixture of bose-einstein condensates, Science 359, 301 (2018).
  * Ferioli  _et al._ [2019] G. Ferioli, G. Semeghini, L. Masi, G. Giusti, G. Modugno, M. Inguscio, A. Gallemí, A. Recati, and M. Fattori, Collisions of self-bound quantum droplets, Physical review letters 122, 090401 (2019).
  * Morera  _et al._ [2020] I. Morera, G. E. Astrakharchik, A. Polls, and B. Juliá-Díaz, Quantum droplets of bosonic mixtures in a one-dimensional optical lattice, Physical Review Research 2, 022008 (2020).
  * Morera  _et al._ [2021] I. Morera, G. E. Astrakharchik, A. Polls, and B. Juliá-Díaz, Universal dimerized quantum droplets in a one-dimensional lattice, Physical Review Letters 126, 023001 (2021).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Reviews of modern physics 80, 885 (2008).
  * Lewenstein  _et al._ [2012] M. Lewenstein, A. Sanpera, and V. Ahufinger,  _Ultracold Atoms in Optical Lattices: Simulating quantum many-body systems_ (Oxford University Press, 2012).
  * Catani  _et al._ [2008] J. Catani, L. De Sarlo, G. Barontini, F. Minardi, and M. Inguscio, Degenerate bose-bose mixture in a three-dimensional optical lattice, Physical Review A 77, 011603 (2008).
  * Gadway  _et al._ [2010] B. Gadway, D. Pertot, R. Reimann, and D. Schneble, Superfluidity of interacting bosonic mixtures in optical lattices, Physical review letters 105, 045303 (2010).
  * Kuklov  _et al._ [2004] A. Kuklov, N. Prokof’ev, and B. Svistunov, Commensurate two-component bosons in an optical lattice: Ground state phase diagram, Physical review letters 92, 050402 (2004).
  * Chen and Yang [2010] P. Chen and M.-F. Yang, Quantum phase transitions in a two-species hard-core boson hubbard model in two dimensions, Physical Review B 82, 180510 (2010).
  * Ozaki  _et al._ [2012] T. Ozaki, I. Danshita, and T. Nikuni, Bose-bose mixtures in an optical lattice: First-order superfluid-insulator transition and elementary excitations, arXiv preprint arXiv:1210.1370 (2012).
  * Yamamoto  _et al._ [2013] D. Yamamoto, T. Ozaki, C. A. S. de Melo, and I. Danshita, First-order phase transition and anomalous hysteresis of bose gases in optical lattices, Physical Review A 88, 033624 (2013).
  * Kato  _et al._ [2014] Y. Kato, D. Yamamoto, and I. Danshita, Quantum tricriticality at the superfluid-insulator transition of binary bose mixtures, Physical Review Letters 112, 055301 (2014).
  * Danshita  _et al._ [2015] I. Danshita, D. Yamamoto, and Y. Kato, Cubic-quintic nonlinearity in superfluid bose-bose mixtures in optical lattices: Heavy solitary waves, barrier-induced criticality, and current-phase relations, Physical Review A 91, 013630 (2015).
  * Wernsdorfer  _et al._ [2010] J. Wernsdorfer, M. Snoek, and W. Hofstetter, Lattice-ramp-induced dynamics in an interacting bose-bose mixture, Physical Review A 81, 043620 (2010).
  * Hu  _et al._ [2011] A. Hu, L. Mathey, E. Tiesinga, I. Danshita, C. J. Williams, and C. W. Clark, Detecting paired and counterflow superfluidity via dipole oscillations, Physical Review A 84, 041609 (2011).
  * Zhang  _et al._ [2012] J. Zhang, C. Shen, and W. Liu, Strong thermalization of the two-component bose-hubbard model at finite temperatures, Physical Review A 85, 013637 (2012).
  * Lundh and Martikainen [2012] E. Lundh and J.-P. Martikainen, Kelvin-helmholtz instability in two-component bose gases on a lattice, Physical Review A 85, 023628 (2012).
  * Shrestha and Ruostekoski [2012] U. Shrestha and J. Ruostekoski, Fragmentation, domain formation and atom number fluctuations of a two-species bose–einstein condensate in an optical lattice, New Journal of Physics 14, 043037 (2012).
  * Jaksch  _et al._ [1998] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Cold bosonic atoms in optical lattices, Physical Review Letters 81, 3108 (1998).
  * Pushkarov  _et al._ [1979] K. I. Pushkarov, D. Pushkarov, and I. Tomov, Self-action of light beams in nonlinear media: soliton solutions, Optical and Quantum Electronics 11, 471 (1979).
  * Pego and Warchall [2002] R. L. Pego and H. A. Warchall, Spectrally stable encapsulated vortices for nonlinear schrödinger equations., Journal of Nonlinear Science 12 (2002).
  * Xi and Saito [2016] K.-T. Xi and H. Saito, Droplet formation in a bose-einstein condensate with strong dipole-dipole interaction, Physical Review A 93, 011604 (2016).
  * Mertes  _et al._ [2007] K. Mertes, J. Merrill, R. Carretero-González, D. Frantzeskakis, P. Kevrekidis, and D. Hall, Nonequilibrium dynamics and superfluid ring excitations in binary bose-einstein condensates, Physical review letters 99, 190402 (2007).
  * Tojo  _et al._ [2010] S. Tojo, Y. Taguchi, Y. Masuyama, T. Hayashi, H. Saito, and T. Hirano, Controlling phase separation of binary bose-einstein condensates via mixed-spin-channel feshbach resonance, Physical Review A 82, 033609 (2010).
  * Mazurenko  _et al._ [2017] A. Mazurenko, C. S. Chiu, G. Ji, M. F. Parsons, M. Kanász-Nagy, R. Schmidt, F. Grusdt, E. Demler, D. Greif, and M. Greiner, A cold-atom fermi–hubbard antiferromagnet, Nature 545, 462 (2017).
  * [33] See the supplemental material for the details.
  * Tylutki  _et al._ [2020] M. Tylutki, G. E. Astrakharchik, B. A. Malomed, and D. S. Petrov, Collective excitations of a one-dimensional quantum droplet, Physical Review A 101, 051601 (2020).
  * Hu and Liu [2020] H. Hu and X.-J. Liu, Collective excitations of a spherical ultradilute quantum droplet, Physical Review A 102, 053303 (2020).
  * Stürmer  _et al._ [2021] P. Stürmer, M. N. Tengstrand, R. Sachdeva, and S. M. Reimann, Breathing mode in two-dimensional binary self-bound bose-gas droplets, Physical Review A 103, 053302 (2021).

Supplemental Material for “Quantum droplet of a two-component Bose gas in an
optical lattice near the Mott insulator transition”

##  .1 The GL action and the coupling constant $u_{12}$

The derivation of the effective sixth-order GL action is described in Ref.
[19], the form of the action being

$\displaystyle S^{\rm{eff}}=\int dt\int
d^{d}x\Biggl{[}\sum_{\alpha=1,2}\biggl{(}$ $\displaystyle i\hbar
K\Psi_{\alpha}^{*}\frac{\partial\Psi_{\alpha}}{\partial
t}-\hbar^{2}W\biggl{|}\frac{\partial\Psi_{\alpha}}{\partial
t}\biggl{|}^{2}-\frac{\hbar^{2}}{2m}|\nabla\Psi_{\alpha}|^{2}+r_{0}|\Psi_{\alpha}|^{2}-\frac{u}{2}|\Psi_{\alpha}|^{4}-\frac{w}{3}|\Psi_{\alpha}|^{6}\biggl{)}$
$\displaystyle-
u_{12}|\Psi_{1}|^{2}|\Psi_{2}|^{2}-w_{12}\left(|\Psi_{1}|^{4}|\Psi_{2}|^{2}+|\Psi_{1}|^{2}|\Psi_{2}|^{4}\right)\Biggl{]}.$
(6)

The GL parameters are given as a function of the parameters
$(zJ/U,\>U_{12}/U,\>\mu_{\alpha}/U)$ of the BH Hamiltonian

$\displaystyle\hat{H}=\sum_{\alpha=1,2}\biggl{[}-J\sum_{\braket{i,j}}(\hat{a}^{\dagger}_{\alpha,j}\hat{a}_{\alpha,i}+\hat{a}^{\dagger}_{\alpha,i}\hat{a}_{\alpha,j})-\mu_{\alpha}\sum_{j}\hat{n}_{\alpha,j}+\frac{U}{2}\sum_{j}\hat{n}_{\alpha,j}(\hat{n}_{\alpha,j}-1)\biggl{]}+U_{12}\sum_{j}\hat{n}_{1,j}\hat{n}_{2,j}$
(7)

through a perturbative expansion, the explicit form being given in Ref. [19].
Here, we give only the expression of the expansion parameter $u_{12}$, which
is important in our study since its negativity realizes the first-order
transition.

Let us treat the hopping term in the BH Hamiltonian as a perturbation. The
ground state of the unperturbed Hamiltonian is given by a product of local
Fock state $|n_{1},n_{2}\rangle=|g,g\rangle$ and the unperturbed energy is

$E_{n_{1},n_{2}}=\sum_{\alpha}\left[-\mu_{\alpha}n_{\alpha}+\frac{U}{2}n_{\alpha}(n_{\alpha}-1)\right]+U_{12}n_{1}n_{2}.$
(8)

In our work, we have focused on $g=1$, namely $n=n_{1}+n_{2}=2$. Using this,
the fourth-order perturbation analysis yields the expansion parameter $u_{12}$
of the GL action as

$\displaystyle u_{12}=a_{0}^{d}Z^{4}J^{4}$
$\displaystyle\Biggl{[}\left(\frac{g+1}{E_{g+1,g}-E_{g,g}}+\frac{g}{E_{g-1,g}-E_{g,g}}\right)\left(\frac{g+1}{(E_{g,g+1}-E_{g,g})^{2}}+\frac{g}{(E_{g,g-1}-E_{g,g})^{2}}\right)$
$\displaystyle+\left(\frac{g+1}{E_{g,g+1}-E_{g,g}}+\frac{g}{E_{g,g-1}-E_{g,g}}\right)\left(\frac{g+1}{(E_{g+1,g}-E_{g,g})^{2}}+\frac{g}{(E_{g-1,g}-E_{g,g})^{2}}\right)$
$\displaystyle-\left(\frac{1}{E_{g+1,g}-E_{g,g}}+\frac{1}{E_{g,g+1}-E_{g,g}}\right)^{2}\frac{(g+1)^{2}}{E_{g+1,g+1}-E_{g,g}}$
$\displaystyle-\left(\frac{1}{E_{g+1,g}-E_{g,g}}+\frac{1}{E_{g,g-1}-E_{g,g}}\right)^{2}\frac{g(g+1)}{E_{g+1,g-1}-E_{g,g}}$
$\displaystyle-\left(\frac{1}{E_{g-1,g}-E_{g,g}}+\frac{1}{E_{g,g+1}-E_{g,g}}\right)^{2}\frac{g(g+1)}{E_{g-1,g+1}-E_{g,g}}$
$\displaystyle-\left(\frac{1}{E_{g-1,g}-E_{g,g}}+\frac{1}{E_{g,g-1}-E_{g,g}}\right)^{2}\frac{g^{2}}{E_{g-1,g-1}-E_{g,g}}\Biggr{]}.$
(9)

When $U\sim U_{12}$, the energy differences $E_{g+1,g-1}-E_{g,g}$ and
$E_{g-1,g+1}-E_{g,g}$ in the fourth and fifth lines in Eq. (9) become small,
so that these negative contributions are resonantly enhanced.

##  .2 The GL analysis in a uniform system

The GL equation

$i\hbar K\frac{\partial\Psi}{\partial
t}=\biggl{[}-\frac{\hbar^{2}}{2m}\nabla^{2}-r_{0}+u_{+}|\Psi|^{2}+w_{+}|\Psi|^{4}\biggl{]}\Psi.$
(10)

is useful to understand the SF-MI transition of the two-component BH model as
discussed thoroughly in Refs.[19, 20]. Here, we summarize several simple
results obtained from the stationary GL equation in a uniform system to
understand the present work.

First, let us consider a uniform solution $\Psi(\bm{r})=\sqrt{n_{0}}$ of the
GL equation, where assume $\Psi_{1}=\Psi_{2}=\Psi$. We also suppose $w_{+}>0$
below to ensure the stability of the solution. Equation (10) then reads

$\left(-r_{0}+u_{+}n_{0}+w_{+}n_{0}^{2}\right)\sqrt{n_{0}}=0,$ (11)

and yields three solutions

$\displaystyle
n_{0}=0,\quad\frac{-u_{+}+\sqrt{D}}{2w_{+}},\quad\frac{-u_{+}-\sqrt{D}}{2w_{+}}$
(12)

with $D=u_{+}^{2}+4r_{0}w_{+}$.

When $u_{+}\geq 0$, the first two solutions of Eq. (12) can satisfy $n_{0}\geq
0$. For $r_{0}<0$, the physical solution is only $n_{0}=0$, which means that
the ground state is the MI. For $r_{0}>0$, $n_{0}$ can have a positive finite
value and the SF becomes the ground state. Thus, the transition takes place at
$r_{0}=0$ and is of the second order type.

Figure 4: The three solutions of Eq. (12) as a function of $r_{0}$, where the
dotted line, the solid red curve, and the dashed blue curve represent the 1st,
2nd and 3rd solutions, respectively.

When $u_{+}<0$, all solutions of Eq. (12) can satisfy $n_{0}\geq 0$. The three
solutions are depicted in Fig. (4) as a function of $r_{0}$. The second
solution becomes a positive finite value when $D>0$, namely $r_{0}\geq-
u_{+}^{2}/(4w_{+})\equiv r_{0}^{\ast}$, being stable branch satisfying
$\partial n_{0}/\partial r_{0}>0$. The last solution becomes positive finite
when $r_{0}^{\ast}\leq r_{0}\leq 0$, corresponding to an unstable branch with
the relation $\partial n_{0}/\partial r_{0}<0$. The first-order transition
point can be determined by means of Maxwell’s construction as
$r_{0}=-3u_{+}^{2}/(16w_{+})\equiv r_{c}$. The MI state with $n_{0}=0$ is the
ground state when $r_{0}<r_{c}$, a metastable when $r_{c}<r_{0}<0$, and an
energetically unstable when $r_{0}>0$. The SF state with $n_{0}\neq 0$ is the
ground state when $r_{0}>r_{c}$ and a metastable state when
$r_{0}^{\ast}<r_{0}<r_{c}$. Substituting $r_{0}=r_{c}$ into the second
solution of Eq. (12), we evaluate the jump $\Delta_{\text{GL}}$ of $n_{0}$ at
the transition point as $\Delta_{\text{GL}}=-3u_{+}/(4w_{+})$.

##  .3 The dimensionless form of the GL equation

In order to scale Eq. (10) into the dimensionless form, we replace the order
parameter as $\Psi\to\tilde{\Psi}/a_{0}$ with the lattice constant $a_{0}$.
Dividing the both sides of Eq. (10) by $U$ of the BH Hamiltonian, we have

$i\hbar\frac{K}{U}\frac{\partial\tilde{\Psi}}{\partial
t}=\biggl{[}-\frac{\hbar^{2}}{2mU}\nabla^{2}-\frac{r_{0}}{U}+\frac{u_{+}}{Ua_{0}^{2}}|\tilde{\Psi}|^{2}+\frac{w_{+}}{Ua_{0}^{4}}|\tilde{\Psi}|^{4}\biggl{]}\tilde{\Psi}..$
(13)

Taking the characteristic time and length scale as $\tau=\hbar K/U$ and
$\xi=\hbar/\sqrt{mU}$, we obtain the dimensionless form of the GL equation

$i\frac{\partial\tilde{\Psi}}{\partial\tilde{t}}=\biggl{[}-\frac{1}{2}\tilde{\nabla}^{2}-\tilde{r}_{0}+\tilde{u}_{+}|\tilde{\Psi}|^{2}+\tilde{w}_{+}|\tilde{\Psi}|^{4}\biggl{]}\tilde{\Psi}..$
(14)

Here, the normalization of the order parameter $N_{\rho}=\int
d^{2}r|\Psi|^{2}$ in the dimensionless form is given by

$\int d^{2}\tilde{r}|\tilde{\Psi}|^{2}=N_{\rho}\frac{a_{0}^{2}}{\xi^{2}}\equiv
N^{\prime}.$ (15)

The factor $a_{0}/\xi$ is written as $a_{0}/\xi=\sqrt{2U/(zJ)}$. We
numerically solve Eq. (14) through the imaginary time propagation by fixing
the norm $N^{\prime}$ by assuming a radial isotropy $\Psi(\bm{r})=\psi(r)$. In
Fig. 5, we show the solutions using the GL parameters
$(u_{+},w_{+})=(-0.151,0.620)$ at $(zJ/U,\mu/U)=(0.14,1.51)$ for some typical
values of $N^{\prime}$.

Figure 5: The Condensate density $\rho=2|\psi(r)|^{2}$ obtained from the
solution of the GL equation for the GL parameters at
$(zJ/U,\mu/U)=(0.14,1.51)$ and $N^{\prime}=100$, 200, 400, 500, 750, 1000, and
2000. The solid curves are numerical solutions, while the dashed curves are
analytical solutions Eq. (5) in the text.

The Gaussian-like profile for a small $N^{\prime}$ transforms to the flat-top
shape with increasing $N^{\prime}$. We also plot the approximate analytic
solution given by Eq.(5) in the text, which can be obtained by neglecting the
first-order derivative $\partial_{r}/r$ in the Laplacian. Here, the radius
$R_{\text{d}}$ is determined from the norm $N^{\prime}$. The approximate
solutions are in reasonably good agreement with the numerical ones, but the
central density of the numerical solution is slightly exceeded from that of
the analytic solution; the latter is exactly given by $2\Delta_{\text{GL}}$.