# Dynamical representations of constrained multicomponent nonlinear
Schrödinger equations in arbitrary dimensions

M. Gulliksson* and M. Ögren

###### Abstract

We present new approaches for solving constrained multicomponent nonlinear
Schrödinger equations in arbitrary dimensions. The idea is to introduce an
artificial time and solve an extended damped second order dynamic system whose
stationary solution is the solution to the time-independent nonlinear
Schrödinger equation. Constraints are often considered by projection onto the
constraint set, here we include them explicitly into the dynamical system. We
show the applicability and efficiency of the methods on examples of relevance
in modern physics applications.

*(Corresponding author) School of Science and Technology, Örebro University, 701 82 Örebro, Sweden

E-mail address<EMAIL_ADDRESS>
**School of Science and Technology, Örebro University, 701 82 Örebro, Sweden,
Hellenic Mediterranean University, P.O. Box 1939, GR-71004, Heraklion, Greece.

E-mail address<EMAIL_ADDRESS>
##  1 Introduction

We will describe and numerically test a family of methods for finding
stationary solutions to coupled multicomponent nonlinear Schrödinger equations
in arbitrary number of spatial dimensions using different representations of
damped dynamical systems.

The nonlinear Schrödinger equation (NLSE) [1] is widely used today as a model,
for example in nonlinear optics [2]; in super-conductivity, modeled with the
related Ginzburg-Landau equation [3]; in models for dark matter [4]; and the
NLSEs have also been used in the description of water waves [5], as well as
for modeling the occurrence of rogue surface waves at sea [6].

Here we present examples of constrained vector-NLSEs [7], often called coupled
Gross-Pitaevskii equations in the settings of a mean-field description of
rotating bosonic atoms in different components. This is currently an active
area of research in experiments on ultra-cold atomic gases [8, 9] in need of
tools for theoretical investigations.

We will limit the presentation to the most common form of NLSEs, i.e. with
cubic nonlinearities, but the nonlinear terms can be more general functionals
of the densities of the components, such is used in modeling, for example of
Tonks gases [10], or superfluid Fermi-Bose mixtures [11], to mention just a
few realizations of different nonlinear couplings. The readers can
straightforwardly modify the presented constraints and dynamic equations to
represent future models of interest.

There are other numerical ways to solve for the stationary solution of the
NLSE numerically. The so called imaginary time dependent SE is common [12].
This means that after discretization in space (e.g., by finite differences,
finite elements, or spectral decomposition) solve a first order damped time
dependent equation numerically, see [13]. Sometimes these methods are called
steepest descent methods When the NLSE is independent of time it is, after
discretization in space, equivalent to a finite dimensional minimization
problem, with nonlinear constraints. Such minimization problems can generally
be solved by a variety of numerical methods including gradient descent
methods, (Quasi-) Newton methods, machine learning techniques etc. [14, 15].

However, as demonstrated in this article, the extended second order damped
dynamical systems presented have many benefits, such as ease of implementation
and low computational complexity.

We begin in the next section with the necessary notation and presentation of
the problem formulated as a minimization problem with constraints. We there
derive the necessary conditions for a solution. In Sec. 3 we introduce the
damped dynamical system to be used to attain a stationary solution that solves
our original constrained minimization problem. This approach to obtain the
solution is an extension of the Dynamical Functional Particle Method (DFPM),
see [17], able to handle nonlinear constraints. Our version of DFPM involves
finding the Lagrange parameters and we derive the linear equations determining
those. We apply DFPM to the NLSE in one dimension in detail in Sec. 4 based on
the general results in the earlier section. We show the applicability,
efficiency, and stability of DFPM in numerical tests on both one- and two-
dimensional NLSEs in Sec. 5. Finally, we summarize and discuss our results in
Sec. 6.

##  2 The multicomponent nonlinear Schrödinger Equation

In this section we present the general $n$-component Nonlinear Schrödinger
equation (NLSE) and constraints. Further, we define the corresponding total
energy and Lagrange functional to be used in later sections.

Let $\mathcal{H}$ denote a Hilbert space of complex valued vector functions
$\vec{\psi}:\mathbb{R}^{m}\longrightarrow\mathbb{C}^{n}$, $m,n\in\mathbb{N}$.
We use the standard notation in physics for the position vector as ${\bf r}$.
The space $\mathcal{H}$ is equipped with the inner product

$\left\langle\vec{\psi}\middle|\vec{\phi}\right\rangle=\int\vec{\psi}^{\dagger}\,\vec{\phi}\,{\rm
d}{\bf r},$ (1)

and the norm

$\|\vec{\psi}\|=\sqrt{\int\vec{\psi}^{\dagger}\,\vec{\psi}\,{\rm d}{\bf r},}$
(2)

where $\dagger$ is transpose conjugation. We assume that the functions in
$\mathcal{H}$ have sufficient smoothness and satisfy either periodic boundary
conditions or $\lim_{\|{\bf r}\|\rightarrow\infty}\vec{\psi}({\bf r})=0$.

The total energy is defined as the convex functional
$\widetilde{E}:\,\mathcal{H}\longrightarrow\mathbb{R}$with

$\widetilde{E}(\vec{\psi})=\int\sum_{i=1}^{m}\nabla_{i}\vec{\psi}^{\dagger}W\nabla_{i}\vec{\psi}+\vec{\psi}^{\dagger}V(x)\vec{\psi}+\vec{\psi}^{\dagger}\Gamma(\vec{\psi})\vec{{\psi}}\,{\rm
d}{\bf r},$ (3)

where $W=\text{diag}(1/(2M_{1}),\ldots,1/(2M_{n}))$ represents different
masses, and $V({\bf r})=\text{diag}(V_{1}({\bf r}),\ldots,V_{n}({\bf r}))$ so
we might have different external potentials for different components. The
matrix

$\Gamma(\vec{\psi})=\left[\begin{array}[]{cccc}\gamma_{11}\left|\psi_{1}\right|^{2}&\gamma_{12}\psi_{1}\bar{\psi}_{2}&\dots&\gamma_{1n}\psi_{1}\bar{\psi}_{n}\\\
\gamma_{12}\psi_{2}\bar{\psi}_{1}&\gamma_{22}\left|\psi_{2}\right|^{2}&\dots&\gamma_{2n}\psi_{2}\bar{\psi}_{n}\\\
\vdots&\vdots&\ddots&\vdots\\\
\gamma_{1n}\psi_{n}\bar{\psi}_{1}&\gamma_{2n}\psi_{n}\bar{\psi}_{2}&\dots&\gamma_{nn}\left|\psi_{n}\right|^{2}\end{array}\right],\gamma_{ij}\in\mathbb{R},$
(4)

specifies the intra- and inter-component interactions, which is repulsive if
$\gamma_{ij}>0$. Further, we assume that there are additional global
constraints $\widetilde{G}_{j}:\,\mathcal{H}\longrightarrow\mathbb{C}$ on the
form

$\widetilde{G}_{j}(\vec{\psi})=c_{j}-\int\vec{{\psi}}^{\dagger}K_{j}\vec{{\psi}}d{\bf
r}=0,\ j=1,\ldots,n_{c},\ c_{j}\in\mathbb{C},$ (5)

where $K_{j}$ is a matrix, to be specified later, not depending on
$\vec{\psi}$. Further, we define the constraint set

$\widetilde{\mathcal{G}}=\left\\{\vec{\psi}:\widetilde{G}_{j}(\vec{\psi})=0,j=1,\ldots,n_{c}\right\\}.$

Denoting the Lagrange multipliers as
$\lambda_{j}\in\mathbb{C},j=1,\ldots,n_{c}$, the Lagrange functional (extended
energy) can be written as

$\widetilde{I}(\vec{\psi},\vec{{\lambda}})=\widetilde{E}(\vec{\psi})+\sum_{j}\lambda_{j}\widetilde{G}_{j}(\vec{\psi}).$
(6)

It is well known [16] that if $\vec{\psi}^{*}$ is a local minimizer of

$\min_{\vec{{\psi}}\in\widetilde{\mathcal{G}}}\widetilde{E}\left(\vec{{\psi}}\right),$
(7)

then the first order necessary condition

$\dfrac{\delta\widetilde{I}}{\delta\vec{{\psi}}^{\dagger}}=\dfrac{\delta\widetilde{E}}{\delta\vec{{\psi}}^{\dagger}}+\sum_{j}\lambda_{j}\dfrac{\delta\widetilde{G}_{j}}{\delta\vec{{\psi}}^{\dagger}}=0$
(8)

is satisfied at $\vec{\psi}^{*}$.

If we define the projection of $\vec{\xi}\in\mathcal{H}$ onto the tangent
space of the constraint set $\widetilde{\mathcal{G}}$ as

$P_{\widetilde{\mathcal{G}}}^{\perp}(\vec{\xi})=\vec{\xi}-\sum_{j=1}^{n_{c}}\alpha_{j}(\vec{\xi})\dfrac{\delta\widetilde{G}_{j}}{\delta\vec{{\psi}}^{\dagger}},$
(9)

where $\alpha_{j}$ are given by the solution of the linear system

$T\vec{\alpha}=\vec{d},$ (10)

where

$t_{ij}=\left\langle\dfrac{\delta\widetilde{G}_{i}}{\delta\vec{{\psi}}^{\dagger}}\,\middle|\dfrac{\delta\widetilde{G}_{j}}{\delta\vec{{\psi}}^{\dagger}}\right\rangle,\,d_{i}=\left\langle\dfrac{\delta\widetilde{G}_{i}}{\delta\vec{{\psi}}^{\dagger}}\,\,\middle|\vec{\xi}\right\rangle,$
(11)

then

$P_{\widetilde{\mathcal{G}}}^{\perp}\left(\sum_{j}\lambda_{j}\dfrac{\delta\widetilde{G}_{j}}{\delta\vec{{\psi}}^{\dagger}}\right)=0,$

and (8) gives the a related necessary condition for a minimum

$P_{\widetilde{\mathcal{G}}}^{\perp}\left(\dfrac{\delta\widetilde{E}}{\delta\vec{{\psi}}^{\dagger}}\right)=0.$
(12)

Note that in (12) the Lagrange parameters are not present. This fact will be
used later.

To be more precise we get from (3) and (5) by using partial integration for
the first term that

$\dfrac{\delta\widetilde{E}(\vec{{\psi}})}{\delta\vec{{\psi}}^{\dagger}}=-W\mathop{}\\!\mathbin{\bigtriangleup}\vec{\psi}+V\vec{\psi}+2\pi\Gamma\left(\vec{{\psi}}\right)\vec{{\psi}},\
\frac{\delta\widetilde{G}_{j}(\vec{{\psi}})}{\delta\vec{\psi}^{\dagger}}=-K_{j}\vec{\psi},$

and we see that (8) is equivalent to the partial differential equation

$-W\mathop{}\\!\mathbin{\bigtriangleup}\vec{\psi}+V\vec{\psi}+2\pi\Gamma\left(\vec{{\psi}}\right)\vec{{\psi}}-\sum_{j}\lambda_{j}K_{j}\vec{{\psi}}=0,$
(13)

where we define the Laplace operator acting on a vector to be

$\mathop{}\\!\mathbin{\bigtriangleup}\vec{\psi}=\left[\begin{array}[]{c}\mathop{}\\!\mathbin{\bigtriangleup}\psi_{1}\\\
\vdots\\\ \mathop{}\\!\mathbin{\bigtriangleup}\psi_{n}\end{array}\right].$

##  3 General formulation of DFPM with constraints

The first step in the Dynamical Functional Particle Method [17], DFPM, is to
introduce an artificial time $\tau$ making all the functions and functionals
in the preceeding section depend on $\tau$. Thus, $\vec{\psi}({\bf r})$
become, say, $\vec{{\Psi}}({\bf r},\tau)$ and we drop all the tildes in the
notation used above to indicate the new dependence on $\tau$, e.g.,
$\widetilde{E}(\vec{\psi})$ will become $E(\vec{\Psi})$ . Then we formulate
the damped dynamical system

$\dfrac{\partial^{2}\vec{\Psi}}{\partial\tau^{2}}+\eta\dfrac{\partial\vec{\Psi}}{\partial\tau}=-\dfrac{\delta
I(\vec{{\Psi}},\vec{{\lambda}})}{\delta\vec{{\Psi}}^{\dagger}},$ (14)

with initial conditions $\vec{\Psi}_{0}=\vec{{\Psi}}({\bf
r},0),\vec{\Phi}_{0}=\frac{\partial\vec{\Psi}}{\partial\tau}({\bf r},0)$, and
where $\eta>0$ is a damping parameter.

We need to satisfy the constraints at the stationary point. This can be done
either by projection onto the constraint set or by formulating additional
damped dynamical systems for the constraints, see below. Assume that the
number of constraints satisfied by projection is $n_{P}$ then
$G_{j}(\vec{{\Psi}}(\tau))=0,j=1,\ldots,n_{P}\leq n_{c}$ and we have the
constraint set $\mathcal{D}(\tau)=$
$\left\\{\vec{\Psi}(\tau):G_{j}(\vec{\Psi}(\tau))=0,j=1,\ldots,n_{P}\right\\}$.
The functional derivatives in (8) should be taken on the constraint set
$\mathcal{D}(\tau)$, i.e.,

$\dfrac{\delta
I(\vec{{\Psi}},\vec{{\lambda}})}{\delta\vec{{\Psi}}^{\dagger}}\left.\right|_{\mathcal{D}}=P_{\mathcal{D}}^{\perp}\left(\dfrac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}+\sum_{j}\lambda_{j}\dfrac{\delta
G_{j}}{\delta\vec{{\Psi}}^{\dagger}}\right),$ (15)

where $P_{\mathcal{D}}^{\perp}$ is the (orthogonal) projection on the tangent
space of the constraint set $\mathcal{D}(\tau)$ defined as in (9). The
dynamical system corresponding to (14) becomes

$\dfrac{\partial^{2}\vec{\Psi}}{\partial\tau^{2}}+\eta\dfrac{\partial\vec{\Psi}}{\partial\tau}=-P_{\mathcal{D}}^{\perp}\left(\dfrac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}+\sum_{j}\lambda_{j}\dfrac{\delta
G_{j}}{\delta\vec{{\Psi}}^{\dagger}}\right),$ (16)

or since $P_{\mathcal{D}}^{\perp}\left(\frac{\delta
G_{j}}{\delta\vec{{\Psi}}^{\dagger}}\right)=0,j=1,\ldots,n_{P}$

$\dfrac{\partial^{2}\vec{\Psi}}{\partial\tau^{2}}+\eta\dfrac{\partial\vec{\Psi}}{\partial\tau}=-P_{\mathcal{D}}^{\perp}\left(\dfrac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}\right)-\sum_{j=n_{P}+1}^{n_{c}}\lambda_{j}P_{\mathcal{D}}^{\perp}\left(\dfrac{\delta
G_{j}}{\delta\vec{{\Psi}}^{\dagger}}\right).$ (17)

The remaining constraints are forced to be satisfied at the stationary
solution by solving additional damped dynamical systems

$\begin{array}[]{c}\dfrac{\partial^{2}G_{n_{P}+1}}{\partial\tau^{2}}+\eta\dfrac{\partial
G_{n_{P}+1}}{\partial\tau}=-k_{n_{P}+1}G_{n_{P}+1},\\\\[11.38109pt]
\dfrac{\partial^{2}G_{n_{P}+2}}{\partial\tau^{2}}+\eta\dfrac{\partial
G_{n_{P}+2}}{\partial\tau}=-k_{n_{P}+2}G_{n_{P}+2},\\\ \vdots\\\
\dfrac{\partial^{2}G_{n_{c}}}{\partial\tau^{2}}+\eta\dfrac{\partial
G_{n_{c}}}{\partial\tau}=-k_{n_{c}}G_{n_{c}},\end{array}$ (18)

where $k_{j}>0$ can be chosen in order to optimize the numerical performance
of the method. Note that by definition these equations ensure that
$G_{j}\rightarrow 0,j=n_{P}+1,\ldots,n_{c}$ when $\tau\rightarrow\infty$.

The unknown Lagrange parameters $\lambda_{j},j=n_{P}+1,\ldots,n_{c}$ are
determined by substituting the known $G_{j},j=n_{P}+1,\ldots,n_{c}$ into (18)
and then using (17) as we will show in detail in Sec. 3.1.

Finally, note that it is possible to have no projection, $n_{P}=0$, as well as
no additional damped systems for the constraints, $n_{P}=n_{c}$.

###  3.1 Finding the Lagrange parameters

First assume that the first $n_{P}$ constraints are satisfied. Further, for a
simpler notation, we use inner products (bra-ket) for the integrals. We will
use (18) to get the equations for the Lagrange parameters
$\lambda_{j},j=n_{P}+1,\ldots,n_{c}$ in (17), i.e., we need the first and
second order derivatives of $G_{j}$ w.r.t. $\tau$. For notational convenience
we drop the brackets for the projection, i.e.,
$P_{\mathcal{D}}^{\perp}(\vec{\xi})=P_{\mathcal{D}}^{\perp}\vec{\xi}$.

From (5) we have that

$\dfrac{\partial
G_{j}}{\partial\tau}=-\left\langle\dfrac{\partial\vec{\Psi}}{\partial\tau}\middle|K_{j}\vec{{\Psi}}\right\rangle-\left\langle\vec{\Psi}\middle|K_{j}\dfrac{\partial\vec{\Psi}}{\partial\tau}\right\rangle,$
(19)

and

$\dfrac{\partial^{2}G_{j}}{\partial\tau^{2}}=-\left\langle\dfrac{\partial^{2}\vec{\Psi}}{\partial\tau^{2}}\middle|K_{j}\vec{{\Psi}}\right\rangle-2\left\langle\dfrac{\partial\vec{\Psi}}{\partial\tau}\middle|K_{j}\dfrac{\partial\vec{\Psi}}{\partial\tau}\right\rangle-\left\langle\vec{{\Psi}}\middle|K_{j}\dfrac{\partial^{2}\vec{\Psi}}{\partial\tau^{2}}\right\rangle$
(20)

and with (19) and (20) we get

$\begin{array}[]{l}\dfrac{\partial^{2}G_{j}}{\partial\tau^{2}}+\eta\dfrac{\partial
G_{j}}{\partial\tau}=-\left\langle\dfrac{\partial^{2}\vec{\Psi}}{\partial\tau^{2}}\middle|K_{j}\vec{{\Psi}}\right\rangle-2\left\langle\dfrac{\partial\vec{\Psi}}{\partial\tau}\middle|K_{j}\dfrac{\partial\vec{\Psi}}{\partial\tau}\right\rangle-\left\langle\vec{{\Psi}}\middle|K_{j}\dfrac{\partial^{2}\vec{\Psi}}{\partial\tau^{2}}\right\rangle\\\
-\eta\left\langle\dfrac{\partial\vec{\Psi}}{\partial\tau}\middle|K_{j}\vec{{\Psi}}\right\rangle-\eta\left\langle\vec{\Psi}\middle|K_{j}\dfrac{\partial\vec{\Psi}}{\partial\tau}\right\rangle=\\\
-\left\langle\dfrac{\partial^{2}\vec{\Psi}}{\partial\tau^{2}}+\eta\dfrac{\partial\vec{\Psi}}{\partial\tau}\middle|K_{j}\vec{\Psi}\right\rangle-2\left\langle\dfrac{\partial\vec{\Psi}}{\partial\tau}\middle|K_{j}\dfrac{\partial\vec{\Psi}}{\partial\tau}\right\rangle-\left\langle\vec{\Psi}\middle|K_{j}\left(\dfrac{\partial^{2}\vec{\Psi}}{\partial\tau^{2}}+\eta\dfrac{\partial\vec{\Psi}}{\partial\tau}\right)\right\rangle.\end{array}$
(21)

With (17) inserted into (21) we get

$\begin{array}[]{l}\dfrac{\partial^{2}G_{j}}{\partial\tau^{2}}+\eta\dfrac{\partial
G_{j}}{\partial\tau}=-\left\langle P_{\mathcal{D}}^{\perp}\dfrac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}-\sum_{k}\lambda_{k}P_{\mathcal{D}}^{\perp}K_{k}\vec{{\Psi}}\middle|K_{j}\vec{{\Psi}}\right\rangle-2\left\langle\dfrac{\partial\vec{\Psi}}{\partial\tau}\middle|K_{j}\dfrac{\partial\vec{\Psi}}{\partial\tau}\right\rangle\\\
-\left\langle\vec{{\Psi}}\middle|K_{j}(P_{\mathcal{D}}^{\perp}\dfrac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}-\sum_{k}\lambda_{k}P_{\mathcal{D}}^{\perp}K_{k}\vec{{\Psi}})\right\rangle.\end{array}$
(22)

We want to solve for the Lagrange parameters so we expand the bra-ket’s in
(22)

$\begin{array}[]{l}\dfrac{\partial^{2}G_{j}}{\partial\tau^{2}}+\eta\dfrac{\partial
G_{j}}{\partial\tau}=-\left\langle P_{\mathcal{D}}^{\perp}\dfrac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}|K_{j}\vec{{\Psi}}\right\rangle+\sum_{k}\lambda_{k}\left\langle
P_{\mathcal{D}}^{\perp}K_{k}\vec{{\Psi}}\middle|K_{j}\vec{{\Psi}}\right\rangle\\\
-2\left\langle\dfrac{\partial\vec{\Psi}}{\partial\tau}\middle|K_{j}\dfrac{\partial\vec{\Psi}}{\partial\tau}\right\rangle-\left\langle\vec{{\Psi}}\middle|K_{j}P_{\mathcal{D}}^{\perp}\dfrac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}\right\rangle+\sum_{k}\lambda_{k}\left\langle\vec{{\Psi}}\middle|K_{j}P_{\mathcal{D}}^{\perp}K_{k}\vec{{\Psi}})\right\rangle,\end{array}$
(23)

and compare with the right hand side of (18). By defining the matrix
$A\in\mathbb{C}^{(n_{c}-n_{P})\times(n_{c}-n_{P})}$ and the vector
$\vec{b}\in\mathbb{C}^{n_{c}-n_{P}}$ with elements

$\displaystyle a_{j-n_{P},k-n_{P}}=\left\langle
P_{\mathcal{D}}^{\perp}K_{k}\vec{{\Psi}}\middle|K_{j}\vec{{\Psi}}\right\rangle+\left\langle\vec{{\Psi}}\middle|K_{j}P_{\mathcal{D}}^{\perp}K_{k}\vec{{\Psi}})\right\rangle,$
(24) $\displaystyle
b_{j-n_{P}}=k_{j}\left(c_{j}-\left\langle\vec{{\Psi}}\middle|K_{j}\vec{{\Psi}})\right\rangle\right)+$
$\displaystyle+\left\langle P_{\mathcal{D}}^{\perp}\dfrac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}\middle|K_{j}\vec{{\Psi}}\right\rangle+2\left\langle\dfrac{\partial\vec{\Psi}}{\partial\tau}\middle|K_{j}\dfrac{\partial\vec{\Psi}}{\partial\tau}\right\rangle+\left\langle\vec{{\Psi}}\middle|K_{j}P_{\mathcal{D}}^{\perp}\dfrac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}\right\rangle,$ (25)

where $j,k=n_{P}+1,\ldots,n_{c}$ the Lagrange parameters
$\vec{\lambda}\in\mathbb{C}^{n_{c}-n_{P}}$ are found by solving the linear
system

$A\vec{\lambda}=\vec{b}.$ (26)

The results in (24) and (25) are the cornerstone of our approach, valid for
arbitrary number of components and dimensions. In the sections to come we
apply it to specific problems and evaluate the resulting numerical
performance. However, we have to restrict ourselves to a limited number of
test cases in one and two dimensions.

##  4 Normalization and rotational constraints for $n$ components in one
dimension

Here we consider one spatial dimension with a constraint on the angular
momentum and on the interval $[-\pi,\pi]$ with periodic boundary conditions
$\vec{\Psi}(-\pi)=\vec{\Psi}(\pi)$, i.e., a ring geometry. From (3) we have
the energy

$E\left(\vec{{\Psi}}\right)=\int_{-\pi}^{\pi}\dfrac{\partial\vec{\Psi}^{\dagger}}{\partial
x}W\dfrac{\partial\vec{\Psi}}{\partial
x}+\vec{\Psi}^{\dagger}V\vec{\Psi}+\pi\vec{\Psi}^{\dagger}\Gamma\left(\vec{{\Psi}}\right)\vec{{\Psi}}\,{\rm
d}x,$ (27)

and

$\frac{\delta
E(\vec{{\Psi}})}{\delta\vec{{\Psi}}^{\dagger}}=-W\frac{\partial^{2}\vec{{\Psi}}}{\partial
x^{2}}+V\vec{\Psi}+2\pi\Gamma\left(\vec{{\Psi}}\right)\vec{{\Psi}}.$ (28)

The constraints are

$G_{j}=c_{j}-\int_{-\pi}^{\pi}\vec{{\Psi}}^{\dagger}K_{j}\vec{{\Psi}}{\rm
d}x=0,j=1,\ldots,n+1,$ (29)

where

$\begin{array}[]{l}\>K_{1}=\left[\begin{array}[]{cccc}1&0&\dots&0\\\
0&0&&\vdots\\\ \vdots&&\ddots&\vdots\\\
0&\dots&\dots&0\end{array}\right],\ldots,\>K_{n}=\left[\begin{array}[]{cccc}0&0&\dots&0\\\
0&0&&\vdots\\\ \vdots&&\ddots&\vdots\\\ 0&\dots&\dots&1\end{array}\right],\\\
\vspace{3mm}K_{n+1}=-i\left[\begin{array}[]{cccc}\frac{\partial}{\partial
x}&0&\dots&0\\\ 0&\frac{\partial}{\partial x}&&\vdots\\\ \vdots&&\ddots&0\\\
0&\dots&0&\frac{\partial}{\partial x}\end{array}\right].\end{array}$ (30)

In our general context we have here $n_{c}=n+1$ and the $K_{j}$ in addition
fulfills $K_{n+1}=L\sum_{k<n+1}K_{k}$ for the (angular-) momentum operator

$L=-i\frac{\partial}{\partial x}.$ (31)

For these constraints we have

$\dfrac{\delta
G_{j}}{\delta\vec{{\Psi}}^{\dagger}}=-K_{j}\vec{\Psi}_{j},\,j=1,\ldots,n+1.$

There are some special cases that are particularly interesting.

We may project all the constraints giving no Lagrange parameters to determine.
However, this requires the constraints to be satisfied, which for nonlinear
constraints usually is a non-trivial computational task. For details, see the
following subsection.

There are two other special cases that for our problem setting are more
interesting from an algorithmic point of view and lends themselves to an
efficient algorithm and numerical tests.

Firstly, it is natural to consider the case of projecting only some or all of
the normalization constraints, $1\leq n_{P}\leq n$, since the projection in
practice is just a trivial normalization.

Secondly, we may consider only dynamically damped constraints, $n_{P}=0$, that
will not force any of the constraints to be satisfied except at the stationary
point, which we have also experienced can lead to less sensitivity on the
initial conditions.

###  4.1 Only projection no dynamically damped constraints

Let us first look at the case $n_{P}=n+1$, i.e., we project all the
constraints. Here, (17) becomes

$\dfrac{\partial^{2}\vec{\Psi}}{\partial\tau^{2}}+\eta\dfrac{\partial\vec{\Psi}}{\partial\tau}=-P_{\mathcal{D}}^{\perp}\left(\dfrac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}\right).$ (32)

In this case we do not have any unknown Lagrange parameters and we need only
to determine the form of the projection $P_{\mathcal{D}}^{\perp}$ in (9). This
we get from (11), i.e.,

$t_{ij}=\left\langle
K_{i}\vec{\Psi}\middle|K_{j}\vec{\Psi}\right\rangle,\,d_{i}=-\left\langle
K_{i}\vec{\Psi}\,\,\middle|\vec{\xi}\right\rangle,$

or in more detail

$\displaystyle
T=\left[\begin{array}[]{cccc}c_{1}&&&\left\langle{\Psi}_{1}\middle|L{\Psi}_{n}\right\rangle\\\
&\ddots&&\vdots\\\
&&c_{n}&\left\langle{\Psi}_{n}\middle|L{\Psi}_{n}\right\rangle\\\\[5.69054pt]
\left\langle{\Psi}_{1}\middle|L{\Psi}_{n}\right\rangle&\cdots&\left\langle{\Psi}_{n}\middle|L{\Psi}_{n}\right\rangle&\sum_{l=1}^{n}\left\langle\Psi_{l}\middle|L^{2}\Psi_{l}\right\rangle\end{array}\right],$
(37)
$\displaystyle\vec{d}=-\left[\begin{array}[]{c}\left\langle\Psi_{1}\middle|\xi_{1}\right\rangle\\\
\vdots\\\ \left\langle\Psi_{n}\middle|\xi_{n}\right\rangle\\\\[5.69054pt]
\sum_{l=1}^{n}\left\langle
L\Psi_{l}\middle|\xi_{l}\right\rangle\end{array}\right].$ (42)

The linear system given by the matrix $T$ and vector $\vec{d}$ in (37) can be
efficiently solved, with order $n$ operations, using sparse Gaussian
elimination.

However, the constraints have to be satisfied. This is trivial in the case of
normalization constraints since any approximate numerical solution can be
normalized without almost any cost. The constraint on the angular momentum is
more complicated to satisfy since it’s nonlinear. However, we believe the idea
of using only projection is principally interesting and thus we describe a
possible formulation. Assume that we have $\vec{\Psi}(x,\tau)$ from the
solution of our dynamical system (32) that is not generally satisfying the
constraints. In order to do so we can solve the least norm minimization
problem

$\min_{\vec{\xi}\in\mathcal{D}(\tau)}\|\vec{\xi}-\vec{\Psi}\|,$ (43)

or in more detail

$\displaystyle\min_{\vec{\xi}}\|\vec{\xi}-\vec{\Psi}\|$ (44)
$\displaystyle\text{subject to }\
c_{i}-\left\langle{\xi}_{i}\middle|{\xi}_{i}\right\rangle=0,i=1,\ldots,n,$
(45) $\displaystyle\hskip
42.67912ptc_{n+1}-\sum_{i=1}^{n}\left\langle{\xi}_{i}\middle|L{\xi}_{i}\right\rangle=0.$
(46)

Finding the optimal $\vec{\xi}(x)$ above is certainly possible but requires an
iterative process that will presumably be much more costly in a numerical
implementation than the other two approaches described below. Therefore, we do
not expect this to be an efficient algorithm and will not include this
approach in the numerical tests.

###  4.2 Projection of one or more normalization constraints

We start by consider the projection of the first $n_{P}$ normalization
constraints. Then the projection is explicitly given as

$P_{\mathcal{D}}^{\perp}(\vec{\xi})=\vec{\xi}-\left[\Psi_{1}\langle\Psi_{1}\middle|\xi_{1}\rangle/c_{1},\Psi_{2}\langle\Psi_{2}\middle|\xi_{2}\rangle/c_{2},\ldots,\Psi_{n_{P}}\langle\Psi_{n_{P}}\middle|\xi_{n_{P}}\rangle/c_{n_{P}},0,\ldots,0\right]^{T}$

and

$\dfrac{\partial^{2}\vec{\Psi}}{\partial\tau^{2}}+\eta\dfrac{\partial\vec{\Psi}}{\partial\tau}=-P_{\mathcal{D}}^{\perp}\left(\dfrac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}\right)-\sum_{j=n_{P}+1}^{n+1}\lambda_{j}P_{\mathcal{D}}^{\perp}\left(K_{j}\vec{\Psi}\right).$
(47)

We can simplify (47) slightly by noting that
$P_{\mathcal{D}}^{\perp}\left(K_{j}\vec{\Psi}\right)=\Psi_{j}\vec{e}_{j},n_{P}+1\leq
j\leq n$ giving

$\dfrac{\partial^{2}\vec{\Psi}}{\partial\tau^{2}}+\eta\dfrac{\partial\vec{\Psi}}{\partial\tau}=-P_{\mathcal{D}}^{\perp}\left(\dfrac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}\right)-\sum_{j=n_{P}+1}^{n}\lambda_{j}\Psi_{j}\vec{e}_{j}-\lambda_{n+1}P_{\mathcal{D}}^{\perp}\left(K_{n+1}\vec{\Psi}\right).$
(48)

In the following lemma we give the details necessary to efficiently calculate
the Lagrange parameters. For notational convenience we drop the brackets for
the projection, i.e.,
$P_{\mathcal{D}}^{\perp}(\vec{\xi})=P_{\mathcal{D}}^{\perp}\vec{\xi}$.

######  Lemma 1.

In (24) and (25) we have

$\left\langle
P_{\mathcal{D}}^{\perp}K_{k}\vec{{\Psi}}\middle|K_{j}\vec{{\Psi}}\right\rangle=\left\langle\vec{{\Psi}}\middle|K_{j}P_{\mathcal{D}}^{\perp}K_{k}\vec{{\Psi}}\right\rangle=$

$\displaystyle\delta_{kj}c_{j},\ n_{P}+1\leq j,k\leq n,$ (49)
$\displaystyle\left\langle{\Psi}_{j}\middle|L{\Psi}_{n}\right\rangle,\ k=n+1,\
j=n_{P}+1,\ldots n$ (50)
$\displaystyle\sum_{l=1}^{n}\left\langle\Psi_{l}\middle|L^{2}\Psi_{l}\right\rangle-\sum_{l=1}^{m}\left\langle\Psi_{l}\middle|L\Psi_{l}\right\rangle^{2}/c_{l},\
j=k=n+1,$ (51)

and

$\left\langle P_{\mathcal{D}}^{\perp}\frac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}\middle|K_{j}\vec{{\Psi}}\right\rangle=\left\langle\vec{{\Psi}}\middle|K_{j}P_{\mathcal{D}}^{\perp}\frac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}\right\rangle=$ (52)
$\displaystyle\left\langle\Psi_{j}\middle|\frac{\delta
E}{\delta\overline{\Psi}_{j}}\right\rangle,j=n_{P}+1,\ldots,n,$ (53)
$\displaystyle\sum_{l=1}^{n}\left\langle\frac{\delta
E}{\delta\overline{\Psi}_{l}}\middle|L\Psi_{l}\right\rangle-\sum_{l=1}^{n_{P}}\left\langle\Psi_{l}\middle|L\Psi_{l}\right\rangle\left\langle\frac{\delta
E}{\delta\overline{\Psi}_{l}}\middle|\Psi_{l}\right\rangle/c_{l},j=n+1.$ (54)

###### Proof.

First consider (49) and note that $K_{k}\vec{{\Psi}}=\Psi_{k}\vec{e}_{k}$ and
thus $P_{\mathcal{D}}^{\perp}K_{k}\vec{{\Psi}}=\Psi_{k}\vec{e}_{k}$ giving
$\left\langle
P_{\mathcal{D}}^{\perp}K_{k}\vec{{\Psi}}\middle|K_{j}\vec{{\Psi}}\right\rangle=\left\langle\Psi_{k}\vec{e}_{k}\middle|\Psi_{j}\vec{e}_{j}\right\rangle=\delta_{kj}c_{j}$.
A similar argument proves that
$\left\langle\vec{{\Psi}}\middle|K_{j}P_{\mathcal{D}}^{\perp}K_{k}\vec{{\Psi}})\right\rangle=\delta_{kj}c_{j}$.

Turning to (50) we first note that $\left\langle
P_{\mathcal{D}}^{\perp}K_{n+1}\vec{{\Psi}}\middle|K_{j}\vec{{\Psi}}\right\rangle=\left\langle
P_{\mathcal{D}}^{\perp}K_{n+1}\vec{{\Psi}}\middle|\Psi_{j}\vec{e}_{j}\right\rangle=\left\langle
L{\Psi}_{j}\middle|\Psi_{j}\right\rangle$ and similarly
$\left\langle\vec{{\Psi}}\middle|K_{j}P_{\mathcal{D}}^{\perp}K_{k}\vec{{\Psi}}\right\rangle=\left\langle\vec{{\Psi}}\middle|K_{j}P_{\mathcal{D}}^{\perp}K_{n+1}\vec{{\Psi}}\right\rangle=\left\langle\vec{{\Psi}}\middle|L\Psi_{j}\vec{e}_{j}\right\rangle=\left\langle\Psi_{j}\middle|L\Psi_{j}\right\rangle$.
By partial integration we get $\left\langle
L{\Psi}_{j}\middle|\Psi_{j}\right\rangle=\left\langle\Psi_{j}\middle|L\Psi_{j}\right\rangle$.

To prove (51) we look in detail on

$\displaystyle\left\langle
P_{\mathcal{D}}^{\perp}K_{n+1}\vec{{\Psi}}\middle|K_{n+1}\vec{{\Psi}}\right\rangle=$
(55)
$\displaystyle\sum_{l=1}^{n_{P}}\left\langle(L\Psi_{l}-\Psi_{l}\langle\Psi_{1}\middle|L\Psi_{l}\rangle/c_{l}\middle|L\Psi_{l}\right\rangle+\sum_{l=n_{P}+1}^{n}\left\langle
L\Psi_{l}\middle|L\Psi_{l}\right\rangle=$ (56)
$\displaystyle\sum_{l=1}^{n}\left\langle
L\Psi_{l}\middle|L\Psi_{l}\right\rangle-\sum_{l=1}^{n_{P}}\left\langle\Psi_{l}\middle|L\Psi_{l}\right\rangle^{2}/c_{l},$
(57)

and

$\displaystyle\left\langle\vec{{\Psi}}\middle|K_{n+1}P_{\mathcal{D}}^{\perp}K_{n+1}\vec{{\Psi}}\right\rangle=$
(58)
$\displaystyle\sum_{l=1}^{n_{P}}\left\langle\Psi_{l}\middle|L(L\Psi_{l}-\Psi_{l}\left\langle\Psi_{1}\middle|L\Psi_{l}\right\rangle/c_{l})\right\rangle+\sum_{l=n_{P}+1}^{n}\left\langle\Psi_{l}\middle|L^{2}\Psi_{l}\right\rangle=$
(59)
$\displaystyle\sum_{l=1}^{n}\left\langle\Psi_{l}\middle|L^{2}\Psi_{l}\right\rangle-\sum_{l=1}^{n_{P}}\left\langle\Psi_{l}\middle|L\Psi_{l}\right\rangle^{2}/c_{l}.$
(60)

Using partial integration we have $\left\langle
L\Psi_{l}\middle|L\Psi_{l}\right\rangle=\left\langle\Psi_{l}\middle|L^{2}\Psi_{l}\right\rangle$
showing the result.

In (53) we again notice that $K_{j}\vec{{\Psi}}=\Psi_{j}\vec{e}_{j}$ and
$P_{\mathcal{D}}^{{\perp}^{\dagger}}K_{j}\vec{{\Psi}}=\Psi_{j}\vec{e}_{j}$ so
it remains to prove that $\left\langle\Psi_{j}\middle|\frac{\delta
E}{\delta\overline{\Psi}_{j}}\right\rangle=\left\langle\frac{\delta
E}{\delta\overline{\Psi}_{j}}\middle|\Psi_{j}\right\rangle$ which can be done
by integration by parts.

Finally we prove (54) by first looking at

$\displaystyle\left\langle P_{\mathcal{D}}^{\perp}\frac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}\middle|K_{n+1}\vec{{\Psi}}\right\rangle=\sum_{l=1}^{n_{P}}\left\langle\frac{\delta
E}{\delta\overline{\Psi}_{l}}-\Psi_{l}\left\langle\Psi_{l}\middle|\frac{\delta
E}{\delta\overline{\Psi}_{l}}\right\rangle/c_{l}\middle|L\Psi_{l}\right\rangle+\sum_{l=n_{P}+1}^{n}\left\langle\frac{\delta
E}{\delta\overline{\Psi}_{l}}\middle|L\Psi_{l}\right\rangle=$ (61)
$\displaystyle\sum_{l=1}^{n}\left\langle\frac{\delta
E}{\delta\overline{\Psi}_{l}}\middle|L\Psi_{l}\right\rangle-\sum_{l=1}^{n_{P}}\left\langle\Psi_{l}\middle|L\Psi_{l}\right\rangle\left\langle\Psi_{l}\middle|\frac{\delta
E}{\delta\overline{\Psi}_{l}}/c_{l}\right\rangle,$ (62)

and then

$\displaystyle\left\langle\vec{{\Psi}}\middle|K_{n+1}P_{\mathcal{D}}^{\perp}\frac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}\right\rangle=\sum_{l=1}^{m}\left\langle\Psi_{l}\middle|L\frac{\delta
E}{\delta\overline{\Psi}_{l}}-L\Psi_{l}\left\langle\Psi_{l}\middle|\frac{\delta
E}{\delta\overline{\Psi}_{l}}\right\rangle\right\rangle+\sum_{l=n_{P}+1}^{n}\left\langle\Psi_{l}\middle|L\frac{\delta
E}{\delta\overline{\Psi}_{l}}\right\rangle=$ (63)
$\displaystyle\sum_{l=1}^{n}\left\langle\Psi_{l}\middle|L\frac{\delta
E}{\delta\overline{\Psi}_{l}}\right\rangle-\sum_{l=1}^{n_{P}}\left\langle\Psi_{l}\middle|\frac{\delta
E}{\delta\overline{\Psi}_{l}}\right\rangle\left\langle\Psi_{l}\middle|L\Psi_{l}\right\rangle/c_{l}.$
(64)

By partial integration we get (54). ∎

Using the lemma above it is easy to see that in (24) and (25) we have

$\displaystyle
A=2\left[\begin{array}[]{cccc}\left\langle{\Psi}_{n_{P}+1}\middle|{\Psi}_{n_{P}+1}\right\rangle&&&\left\langle{\Psi}_{n_{P}+1}\middle|L{\Psi}_{n}\right\rangle\\\
&\ddots&&\vdots\\\
&&\left\langle{\Psi}_{n}\middle|{\Psi}_{n}\right\rangle&\left\langle{\Psi}_{n}\middle|L{\Psi}_{n}\right\rangle\\\
\left\langle{\Psi}_{n_{P}+1}\middle|L{\Psi}_{n}\right\rangle&\cdots&\left\langle{\Psi}_{n}\middle|L{\Psi}_{n}\right\rangle&\sum_{l=1}^{n}\left\langle\Psi_{l}\middle|L^{2}\Psi_{l}\right\rangle-\sum_{l=1}^{n_{P}}\left\langle\Psi_{l}\middle|L\Psi_{l}\right\rangle^{2}/c_{l}\end{array}\right],$
(69)
$\displaystyle\vec{b}=2\left[\begin{array}[]{c}k_{n_{P}+1}G_{n_{P}+1}/2+\left\langle\Psi_{n_{P}+1}\middle|\dfrac{\delta
E}{\delta\overline{\Psi}_{1}}\right\rangle-\left\langle\dfrac{\partial\Psi_{1}}{\partial\tau}\middle|\dfrac{\partial\Psi_{1}}{\partial\tau}\right\rangle\\\
\vdots\\\ k_{n}G_{n}/2+\left\langle\Psi_{n}\middle|\dfrac{\delta
E}{\delta\overline{\Psi}_{n}}\right\rangle-\left\langle\dfrac{\partial\Psi_{n}}{\partial\tau}\middle|\dfrac{\partial\Psi_{n}}{\partial\tau}\right\rangle\\\
k_{n+1}G_{n+1}/2+\sum_{l=1}^{n}\left\langle\dfrac{\delta
E}{\delta\overline{\Psi}_{l}}\middle|L\Psi_{l}\right\rangle-\left\langle\dfrac{\partial\Psi_{l}}{\partial\tau}\middle|L\dfrac{\partial\Psi_{l}}{\partial\tau}\right\rangle-\sum_{l=1}^{n_{P}}\left\langle\Psi_{l}\middle|L\Psi_{l}\right\rangle\left\langle\dfrac{\delta
E}{\delta\overline{\Psi}_{l}}\middle|\Psi_{l}\right\rangle/c_{l}\end{array}\right].$
(74)

The system (26) with the matrix in (69) can be solved efficiently using sparse
Gaussian elimination with the number of operations of order $n-n_{p}$.

###  4.3 Projection of all normalization constraints

Let us consider $n_{P}=n$, i.e., we project all normalization constraints but
not the last constraint on the angular momentum. This will give only one
unknown Lagrange parameter and in practice a trivial projection of the
normalization constraints $\Psi_{j}=\sqrt{c_{j}}\Psi_{j}/\|\Psi_{j}\|$. The
projection is explicitly given as

$P_{\mathcal{D}}^{\perp}(\vec{\xi})=\vec{\xi}-\left[\Psi_{1}\langle\Psi_{1}\middle|\xi_{1}\rangle/c_{1},\Psi_{2}\langle\Psi_{2}\middle|\xi_{2}\rangle/c_{2},\ldots,\Psi_{n}\langle\Psi_{n}\middle|\xi_{n}\rangle/c_{n}\right]^{T},$
(75)

and

$\dfrac{\partial^{2}\vec{\Psi}}{\partial\tau^{2}}+\eta\dfrac{\partial\vec{\Psi}}{\partial\tau}=-P_{\mathcal{D}}^{\perp}\left(\dfrac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}\right)-\lambda_{n+1}P_{\mathcal{D}}^{\perp}\left(K_{n+1}\vec{\Psi}\right).$
(76)

The Lagrange parameter is

$\lambda_{n+1}=\dfrac{k_{n+1}G_{n+1}/2+\sum_{l=1}^{n}\left\langle\dfrac{\delta
E}{\delta\overline{\Psi}_{l}}\middle|L\Psi_{l}\right\rangle-\left\langle\Psi_{l}\middle|L\Psi_{l}\right\rangle\left\langle\dfrac{\delta
E}{\delta\overline{\Psi}_{l}}\middle|\Psi_{l}\right\rangle/c_{l}-\left\langle\dfrac{\partial\Psi_{l}}{\partial\tau}\middle|L\dfrac{\partial\Psi_{l}}{\partial\tau}\right\rangle}{\sum_{l=1}^{n}\left\langle\Psi_{l}\middle|L^{2}\Psi_{l}\right\rangle-\left\langle\Psi_{l}\middle|L\Psi_{l}\right\rangle^{2}/c_{l}}.$
(77)

###  4.4 Only dynamically damped constraints

No projection means that $P_{\mathcal{D}}^{\perp}=I$ in (17), $n_{P}=0$, and
the number of unknown Lagrange parameters is $n+1$ and thus

$\dfrac{\partial^{2}\vec{\Psi}}{\partial\tau^{2}}+\eta\dfrac{\partial\vec{\Psi}}{\partial\tau}=-\dfrac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}-\sum_{j=1}^{n+1}\lambda_{j}\dfrac{\delta
G_{j}}{\delta\vec{{\Psi}}^{\dagger}},$ (78)

where

$\dfrac{\delta
G_{j}}{\delta\vec{{\Psi}}^{\dagger}}=-\Psi_{j}\vec{e}_{j},j=1,\ldots,n,\dfrac{\delta
G_{n+1}}{\delta\vec{{\Psi}}^{\dagger}}=-\sum_{j=1}^{n}L\Psi_{j}\vec{e}_{j}.$

The Lagrange parameters are found by solving the $(n+1)\times(n+1)$ linear
system $A\vec{\lambda}=\vec{b}$ where we from (69) and (74) have

$\displaystyle
A=\left[\begin{array}[]{cccc}\left\langle{\Psi}_{1}\middle|{\Psi}_{1}\right\rangle&&&\left\langle{\Psi}_{1}\middle|{\Psi}_{n}\right\rangle\\\
&\ddots&&\vdots\\\
&&\left\langle{\Psi}_{n}\middle|{\Psi}_{n}\right\rangle&\left\langle{\Psi}_{n}\middle|L{\Psi}_{n}\right\rangle\\\
\left\langle{\Psi}_{1}\middle|L{\Psi}_{n}\right\rangle&\cdots&\left\langle{\Psi}_{n}\middle|L{\Psi}_{n}\right\rangle&\sum_{l=1}^{n}\left\langle\Psi_{l}\middle|L^{2}\Psi_{l}\right\rangle\end{array}\right],$
(83)
$\displaystyle\vec{b}=\left[\begin{array}[]{c}k_{1}G_{1}/2+\left\langle\Psi_{1}\middle|\dfrac{\delta
E}{\delta\overline{\Psi}_{1}}\right\rangle-\left\langle\dfrac{\partial\Psi_{1}}{\partial\tau}\middle|\dfrac{\partial\Psi_{1}}{\partial\tau}\right\rangle\\\
\vdots\\\ k_{n}G_{n}/2+\left\langle\Psi_{n}\middle|\dfrac{\delta
E}{\delta\overline{\Psi}_{n}}\right\rangle-\left\langle\dfrac{\partial\Psi_{n}}{\partial\tau}\middle|\dfrac{\partial\Psi_{n}}{\partial\tau}\right\rangle\\\\[8.53581pt]
k_{n+1}G_{n+1}/2+\sum_{l=1}^{n}\left\langle\dfrac{\delta
E}{\delta\overline{\Psi}_{l}}\middle|L\Psi_{l}\right\rangle-\left\langle\dfrac{\partial\Psi_{l}}{\partial\tau}\middle|L\dfrac{\partial\Psi_{l}}{\partial\tau}\right\rangle\end{array}\right].$
(88)

For $n=1$ we have

$\displaystyle
A=\left[\begin{array}[]{cc}\left\langle{\Psi}_{1}\middle|{\Psi}_{1}\right\rangle&\left\langle{\Psi}_{1}\middle|L{\Psi}_{1}\right\rangle\\\
\left\langle{\Psi}_{1}\middle|L{\Psi}_{1}\right\rangle&\left\langle\Psi_{1}\middle|L^{2}\Psi_{1}\right\rangle\end{array}\right],$
(91)
$\displaystyle\vec{b}=\left[\begin{array}[]{c}k_{1}G_{1}/2+\left\langle\Psi_{1}\middle|\dfrac{\delta
E}{\delta\overline{\Psi}_{1}}\right\rangle-\left\langle\dfrac{\partial\Psi_{1}}{\partial\tau}\middle|\dfrac{\partial\Psi_{1}}{\partial\tau}\right\rangle\\\
k_{2}G_{2}/2+\left\langle\dfrac{\delta
E}{\delta\overline{\Psi}_{1}}\middle|L\Psi_{1}\right\rangle-\left\langle\dfrac{\partial\Psi_{1}}{\partial\tau}\middle|L\dfrac{\partial\Psi_{1}}{\partial\tau}\right\rangle\end{array}\right],$
(94)

which agree with (B.10) in [18].

##  5 Numerical tests

In this section we present numerical calculations for some representative
cases in 1D and 2D.

We want to solve (17) through an iterative method in artificial time and start
by writing (17) as the first order system

$\displaystyle\frac{\partial\vec{\Psi}}{\partial\tau}=\vec{\Phi},$ (95)
$\displaystyle\frac{\partial\vec{\Phi}}{\partial\tau}=-\eta\vec{\Phi}+\vec{F}(\vec{\Psi}),$

where

$\vec{F}(\vec{\Psi})=-P_{\mathcal{D}}^{\perp}\left(\dfrac{\delta
E}{\delta\vec{{\Psi}}^{\dagger}}\right)-\sum_{j=n_{P}+1}^{n_{c}}\lambda_{j}P_{\mathcal{D}}^{\perp}\left(\dfrac{\delta
G_{j}}{\delta\vec{{\Psi}}^{\dagger}}\right).$ (96)

We assume that we have $N$ discretization points in space and store the
numerical values of $\vec{\Psi},\vec{\Phi},\vec{F}$ at iteration $k$ as
matrices $\mathbf{U},\mathbf{V},\mathbf{F}\in\mathbb{R}^{N\times n}$,
respectively. Then for some initial conditions $\mathbf{U}_{0},\mathbf{V}_{0}$
we can apply the symplectic Euler algorithm, see [19], to (95) to get the
following semi-implicit method

$\displaystyle\mathbf{V}_{k+1}=(I-\eta\,\Delta\tau)\mathbf{V}_{k}+\Delta\tau\,\mathbf{F}_{k},$
(97)
$\displaystyle\mathbf{U}_{k+1}=\mathbf{U}_{k}+\Delta\tau\,\mathbf{V}_{k+1}\>,\,k=0,1,...,$

where $\Delta\tau$ is the step in artificial time. Note that the symplectic
Euler above has a dual version where $\mathbf{U}_{k}$ is updated before
$\mathbf{V}_{k}$ but with the same numerical properties (convergence and
efficiency) [19].

Figure 1: (Color online) Stationary solution of (101), (102) with projection
for the normalization constants, see Sec. 5.1.1. Modulus (left) and phases
(mid), of $\Psi_{j}(x),\>j=1,2$. Compare with FIG. 4 (most right) of [20].
Convergence (right) for $\vec{\Psi}$ and the constraint $G_{3}$ of (99), see
inset legend. We have used the numerical parameters $\eta=k_{3}=2$ and
$\Delta\tau=0.99\Delta x\simeq 6.2\cdot 10^{-3}$. Physical parameters are set
to $\gamma=\gamma_{11}=\gamma_{12}=\gamma_{22}=23$, $c_{1}=0.96$,
$c_{2}=0.04$, and $c_{3}=0.45$, which is the same as in [20].

###  5.1 Multicomponent NLSE in 1D under rotation in a ring, applications to
vector solitons

In these examples we want to find normalized solutions $\vec{\Psi}(x)$ of (28)
where $W=\text{diag}(1,\ldots,1)$ and $V=\text{diag}(0,\ldots,0)$, i.e.,

$\frac{\delta
E(\vec{{\Psi}})}{\delta\vec{{\Psi}}^{\dagger}}=-\frac{\partial^{2}\vec{{\Psi}}}{\partial
x^{2}}+2\pi\Gamma\left(\vec{{\Psi}}\right)\vec{{\Psi}}=0,$ (98)

with a constrained angular momentum.

####  5.1.1 Projection of all the normalization constraints

From Eqs. (29), (30), and (31) we are left with one dynamic constraint for the
angular momentum

$G_{n+1}=c_{n+1}+i\sum_{j=1}^{n}\int\limits_{-\pi}^{\pi}\overline{\Psi}_{j}\frac{\partial}{\partial
x}\Psi_{j}\,{\rm d}x=0.$ (99)

We obtain the Lagrange parameter $\lambda_{n+1}$ in each time step for the
angular momentum from (77), with $L$ from (31) and $\delta
E(\vec{{\Psi}})/\delta\vec{{\Psi}}^{\dagger}$ from (98). The normalization of
each component, $\int\left|\Psi_{j}\right|^{2}dx=c_{j}$, is here fulfilled by
(trivial) projection, $\Psi_{j}=\sqrt{c_{j}}\Psi_{j}/\|\Psi_{j}\|$, in each
time step.

####  5.1.2 Only dynamically damped constraints

From Eqs. (29), (30), and (31) we end up with $n+1$ dynamic constraints

$G_{j}=c_{j}-\int\limits_{-\pi}^{\pi}\left|\Psi_{j}\right|^{2}\,{\rm d}x=0,\
j=1,\>...,n,\
G_{n+1}=c_{n+1}+i\sum_{j=1}^{n}\int\limits_{-\pi}^{\pi}\overline{\Psi}_{j}\frac{\partial}{\partial
x}\Psi_{j}\,{\rm d}x=0.$ (100)

We obtain the corresponding Lagrange parameters $\vec{\lambda}$ from the
linear system given by the matrix $A$ in (83) and the vector $\vec{b}$ in (88)
by solving $A\vec{\lambda}=\vec{b}$ in each time step.

####  5.1.3 Numerical benchmarking for $n=2$ components

In the special case of two components there is a rich literature on
applications in for example optics and cold atomic gases, see, e.g., [20, 21,
22] and references therein.

For the ease of the readers familiar with those applications, we here
explicitly write the corresponding dynamic equations (95) and (96) for the
stationary equation (98) in the case of two components. With
$\vec{\Psi}=\left[\Psi_{1},\Psi_{2}\right]^{T}$ and
$\vec{\Phi}=\left[\partial\Psi_{1}/\partial\tau,\partial\Psi_{2}/\partial\tau\right]^{T}$,
we have

$\frac{\partial^{2}}{\partial x^{2}}\left[\begin{array}[]{c}\Psi_{1}\\\
\Psi_{2}\end{array}\right]-2\pi\left[\begin{array}[]{cc}\gamma_{11}\left|\Psi_{1}\right|^{2}&\gamma_{12}\Psi_{1}\overline{\Psi}_{2}\\\
\gamma_{21}\overline{\Psi}_{1}\Psi_{2}&\gamma_{22}\left|\Psi_{2}\right|^{2}\end{array}\right]\left[\begin{array}[]{c}\Psi_{1}\\\
\Psi_{2}\end{array}\right]$ $+\lambda_{1}\left[\begin{array}[]{c}\Psi_{1}\\\
0\end{array}\right]+\lambda_{2}\left[\begin{array}[]{c}0\\\
\Psi_{2}\end{array}\right]+\lambda_{3}L\left[\begin{array}[]{c}\Psi_{1}\\\
\Psi_{2}\end{array}\right]=0,$ (101)

where $\lambda_{1}=\lambda_{2}=0$, and (77)

$\lambda_{3}=\dfrac{k_{3}G_{3}/2+\sum_{l=1}^{2}\left\langle\dfrac{\delta
E}{\delta\overline{\Psi}_{l}}\middle|L\Psi_{l}\right\rangle-\left\langle\Psi_{l}\middle|L\Psi_{l}\right\rangle\left\langle\dfrac{\delta
E}{\delta\overline{\Psi}_{l}}\middle|\Psi_{l}\right\rangle/c_{l}-\left\langle\Phi_{l}\middle|L\Phi_{l}\right\rangle}{\sum_{l=1}^{2}\left\langle\Psi_{l}\middle|L^{2}\Psi_{l}\right\rangle-\left\langle\Psi_{l}\middle|L\Psi_{l}\right\rangle^{2}/c_{l}},$
(102)

for the case of projecting the two normalization constraints. The case of only
dynamically damped constraints, the three Lagrange parameters are in each time
step given from the linear system (83), (88)

$\left[\begin{array}[]{ccc}\left\langle{\Psi}_{1}\middle|{\Psi}_{1}\right\rangle&0&\left\langle{\Psi}_{1}\middle|L{\Psi}_{1}\right\rangle\\\
0&\left\langle{\Psi}_{2}\middle|{\Psi}_{2}\right\rangle&\left\langle{\Psi}_{2}\middle|L{\Psi}_{2}\right\rangle\\\
\left\langle{\Psi}_{1}\middle|L{\Psi}_{1}\right\rangle&\left\langle{\Psi}_{2}\middle|L{\Psi}_{2}\right\rangle&\left\langle{\Psi}_{1}\middle|L^{2}{\Psi}_{1}\right\rangle+\left\langle{\Psi}_{2}\middle|L^{2}{\Psi}_{2}\right\rangle\end{array}\right]\left[\begin{array}[]{c}\lambda_{1}\\\
\lambda_{2}\\\ \lambda_{3}\end{array}\right]$
$=\left[\begin{array}[]{c}k_{1}G_{1}/2+\left\langle\Psi_{1}\middle|\dfrac{\delta
E}{\delta\overline{\Psi}_{1}}\right\rangle-\|\Phi_{1}\|^{2}\\\
k_{2}G_{2}/2+\left\langle\Psi_{2}\middle|\dfrac{\delta
E}{\delta\overline{\Psi}_{2}}\right\rangle-\|\Phi_{2}\|^{2}\\\
k_{3}G_{3}/2+\left\langle\Psi_{1}\middle|\dfrac{\delta
E}{\delta\overline{\Psi}_{1}}\right\rangle+\left\langle\Psi_{2}\middle|\dfrac{\delta
E}{\delta\overline{\Psi}_{2}}\right\rangle-\left\langle\Phi_{1}\middle|L\Phi_{1}\right\rangle+\left\langle\Phi_{2}\middle|L\Phi_{2}\right\rangle\end{array}\right].$
(103)

Above $L$ is the angular momentum operator (31), $G_{1}$, $G_{2}$, $G_{3}$ are
taken from (100), and the functional derivatives of the energy are given by
(98).

To compare our numerical results in 1D with the literature, we first
reproduced some well-known solutions with (101), which can be given
analytically for certain parameter values [20]. We perform the calculations
with both the dynamic methods presented here for obtaining the Lagrange
parameters. We present the numerical results for the case of projecting all
the normalization constraints, described in Sec. 5.1.1, in Fig. 1, and for the
case of all the constraints being dynamic, described in Sec. 5.1.2, in Fig. 2.

Note that even if the physical quantities turns out the same, left and mid
plots of Figs. 1 and 2, the convergence for the constraints seen in the right
plots are qualitatively different. In order to reach even higher accuracy for
the angular momentum constraint $G_{3}$, one need a smaller $\Delta x$.

Figure 2: (Color online) Stationary solution of (101), (103) with fully
dynamic constraints, see Sec. 5.1.2. Everything is the same as for Fig. 1,
except that we here have $3$ damping constants for the $3$ dynamic constraints
(100), i.e., $\eta=k_{1}=k_{2}=k_{3}=2$, and a qualitatively different
convergence of the constraints, $G_{1},G_{2}$, see the right plot.

####  5.1.4 Computational complexity for $n\leq 10$ components

Theoretical and experimental studies of physical systems described by the NLSE
with more than two components is a rapidly growing field, see, e.g., [7, 9]
and references therein. In this light, efficient numerical methods to treat
such systems for any physical parameters, as we have formulated in this
article, are timely. In this section we evaluate numerically the computational
complexity of the presented dynamic methods with up to ten components.

If one or several components represent a very small part of the total density,
they will not influence the other components much, which may result in
numerical results not showing general behaviour of our method but rather
numerical artifacts. Hence, given that we want to test up to $n\leq 10$
components numerically here, we set a lower bound for each density to
$1/(2n)$, such that no component will contain less than $5\%$ of the total
density.

Another possible qualitative problem when evaluating different number of
components with similar parameter values that we want to decrease, is that
many aspects of the solutions are dictated by number theoretical effects, as
have been discussed for $n=2$ components in [21]. Therefore we add to
component $k$ a term $a^{k}$ from a geometric series with $n$ terms and sum
$1-n\cdot 1/(2n)=1/2$, such that
$\sum_{k=1}^{n}c_{k}=\sum_{k=1}^{n}1/(2n)+a^{k}=1$, i.e., $0<a<1$ is a (real)
solution to $a^{n+1}+3a/2+1/2=0$. Then we obtain for example: $c_{1}=1$ for
$n=1$; $c_{1}=0.6160,\>c_{2}=0.3840$ for $n=2$; and
$c_{1}=0.5092,\>c_{2}=0.2840,\>c_{3}=0.2068$ for $n=3$.

The strength of the parameters for the nonlinear coupling were here set to
$\gamma_{ii}=10$ and $\gamma_{ij}=5,\>i\neq j$ (which do not give spatial
separations for the different components).

Initial conditions used for the calculations in all Figs. 1, 2 and 3 are
$\Psi_{j}(x,0)=$ $\sqrt{c_{j}/(4\pi)}(1+\exp(ix))$, for each component
$j=1,2,\ …,n$, such that $\langle\vec{\Psi}\mid\vec{\Psi}\rangle=1$ and
$\langle\vec{\Psi}\mid L\vec{\Psi}\rangle=0.5$ for $\tau=0$. As is natural, we
observed that the convergence is faster for calculations with $\ell_{0}$ close
to the initial value $\ell=0.5$, but the number of iterations that we use to
calculate the computational complexity is the average of the $9$ different
values $\ell_{0}=0.1,\ 0.2,\ …,0.9$. For one case of the data reported in the
green curve of Fig. 3, we did not reach any convergence at all for the above
described choices, which also represents a specific indication of our general
obtained experience that the SHAKE and RATTLE [23, 24] method is further more
sensitive to the initial conditions than the methods with dynamic constraints.
Instead, for this case we both reduced $\Delta\tau=0.9\Delta x$ and used
$\Psi_{2}(x,0)=\sqrt{c_{2}/(4\pi)}(1+\exp(1.1ix))$ in order to reach
convergence.

In addition to the obtained (averaged over $c_{n+1}=\ell_{0}$) number of
iterations $k$ needed for the norm
$||\vec{\Psi}_{k}-\vec{\Psi}_{k-1}||<10^{-6}$, we multiplied with an overall
factor $n$ describing the complexity in setting up the equations to obtain the
Lagrange parameters of each method, see Eqs. (77), and (83), (88),
respectively. For the SHAKE and RATTLE we use an additional factor $4$ coming
from the additional (half-) step in updating the numerical solution together
with $3$ Newton iterations, which was the observed minimum number of
iterations in the inner solver for the Lagrange parameters. The time step in
use for the benchmarking in Fig. 3, is the seemingly maximal possible
$\Delta\tau<\Delta x$ (we set $\Delta\tau=0.99\Delta x$) However, we have not
been able to justify this theoretically. The damping parameters are in all
cases set to $\eta=k_{j}=2$, which gives good performance, but we stress that
we did not perform any individual optimization of the damping parameters for
each case. For the SHAKE and RATTLE method [23, 24] one needs in addition
initial conditions for the Lagrange parameters, which we set to
$\lambda_{j}=1$ with $j=1,2$ for $n=1$, and $j=1,2,3$ for $n=2$. The
substantially larger errorbars for some of the $n$-values in Fig. 3, is due to
number theoretical effects [21], even if we have chosen parameters as to
reduce such effects.

Figure 3: (Color online) Benchmarking of the computational complexity for
constrained $n$-component solutions of (98). The solid (blue) data is for Only
dynamically damped constraints (Sec. 5.1.2), dashed (red) data (slightly
translated to the right of the integer values $n$ for improved visibility) is
for Projection of all the normalization constraints (Sec. 5.1.1), dashed-
dotted (green) data is for SHAKE and RATTLE [24] with available algorithms for
$n=1$ and $n=2$ taken from [22]. The computational complexity for each method
has been calculated by averaging over different $\ell$ values as described in
the main text, and the mean (curves) and standard deviation (half the total
length of the errorbars) are formed from the number of iterations until
convergence, multiplied with different factors describing the three different
algorithms under study here, see main text.

###  5.2 NLSE in 2D under rotation in a harmonic oscillator potential,
applications to vortex structures

Theoretical and experimental studies of quantized vortices in, e.g., liquid
Helium, superconductors, and cold atomic gases, have been undertaken for many
years. Solving the NLSE under rotation is a common but non-trivial such
example, see, e.g., [12, 25, 26, 27, 28, 29] and references therein.

In the 2D examples presented here we want to find $\Psi(x_{1},x_{2})$ where

$E\left(\Psi\right)=\frac{1}{2}\int\limits_{\mathbb{R}^{2}}\left|\dfrac{\partial\Psi}{\partial
x_{1}}\right|^{2}+\left|\dfrac{\partial\Psi}{\partial
x_{2}}\right|^{2}+\left(x_{1}^{2}+x_{2}^{2}\right)\left|\Psi\right|^{2}+g\left|\Psi\right|^{4}\,{\rm
d}x_{1}{\rm d}x_{2},$ (104)

gives the stationary equation

$\frac{\delta
E(\Psi)}{\delta\Psi^{\dagger}}=-\dfrac{1}{2}\left(\frac{\partial^{2}\Psi}{\partial
x_{1}^{2}}+\frac{\partial^{2}\Psi}{\partial
x_{2}^{2}}\right)+\dfrac{1}{2}\left(x_{1}^{2}+x_{2}^{2}\right)\Psi+g\left|\Psi\right|^{2}\Psi=0.$
(105)

The dynamic constraints are

$G_{1}=c_{1}-\int\limits_{\mathbb{R}^{2}}\left|\Psi\right|^{2}\,{\rm
d}x_{1}{\rm d}x_{2}=0,\
G_{2}=c_{2}-\int\limits_{\mathbb{R}^{2}}\Psi^{\dagger}L\Psi\,{\rm d}x_{1}{\rm
d}x_{2}=0,$ (106)

with the angular momentum operator

$L=-i(x_{1}\frac{\partial}{\partial x_{2}}-x_{2}\frac{\partial}{\partial
x_{1}}).$ (107)

We can generally attain the corresponding Lagrange parameters $\vec{\lambda}$
from the linear system given by the matrix $A$ in (83) and the vector $b$ in
(88) by solving $A\vec{\lambda}=b$.

Figure 4: (Color online) Stationary solution of (108), density (left), and
phase (right), of $\Psi(x_{1},x_{2})$. Compare with Figure 1 of [12]. We used
the damping $\eta=4$ and $\Delta\tau=0.1$.

####  5.2.1 The groundstate of a rapidly rotating Bose-Einstein condensate

To compare our numerical results in 2D with the literature, we first
reproduced "A simple but complete example" from the dimensionless GPE in Sec.
6 of [12]. In this example the rotational velocity $\Omega$ is kept constant,
so we used the following dynamical reformulation of (105)

$\frac{\partial^{2}\Psi}{\partial\tau^{2}}+\eta\frac{\partial\Psi}{\partial\tau}=\frac{1}{2}\left(\frac{\partial^{2}\Psi}{\partial
x_{1}^{2}}+\frac{\partial^{2}\Psi}{\partial
x_{2}^{2}}\right)-\frac{1}{2}\left(x_{1}^{2}+x_{2}^{2}\right)\Psi-500\left|\Psi\right|^{2}\Psi+\lambda_{1}\Psi+\lambda_{2}L\Psi,$
(108)

where $\lambda_{2}=\Omega$ (constant),
$\lambda_{1}=\mu=E_{\textnormal{{rot}}}+\int\int
250\left|\Psi\right|^{4}+\Psi^{\dagger}\left(\frac{\partial^{2}\Psi}{\partial\tau^{2}}+\eta\frac{\partial\Psi}{\partial\tau}\right)dx_{1}dx_{2}$,
and with projection $\Psi=\Psi/\left\|\Psi\right\|$ in each time step to
fulfill the normalization constraint. Here $E_{\textnormal{{rot}}}$ is the
energy in the rotating frame

$E_{\textnormal{{rot}}}\left(\Omega\right)=\int_{\mathbb{R}^{2}}\frac{1}{2}\left|\nabla\Psi\right|^{2}+\frac{1}{2}\left(x_{1}^{2}+x_{2}^{2}\right)\left|\Psi\right|^{2}+250\left|\Psi\right|^{4}-\Omega\Psi^{*}L\Psi\,dx_{1}dx_{2}.$
(109)

As in sec. 6 of [12] we use a spatial grid $-10\leq x_{1},x_{2}\leq 10$ with
$129$ points in each dimension, and the initial condition [Eq. (3.23) of [12],
with $\gamma_{x}=\gamma_{y}=1$],

$\Psi\left(x_{1},x_{2},0\right)=\frac{\left[\left(1-\Omega\right)+\Omega\left(x_{1}+ix_{2}\right)\right]\exp\left(-x_{1}^{2}/2-x_{2}^{2}/2\right)}{\left\|\left[\left(1-\Omega\right)+\Omega\left(x_{1}+ix_{2}\right)\right]\exp\left(-x_{1}^{2}/2-x_{2}^{2}/2\right)\right\|}.$
(110)

Performing the numerical calculation solving (108), we obtain for $\Omega=0.5$
a stationary solution with a density $\rho=\left|\Psi\right|^{2}$, which is
5-fold symmetric, depicted in Fig. 4.

Figure 5: (Color online) Nucleation of a central vortex studied with dynamical
constraints. The upper panels shows the density
$\rho=\left|\Psi(x_{1},x_{2})\right|^{2}$ for angular momentum
$c_{2}=\ell=0.2,\>0.8,\>1$, and the middle panels show the phase of $\Psi$,
for stationary solutions of (108) with $\mu$ and $\Omega$ being Lagrange
parameters. In the lower panels the corresponding numerical convergence of the
DFPM and the dynamical constraints are illustrated. The (black) dashed curve
shows $\left\|\Psi_{k}-\Psi_{k-1}\right\|$ for iteration step $k$, the (blue)
solid curve shows the normalization constraint
$\left|G_{1}\left(\Psi_{k}\right)\right|$, and the (green) dotted-dashed curve
shows the angular momentum constraint
$\left|G_{2}\left(\Psi_{k}\right)\right|$. For this figure we used the
numerical parameters $\eta=k_{1}=k_{2}=2$ and $\Delta\tau=0.05$.

####  5.2.2 Transition of a non rotating Bose-Einstein condensate forming a
central vortex

After the first test of DFPM without dynamical constraints above, we now test
the method with two dynamical constraints on a well-known example, the
transition of a distant vortex approaching the center
$\left(x_{1},x_{2}\right)=\left(0,0\right)$ for $0<\ell\leq 1$. Compare our
Fig. 5, for example with Figure 3 in [28], or Fig. 3 in [27]. This means we
are now using Eq. (108) with the two Lagrange multipliers evaluated
dynamically, corresponding to Eqs. (91), (94), and the same initial condition,
Eq. (110) with $\Omega=0.1$, for all the three different constrained
$\ell$-values in Fig. 5. The densities (upper plots) and phases (middle) of
the stationary solution to (108) behaves as expected, and the convergence of
the iteration and the constraints (106) are shown in the lower plots.

In conclusion we find the presented method with dynamic constraints an
effective tool to calculate vortex structures in 2D with no assumptions on the
strength of the nonlinearity or the form of the external potential. For
example we have also successfully tested to calculate giant vortices for an
additional quartic term in the potential [25], and rotating states for an
asymmetric potential [26]. We have kept the damping parameters and the time
step constant throughout this subsection to illustrate the robustness of the
method. However, we noted that the numerical performence can be substantially
improved by finetuning those numerical parameters.

##  6 Conclusions

We have set up a general dynamical formulation Eqs. (17), (18), (24), (25) for
obtaining stationary solutions for nonlinear equations with $n_{c}$ global
constraints, where $0\leq n_{P}\leq n_{c}$ of them are projected and
$n_{c}-n_{P}$ of them are included in the (extended) dynamical formulation.

We presented detailed derivations for the case of rotation in 1D system, with
examples of treating the constraints only with projection, Sec. 4.1,
projection of some, Sec. 4.2, and all, Sec. 4.3, normalization constraints,
and finally with only dynamic constraints, Sec. 4.4.

In the final part of the article we have given detailed formulae for some 1D,
Figs. 1, and 2, and 2D, Figs. 4, and 5, examples that we have evaluated
numerically and compared to the literature. We also evaluated the
computational performance numerically with up to ten components, Fig. 3.

In conclusion, the dynamical formulation of complicated global constraints
have the following main advantages:

  * •

Simpler implementations

  * •

Less computational complexity

  * •

Less sensitivity to initial conditions

We do not report any disadvantages with the methods studied. It may appear as
that we need to chose additional suitable numerical parameter values (the
$k_{j}$:s in (18)) for the damping of the constraints. However, with
projection methods one needs instead to chose suitable initial values for the
Lagrange parameters.

## References

  * [1] Fibich G 2015, “The Nonlinear Schrödinger Equation”, Springer, New York. ISBN 978-3-319-12747-7.
  * [2] Ögren M, Abdullaev F Kh, and Konotop V V 2017, ”Solitons in a PT-symmetric $\chi^{(2)}$ coupler”, Opt. Lett. 42, 4079.
  * [3] Sørensen M P, Pedersen N F, and Ögren M 2017, “The dynamics of magnetic vortices in type II superconductors with pinning sites studied by the time dependent Ginzburg-Landau model”, Physica C 533, 40.
  * [4] Chavanis P-H 2011, “Mass-radius relation of Newtonian self-gravitating Bose-Einstein condensates with short-range interactions. I. Analytical results”, Phys. Rev. D 84, 043531.
  * [5] Rozenman G G, Shemer L, and Arie A 2020, “Observations of accelerating solitary wavepackets”, Phys. Rev. E 101, 050201(R).
  * [6] Dysthe K, Krogstad H E, and Müller P 2008, “Oceanic Rogue Waves”, Annu. Rev. Fluid Mech. 40, 287.
  * [7] Feng B, F, 2014, "General N-soliton solution to a vector nonlinear Schrödinger equation", J. Phys. A: Math. Theor. 47, 355203.
  * [8] Beattie S, Moulder S, Fletcher R J, and Hadzibabic Z 2013, “Persistent Currents in Spinor Condensates”, Phys. Rev. Lett. 110, 025301.
  * [9] Lannig S, Schmied C M, Prüfer M, Kunkel P, Strohmaier R, Strobel H, Gasenzer T, Kevrekidis P G, and Oberthaler M K, 2020, "Collisions of Three-Component Vector Solitons in Bose-Einstein Condensates", Phys. Rev. Lett. 125, 170401.
  * [10] Ögren M, Kavoulakis G M, and Jackson A D 2005, “Solitary waves in elongated clouds of strongly-interacting bosons”, Phys. Rev. A 72, 021603(R).
  * [11] Ögren M and Kavoulakis G M 2020, “Rotational properties of superfluid Fermi-Bose mixtures in a tight toroidal trap”, Phys. Rev. A .102, 013323.
  * [12] Antoine X and Duboscq R, 2014, "GPELab, a Matlab Toolbox to Solve Gross-Pitaevskii Equations I: Computation of Stationary Solutions, Computer Physics Communications" 185 (11), pp. 2969-2991.
  * [13] Smyrlis G and Zisis V 2004, “Local convergence of the steepest descent method in Hilbert spaces”, Journal of Mathematical Analysis and Applications , 300, 436.
  * [14] Nocedal J and Wright S J 2006, “Numerical Optimization”, 2nd ed. Springer, ISBN-13 978-0387-30303-1.
  * [15] Sra S, Nowozin S, and Wright S J 2012, "Optimization from Machine Learning", MIT Press, ISBN 9780262016469
  * [16] Boron L F and Zeidler E 1984, ”Nonlinear Functional Analysis and its Applications: III: Variational Methods and Optimization”, Springer New York.
  * [17] Gulliksson M, Ögren M, Oleynik A, and Zhang Y 2018, “Damped Dynamical Systems for Solving Equations and Optimization Problems”. In: Sriraman B. (eds) Handbook of the Mathematics of the Arts and Sciences. Springer, Cham, ISBN 978-3-319-70658-0.
  * [18] Ögren M and Gulliksson M 2020, "A numerical damped oscillator approach to constrained Schrödinger equations", Eur. J. Phys. 41, 065406.
  * [19] Hairer E, Lubich C, and Wanner G 2006, “Geometric Numerical Integration”, 2nd ed. Springer, ISBN 978-3-540-30666-5.
  * [20] Wu Z and Zaremba E 2013, "Mean-field yrast spectrum of a two-component Bose gas in ring geometry: Persistent currents at higher angular momentum", Phys. Rev. A 88, 063640.
  * [21] Roussou A, Smyrnakis J, Magiropoulos M, Efremidis N K, Kavoulakis G M, Sandin P, Ögren M, and Gulliksson M, 2018, "Excitation spectrum of a mixture of two Bose gases confined in a ring potential with interaction asymmetry.", New J. Phys. 20, 045006.
  * [22] Sandin P, Ögren M, and Gulliksson M 2016, “Numerical solution of the stationary multicomponent nonlinear Schrödinger equation with a constraint on the angular momentum”, Phys. Rev. E 93, 033301.
  * [23] Andersen H C 1983, "Rattle: A ’velocity’ version of the shake algorithm for molecular dynamics calculations", Journal of Computational Physics, 52 (1):24 – 34.
  * [24] McLachlan R, Modin K, Verdier O, and Wilkins M 2014, “Geometric generalisations of SHAKE and RATTLE”, Foundations of Computational Mathematics 14, 339.
  * [25] Jackson A D, Kavoulakis G M, and Lundh E, 2004, "Phase diagram of a rotating Bose-Einstein condensate with anharmonic confinement", Phys. Rev. A 69, 053619.
  * [26] Linn M, Niemeyer M, and Fetter A L, 2001, "Vortex stabilization in a small rotating asymmetric Bose-Einstein condensate", Phys. Rev. A 64, 023602.
  * [27] Kavoulakis G M, Mottelson B, and Pethick C J, 2000, "Weakly interacting Bose-Einstein condensates under rotation", Phys. Rev. A 62, 063605.
  * [28] Butts D A and Rokhsar D S, 1999, "Predicting signatures of rotating Bose-Einstein condensates", Nature 397, 327.
  * [29] Bao W, Wang H, and Markowich P A, 2005, "Ground, symmetric and central vortex states in rotating Bose-Einstein condensates", Communications in Mathematical Sciences, 3(1):57-88.