The $J$-method for the Gross-Pitaevskii eigenvalue problem

TL;DR

This paper analyzes the $J$-method for solving the Gross-Pitaevskii eigenvalue problem, demonstrating its global and local convergence properties, and showcasing its effectiveness in approximating excited states and localized phenomena.

Contribution

It provides a rigorous convergence analysis of the $J$-method in a Hilbert space framework, including a damping modification and applications to complex physical states.

Findings

Proves global convergence for a damped $J$-method.

Establishes local linear convergence near eigenfunctions.

Demonstrates effectiveness in numerical experiments with localized states.

Abstract

This paper studies the $J$-method of [E. Jarlebring, S. Kvaal, W. Michiels. SIAM J. Sci. Comput. 36-4:A1978-A2001, 2014] for nonlinear eigenvector problems in a general Hilbert space framework. This is the basis for variational discretization techniques and a mesh-independent numerical analysis. A simple modification of the method mimics an energy-decreasing discrete gradient flow. In the case of the Gross-Pitaevskii eigenvalue problem, we prove global convergence towards an eigenfunction for a damped version of the $J$-method. More importantly, when the iterations are sufficiently close to an eigenfunction, the damping can be switched off and we recover a local linear convergence rate previously known from the discrete setting. This quantitative convergence analysis is closely connected to the~$J$-method's unique feature of sensitivity with respect to spectral shifts. Contrary to classical gradient flows, this allows both the selective approximation of excited states as well as the amplification of convergence beyond linear rates in the spirit of the Rayleigh quotient iteration for linear eigenvalue problems. These advantageous convergence properties are demonstrated in a series of numerical experiments involving exponentially localized states under disorder potentials and vortex lattices in rotating traps.

Equations

$$\mathcal{A}(v)=A(v,v),$$

$$\left[-\Delta +W-\Omega {\mathcal{L}}_z+\kappa |u^{\ast }{|}^2\right]u^{\ast }={\lambda }^{\ast }{u}^{\ast }.$$

$$\mathcal{A}:V\to V^{\ast }$$

$$\mathcal{A}(v)=A(v,v),$$

$$A(v,w_1+w_2)=A(v,w_1)+A(v,w_2)\text{mbox}andA(v,\alpha w_1)=\alpha A(v,w_1).$$

$$A(u,v)=A(\alpha u,v)$$

$$(\cdot ,\cdot {)}_H:H\times H\to \mathbb{R}.$$

$$⟨\mathcal{I}u,\cdot {⟩}_{V^{\ast },V}=(u,\cdot {)}_H$$

$$⟨\mathcal{A}(u^{\ast }),v⟩=⟨{\lambda }^{\ast }\mathcal{I}{u}^{\ast },v⟩={\lambda }^{\ast }(u^{\ast },v{)}_H\text{for all test functions }v\in V.$$

$$(v,w{)}_H:=\Re \left({\int }_{\mathcal{D}}v\overline{w}\text{d}x\right)$$

$$(v,w{)}_V:=\Re \left({\int }_{\mathcal{D}}\nabla v\cdot \overline{\nabla w}\text{d}x\right).$$

$$⟨A(u,v),w⟩:=\Re \left({\int }_{\mathcal{D}}{\nabla }_{\text{mbox}\mathrm{R}}v\cdot \overline{{\nabla }_{\text{mbox}\mathrm{R}}w}+W_{\text{mbox}\mathrm{R}}v\overline{w}\text{d}x+\frac{\kappa }{\Vert u{\Vert }^2}{\int }_{\mathcal{D}}|u{|}^{2}v\overline{w}\text{d}x\right).$$

$${\nabla }_{\text{mbox}\mathrm{R}}v:=\nabla v+\mathrm{i}\frac{\Omega }{2}{R}^{\top }v$$

$$\underset{h\in X\setminus \{0\},h\to 0}{\lim }\frac{\Vert F(u+h)-F(u)-F^{\prime }(u)h{\Vert }_Y}{\Vert h{\Vert }_X}=0.$$

$$\underset{t\to 0}{\lim }\frac{\Vert A(u+tu,\cdot )-A(u,\cdot )-0{\Vert }_{\mathcal{L}(V,V^{\ast })}}{t\Vert u{\Vert }_V}=\underset{t\to 0}{\lim }\frac{\Vert A(u,\cdot )-A(u,\cdot ){\Vert }_{\mathcal{L}(V,V^{\ast })}}{t\Vert u{\Vert }_V}=0.$$

$${\mathcal{A}}^{\prime }(u;w)$$

$$=\underset{t\to 0}{\lim }\frac{A(u+tw,u)-A(u,u)}{t}+\underset{t\to 0}{\lim }\frac{A(u+tw,u+tw)-A(u+tw,u)}{t}$$

$$={\partial }_{1}A(u,u;w)+A(u,w).$$

$${\mathcal{A}}^{\prime }(u;u)={\partial }_{1}A(u,u;u)+A(u,u)=A(u,u)=\mathcal{A}(u).$$

$$J(u):={\mathcal{A}}^{\prime }(u;\cdot ).$$

$$⟨J(u^{\ast })u^{\ast },v⟩={\lambda }^{\ast }(u^{\ast },v{)}_H\text{mbox}forallv\in V.$$

$$⟨(J(v)+\sigma \mathcal{I})u,w⟩=⟨F,w⟩$$

$$\psi (v):=(J(v)+\sigma \mathcal{I}{)}^{-1}\mathcal{I}v.$$

$$\varphi (v):=\frac{\psi (v)}{\Vert \psi (v)\Vert }.$$

$$({\lambda }^{\ast }+\sigma )\psi (u^{\ast })$$

$$=(J(u^{\ast })+\sigma \mathcal{I}{)}^{-1}(J(u^{\ast })+\sigma \mathcal{I})u^{\ast }=u^{\ast }.$$

$$u^{n+1}=\varphi (u^n)=\frac{\psi (u^n)}{\Vert \psi (u^n)\Vert }={\alpha }_n(J(u^n)+\sigma \mathcal{I}{)}^{-1}\mathcal{I}{u}^n$$

$${\tilde{u}}^{n+1}=(1-{\tau }_n)u^n+{\tau }_n{\gamma }_n{J}_{\sigma }(u^n{)}^{-1}\mathcal{I}{u}^n,u^{n+1}=\frac{{\tilde{u}}^{n+1}}{\Vert {\tilde{u}}^{n+1}\Vert }$$

$$({\tilde{u}}^{n+1}-u^n,u^n{)}_H=0.$$

$$1=\Vert u^n{\Vert }^2=(u^n,{\tilde{u}}^{n+1}{)}_H\le \Vert {\tilde{u}}^{n+1}\Vert .$$

$$u^{n+1}:=\varphi (u^n)$$

$$\Vert u^n-u^{\ast }{\Vert }_V\le C_{\epsilon }|{\rho }^{\ast }+\epsilon {|}^n\Vert u^0-u^{\ast }{\Vert }_V\text{mbox}forall{u}^0\in S_{\epsilon }$$

Full text

P. Henning acknowledges funding by the Swedish Research Council (grant 2016-03339). The work of D. Peterseim is part of a project that has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 865751).

1. Introduction

This paper studies eigenvalues of nonlinear differential operators $\mathcal{A}$ on a real Hilbert space $V$ of the form

[TABLE]

where $A\colon V\times V\rightarrow V^{*}$ is a continuous, bounded mapping that is invariant under scaling of the first argument and real-linear in the second argument. A well-known example that can be cast in such a format is the Gross-Pitaevskii eigenvalue problem (GPEVP) in a bounded domain $\mathcal{D}\subset\mathbb{R}^{d}$ (with $d=2,3$). In the strong form the GPEVP reads

[TABLE]

Here, $W\in L^{\infty}(\mathcal{D},\mathbb{R})$ denotes a non-negative and space-dependent potential, $\Omega\in\mathbb{R}$ is the angular velocity, $\mathcal{L}_{z}=-\mathrm{i}\left(x\partial_{y}-y\partial_{x}\right)$ is the $z$-component of the angular momentum, and $\kappa\geq 0$ regulates the nonlinearity of the problem. This problem is related, e.g., to the modeling of Bose-Einstein condensates [Bos24, Ein24, DGPS99, PS03]. In this context, a solution $u^{*}$ represents a stationary quantum state of the condensate, $|u^{*}|^{2}$ is the corresponding density, and $\lambda^{*}$ the so-called chemical potential.

The numerical solution of the GPEVP (1) has been studied extensively in recent years. Popular methods for solving the nonlinear eigenvalue problem are for example the Self Consistent Field Iteration (SCF), cf. [Can00, CLB00a, CLB00b, DC07, UJR18], which requires to solve a linearized eigenvalue problem in each iteration. Algorithms that belong to the SCF class are for instance the Roothaan algorithm [Roo51] or the optimal damping algorithm proposed in [CLB00a]. Another important class of iterative methods is based on gradient flows for the energy functional associated with (1), where we mention the Discrete Normalized Gradient Flow (DNGF), cf. [BD04, BCL06, BS08, BWM05], which is based on an implicit Euler discretization of the $L^{2}$-gradient flow. Improvements of the approach by using conjugated gradients were proposed in [ALT17]. Here we also mention the Projected Sobolev Gradient Flows (PSGFs), cf. [DK10, GRPG01, HSW19, HP20, KE10, RSB14, RSSL09, Zha19], which form a subclass of the gradient flow methods. Sobolev gradients are the Riesz representants of the Fréchet derivative of the energy functional in a suitable Hilbert space. With this, PSGFs are based on computing such Sobolev gradients and then projecting them onto the tangent space of the normalization constraint. An improvement of PSGF by using Riemannian conjugate gradients was suggested in [DP17]. Another strategy to solve the GPEVP involves a direct minimization of the energy functional, cf. [BT03, CORT09], which means that (1) is written as a nonlinear saddle point problem which is then solved by a Newton-type method. Global convergence results for solving the GPEVP are very rare in the literature. To the best of our knowledge, global convergence had been so far only established for a damped PSGF suggested in [HP20], where corresponding analytical results can be found in [HP20, Zha19].

The SCF and PSGF approaches are typically based on the simple linearization of the partial differential operator by replacing the density $|u^{*}|^{2}$ in (1) by the density of some given approximation of $u^{*}$. In the finite-dimensional case (after spatial discretization) these approaches can be interpreted as generalizations of the inverse iteration (inverse power method). Linear convergence is observed empirically which is in agreement with the convergence theory for linear matrix eigenvalue problems. However, in the presence of clustered eigenvalues, which are typically connected to interesting physical phenomena, linear convergence can be very slow and shifting (as in the shifted inverse iteration or the Rayleigh quotient iteration) is a well-established technique for the acceleration of matrix eigenvalue solvers. The problem is that the aforementioned schemes seem to be unable to achieve any speed-up by a sophisticated shift of the nonlinear operator. In [JKM14], Jarlebring, Kvaal, and Michiels observed that this insensitivity towards spectral shifts is the result of an unsuitable choice of linearization. The authors propose a natural but conceptually different linearization using the Jacobian of the nonlinear operator $\mathcal{A}$. The resulting approach does not belong to any of the classes above and we will refer to it as the $J$-method.

This paper generalizes the $J$-method and its numerical analysis to an abstract Hilbert space setting. Hence, the resulting iteration scheme is based on a variational formulation which is mesh-independent (cf. [AF20] for a similar approach in electromagnetism). This then allows the application of any spatial discretization, including finite elements and spectral methods. As discovered in  [JKM14], the correct choice of the linearization in the $J$-method does not only lead to a competitive nonlinear eigenvalue solver, it also allows major theoretical progress such as a proof of local linear convergence based on a version of Ostrowski’s theorem (see Section 3 for the corresponding result the Hilbert space setting). The obtained linear convergence rate, which is essentially $|\lambda^{*}+\sigma|/|\mu+\sigma|$, depends on the spectral gap of the $\sigma$-shifted Jacobian, where $\mu+\sigma$ is the second-smallest eigenvalue of the $\sigma$-shifted Jacobian around $u^{*}$. Our results generalize the findings of [JKM14] from the matrix case to the abstract setting. In the case of the GPEVP, which is discussed in Section 4, this marks the first quantified convergence result. Moreover, an adaptive choice of the shifts in the spirit of the Rayleigh quotient iteration amplifies convergence beyond the linear rate in representative numerical experiments, see Section 6. At the same time, the sensitivity to spectral shifts facilitates the computation of excited states.

As usual for nonlinear problems, the quantitative convergence results are of local nature in the sense that they require a sufficiently accurate initial approximation. In practical computations this may be a rather challenging task. A simple damping strategy inspired by an energy-dissipative gradient flow [HP20] resolves this problem. In the case of the GPEVP we prove global convergence of the method towards an eigenfunction with a guaranteed decrease of energy, cf. Section 5.

Altogether, the combination of the $J$-method with damping for globalization and shifting for acceleration provides a powerful methodology for the simulation and analysis of nonlinear PDE eigenvector problems such as the GPEVP.

2. Definition of the $J$-method in Hilbert spaces

2.1. Nonlinear eigenvector problem

We consider a nonlinear differential operator

[TABLE]

that maps a real Hilbert space $V$ into its dual space $V^{*}$. In particular, $V$ is equipped with an $\mathbb{R}$-inner product and every $F\in V^{*}$ satisfies $\langle F,v\rangle_{V^{*},V}\in\mathbb{R}$ for all $v\in V$. Note that $V$ is still allowed to contain complex-valued functions. We assume that $\mathcal{A}$ has the form

[TABLE]

where $A\colon V\times V\rightarrow V^{*}$ is continuous, bounded, and (real-)linear in the second argument. Recall that $A(v,\,\cdot\,)$ being real-linear (cf. [BCS17, Def. 3.50]) means that for all $v,w_{1},w_{2}\in V$ and $\alpha\in\mathbb{R}$ we have

[TABLE]

Furthermore, we assume that $A(\,\cdot\,,\cdot\,)$ is sufficiently smooth on $V\setminus\{0\}\times V$ (in particular real-Fréchet differentiable in both arguments) and that it is real-scaling invariant in the first argument, i.e.,

[TABLE]

for all $u,v\in V$ and all $\alpha\in\mathbb{R}\setminus\{0\}$.

For the weak formulation of the eigenvalue problem with a nonlinearity in the eigenfunction, we introduce another real Hilbert space $H$ (the pivot space) with $\mathbb{R}$-inner product

[TABLE]

such that $V,H,V^{*}$ forms a Gelfand triple [Zei90, Ch. 23.4]. The corresponding $H$-norm $\|\cdot\|:=\|\cdot\|_{H}$ will be used for the normalization condition of the eigenfunction. The corresponding embedding $\mathcal{I}\colon V\to V^{*}$ is defined via

[TABLE]

for $u\in V$. For brevity, we shall write $\langle\,\cdot\,,\,\cdot\rangle$ for the canonical duality pairing in $V$.

The goal of this paper is to solve the corresponding PDE eigenvalue problem: find an eigenfunction $u^{*}\in V$ with $\|u^{*}\|=1$ and an eigenvalue $\lambda^{\ast}\in\mathbb{R}$ such that $\mathcal{A}(u^{*})=\lambda^{*}\mathcal{I}u^{*}$. This is equivalent to the variational formulation

[TABLE]

Throughout the paper, we denote by $\overline{\alpha}$ the complex conjugate of a complex number $\alpha\in\mathbb{C}$, by $\Re(v)$ the real part of $v$, and by $\Im(v)$ the imaginary part.

Example 2.1**.**

We are particularly interested in the GPEVP (1). Here we consider $H=L^{2}(\mathcal{D}):=L^{2}(\mathcal{D},\mathbb{C})$ as a real Hilbert space equipped with the inner product

[TABLE]

and analogously $V=H^{1}_{0}(\mathcal{D}):=H^{1}_{0}(\mathcal{D},\mathbb{C})$ as a real Hilbert space equipped with

[TABLE]

The (real) dual space is denoted by $V^{*}=H^{-1}(\mathcal{D}):=H^{-1}(\mathcal{D},\mathbb{C})$. The corresponding nonlinear operator reads

[TABLE]

Here, we have $W_{\mbox{\rm\tiny R}}(x)=W(x)-\tfrac{1}{4}\,\Omega^{2}\,|x|^{2}$ and the rotational gradient $\nabla_{\mbox{\rm\tiny R}}$ is given by

[TABLE]

for the divergence-free vector field $R(x,y,z):=(y,-x,0)$ if $d=3$ and $R(x,y):=(y,-x)$ if $d=2$. Note that the nonlinear term is multiplied by $\|u\|^{-2}$ in order to achieve the assumed scaling invariance of $A$ in the first component.

2.2. Linearization and the $J$-operator

Recall that the real-Fréchet derivative of some $F\colon X\rightarrow Y$ in a point $u\in X$ is a bounded and $\mathbb{R}$-linear operator $F^{\prime}(u)\colon X\rightarrow Y$ with the property

[TABLE]

We write $F^{\prime}(u;h):=F^{\prime}(u)h$ for the Fréchet derivative in $u$ in direction of $h$. In the following, we speak about the Fréchet derivatives of operators, where we always mean the real-Fréchet derivative as stated above. Moreover, we neglect writing the supplement $h\in X\setminus\{0\}$ beneath the limit. For the operator $A$, we denote the partial Fréchet derivative with respect to the (nonlinear) first component by $\partial_{1}A\colon V^{3}\rightarrow V^{*}$. Because of the scaling invariance of $A$, i.e., $A(u,\cdot\,)=A(\alpha u,\cdot\,)$, we have $\partial_{1}A(u,\cdot\,;u)=0$, which can be seen by

[TABLE]

The Fréchet derivative of $\mathcal{A}$ in $u\in V$ and in direction $w\in V$ can be expressed as

[TABLE]

Due to $\partial_{1}A(u,\cdot\,;u)=0$, we conclude

[TABLE]

For an element $u\in V$ we now define the operator $J(u)\colon V\rightarrow V^{*}$ by

[TABLE]

Hence, using $\langle J(u^{\ast})u^{\ast},v\rangle=\langle\mathcal{A}^{\prime}(u^{\ast},u^{\ast}),v\rangle=\langle\mathcal{A}(u^{\ast}),v\rangle$, we observe that the eigenvalue problem (2) can be rewritten as: find $u^{\ast}\in V$ with $\|u^{\ast}\|=1$ and $\lambda^{\ast}\in\mathbb{R}$ such that

[TABLE]

2.3. The shifted $J$-method

The shifted $J$-method is defined by applying the inverse power iteration to the reformulated eigenvalue problem (5) based on the linearization $J$ and some spectral shift $\sigma$. For that, we consider a function $v\in V$ and $\sigma\in\mathbb{R}$ such that $J(v)+\sigma\mathcal{I}\colon V\to V^{*}$ is invertible. In practice, one may either choose $\sigma$ large such that $J(v)+\sigma\mathcal{I}$ is coercive (see Section 5.1 for the example of the GPEVP) or close to the negative of the target eigenvalue (as it is done in the numerical experiments of Section 6).

For $F\in V^{*}$ we define $u=(J(v)+\sigma\mathcal{I})^{-1}F\in V$ as the unique solution of the variational problem

[TABLE]

for all $w\in V$. Considering a fixed shift $\sigma$ such that $-\sigma$ is no eigenvalue of $J(v)$, we define $\psi\colon V\to V$ by

[TABLE]

Including an additional normalization step finally leads to the operator $\phi\colon V\to V$,

[TABLE]

Now consider a normalized eigenpair $(u^{*},\lambda^{*})$ of the eigenvalue problem (2). Then, $\psi(u^{*})$ satisfies due to $J(u^{*})u^{*}=\mathcal{A}(u^{*})=\lambda^{*}\mathcal{I}u^{*}$,

[TABLE]

As a consequence, we observe that $u^{\ast}$ is a fixed point of $\phi$, i.e., $u^{\ast}=\phi(u^{\ast})$. This motivates the following iteration scheme: given $u^{0}\in V$ with $\|u^{0}\|=1$ and a shift $\sigma$ such that $J(u^{0})+\sigma\mathcal{I}$ is invertible, compute iterates $u^{n+1}$ for $n\geq 0$ by

[TABLE]

with the normalization factor $\alpha_{n}=1/\|(J(u^{n})+\sigma\mathcal{I})^{-1}\mathcal{I}u^{n}\|$. For the matrix case, the convergence of this scheme was analyzed in [JKM14]. We emphasize that this scheme is well-defined if the shift $\sigma$ guarantees the invertibility of $J(u^{n})+\sigma\mathcal{I}$. This property has to be checked within the specific application and will be proved for the GPEVP in Section 5.1.

2.4. The damped $J$-method

In many applications, eigenvalue problems (with a nonlinearity in the eigenfunction) can be equivalently formulated as a constrained energy minimization problem. In this case, the radius of convergence can be considerably enhanced by introducing a variable damping parameter $\tau_{n}$. This damping parameter (or step size) is adaptively selected such that the energy is optimally minimized in each iteration. In [HP20], a geometric justification of this approach was given, by interpreting it as the discretization of a certain projected Sobolev gradient flow, where the inner product (with which respect the Sobolev gradient is computed) is based on a repeated linearization of the differential operator. Even though $\langle(J(v)+\sigma\mathcal{I})\,\cdot\,,\cdot\,\rangle$ typically does not define an inner product (and hence does not allow a geometric interpretation), it is still possible to generalize the $J$-method defined in (6) in the spirit of the results in [HP20]. In Section 5 below, we shall give a rigorous justification of this approach by proving global convergence of the damped $J$-method in the context of the GPEVP.

The damping strategy aims to find a (optimized) linear combination of $u^{n}$ and $(J(u^{n})+\sigma\mathcal{I})^{-1}\mathcal{I}u^{n}$. For brevity we introduce $J_{\sigma}(u):=J(u)+\sigma\mathcal{I}$ and assume for a moment that the shift $\sigma$ is chosen such that $J_{\sigma}(u^{n})$ remains coercive.

One step of the damped $J$-method then reads as follows: given $u^{n}\in V$ with $\|u^{n}\|=1$ and a step size $\tau_{n}>0$ compute

[TABLE]

with $\gamma_{n}^{-1}:=(J_{\sigma}(u^{n})^{-1}\mathcal{I}u^{n},u^{n})_{H}$. This choice provides the important $H$-orthogonality

[TABLE]

Due to the normalization in (7), we recover the original $J$-method (6) for $\tau_{n}=1$. Further, the $H$-orthogonality (8) implies that the mass (i.e., the $H$-norm) of $\tilde{u}^{n+1}$ is greater or equal to $1$ (even strictly larger unless ${u}^{n+1}=u^{n}$). This can be seen from

[TABLE]

For the example of the GPEVP we will show in Section 5 that the iteration (7) is well-defined and that it converges to an eigenstate if the step size $\tau_{n}$ is sufficiently small.

To summarize, we have introduced a damped version of the $J$-method (7) (with typical choice $\sigma=0$ or $\sigma$ large) which intends to ensure global convergence and a shifted version (6) to accelerate the convergence (with typical choice $\sigma$ close to $-\lambda^{\ast}$). Note, however, that we do not intend to use damping and shifting simultaneously as this seems hard to control numerically. Instead, we propose to apply damping for globalization of convergence and then switch to shifting if a certain accuracy is obtained. This practice then provides a powerful methodology as illustrated in Section 6.

3. Abstract local convergence of the shifted $J$-method

The sensitivity of the $J$-method to spectral shifts facilitates the numerical approximation of excited states. To see this, we transfer the local convergence result presented in [JKM14] to the Hilbert space setting. This shows that the shifted $J$-method converges to an eigenfunction $u^{\ast}$ of $\mathcal{A}$ if the starting function and the shift are sufficiently close to the eigenpair $(u^{\ast},\lambda^{\ast}$) of interest. In this section, we consider the undamped version of the $J$-method, i.e., we consider the iteration $u^{n+1}=\phi(u^{n})$, including an arbitrary shift $\sigma\in\mathbb{R}$ and assuming that $J_{\sigma}(u^{n})$ is invertible. In order to prove local convergence with a linear rate that depends on the shift, we make use of the following well-known lemma that is a version of the Ostrowski theorem in Banach spaces.

Proposition 3.1**.**

Let $\phi\colon V\rightarrow V$ be a mapping on the Banach space $V$ that is Fréchet-differentiable at a point $u^{\ast}\in V$ such that the Fréchet-derivative $\phi^{\prime}(u^{\ast})\colon V\rightarrow V$ is a bounded linear operator with spectral radius $\rho^{\ast}:=\rho(\hskip 2.0pt\phi^{\prime}(u^{\ast})\hskip 2.0pt)<1.$ Then there is an open neighborhood $S$ of $u^{\ast}$, such that for all starting values $u^{0}\in S$ we have that the fixed point iterations

[TABLE]

converge strongly in $V$ to $u^{\ast}$, i.e., $\|u^{n}-u^{\ast}\|_{V}\rightarrow 0$ for $n\rightarrow\infty$. Furthermore, for every $\varepsilon>0$ there exist a neighborhood $S_{\varepsilon}$ of $u^{\ast}$ and a (finite) constant $C_{\varepsilon}>0$ such that

[TABLE]

and $n\geq 1$. Hence, $\rho^{\ast}$ defines the asymptotic linear convergence rate of the fixed point iteration.

Proof.

The proof is based on the observation that under the assumptions of the lemma and for each $\varepsilon>0$, there exists an induced operator norm $\|\cdot\|_{\varepsilon}$ such that $\|\phi^{\prime}(u^{\ast})\|_{\varepsilon}\leq\rho^{\ast}+\varepsilon$, cf. [AMR88, Lem. 4.3.7]. With this result at hand, we can exploit the norm equivalence together with the definition of Fréchet derivatives to conclude the desired local convergence rate. This simple argument is elaborated e.g. in [Shi81, Fin99] and generalizes the previous findings obtained in [Kit66]. ∎

In the finite dimensional case and under some additional assumptions on the structure of $\phi^{\prime}(u^{\ast})$ it can be shown that $\|u^{n}-u^{\ast}\|_{V}\leq C\hskip 2.0pt|\rho^{\ast}|^{n}\|u^{0}-u^{\ast}\|_{V}$, cf. [OR70, Th. 10.1.3 and 10.1.4; NR 10.1-5].

In the spirit of Proposition 3.1 we need to study the spectrum of $\phi^{\prime}(u^{\ast})$ to conclude local convergence.

3.1. Derivative of $\phi$

Since we interpret $\psi$ and $\phi$ as operators from $V$ to $V$, we have likewise for the first Fréchet derivative in $v\in V\setminus\{0\}$, $\psi^{\prime}(v\,;\,\cdot\,),\ \phi^{\prime}(v\,;\,\cdot\,)\colon V\to V.$ In order to apply Proposition 3.1, we need to compute the Fréchet derivative of $\phi$ in $u^{\ast}$. As a first step, we compute the derivative of the mapping $v\mapsto\|\psi(v)\|$, leading to

[TABLE]

or, more compactly, $\partial_{v}\|\psi(v)\|=(\psi^{\prime}(v),\phi(v))_{H}$. In order to compute $\phi^{\prime}(v)$, we write $\psi^{\prime}(v)$ in the form

[TABLE]

and conclude that

[TABLE]

For $v=u^{\ast}$ we can exploit $\phi(u^{\ast})=u^{\ast}$, $\|u^{\ast}\|=1$, and $u^{*}=(\lambda^{*}+\sigma)\,\psi(u^{\ast})$. This implies $\|\psi(u^{\ast})\|=|\lambda^{\ast}+\sigma|^{-1}$ and thus,

[TABLE]

To compute $\psi^{\prime}(u^{\ast})$, in turn, we use the identity

[TABLE]

Note that $J^{\prime}(v)\colon V\times V\rightarrow V^{*}$ denotes here the Fréchet derivative of $J$ in a fixed element $v\in V\setminus\{0\}$. We conclude

[TABLE]

Next, we want to verify that $J^{\prime}(u^{\ast};\,\cdot\,)\psi(u^{\ast})=0$. For that, we consider the definition of $J^{\prime}$ which yields

[TABLE]

Here we also used that in a neighborhood of $u^{*}$ (which excludes the zero element due to $\|u^{\ast}\|=1$) the operator $J$ is continuous. In summary, we proved $J^{\prime}(u^{\ast};w)\,u^{\ast}=0$ for all $w\in V$. Since $u^{*}$ equals $\psi(u^{\ast})$ up to a multiplicative constant, we conclude that $J^{\prime}(u^{\ast};\cdot)\,\psi(u^{\ast})=0$ and hence, $\psi^{\prime}(u^{\ast})=J_{\sigma}(u^{\ast})^{-1}\mathcal{I}$. Plugging this into previous estimates, we directly obtain

[TABLE]

3.2. Local convergence results

Since the $J$-method is of the from $u^{n+1}=\phi(u^{n})$, we recall from Proposition 3.1 that the local convergence rate is strongly connected to the spectrum of $\phi^{\prime}(u^{\ast})$. Besides the modified expression in (10), the analysis of this spectrum also requires the following orthogonality result.

Lemma 3.2**.**

Let $(u,\lambda)\in V\times\mathbb{R}$ be an eigenpair of $J(u^{*})$. Further, let $(w,\mu)\in V\times\mathbb{R}$ be an adjoint-eigenpair of $J(u^{*})$ with $\mu\neq\lambda$, i.e.,

[TABLE]

for all $v\in V$. Then, it holds that $(u,w)_{H}=0$.

Proof.

By definition we know that $\langle J(u^{*})\,u,v\rangle=\lambda\,(u,v)_{H}$ for all $v\in V$. Thus, with the test function $v=w$ we have

[TABLE]

The assumption ${\mu}\neq{\lambda}$ then directly implies the $H$-orthogonality of primal and adjoint eigenfunctions. ∎

With this, we are prepared to study the spectrum of the linear operator $\phi^{\prime}(u^{*})\colon V\to V$ to make conclusions on the convergence rate using Proposition 3.1. The corresponding eigenvalue problem reads: find $z\in V\setminus\{0\}$ and $\mu\in\mathbb{R}$ so that

[TABLE]

for all $v\in V$. The subsequent lemma relates the spectra of $\phi^{\prime}(u^{*})$ and  $J(u^{\ast})$. Recall that $J(u^{\ast})$ is a linear, bounded operator over $\mathbb{R}$. Thus, its (real) spectrum coincides with the (real) spectrum of the corresponding adjoint operator [HR11].

Lemma 3.3**.**

Assume that the (real) eigenvalues of $J(u^{\ast})$ are countable and that there exists a basis of corresponding adjoint-eigenfunctions of $V$. Let the eigenvalue of interest $\lambda^{*}$ be simple and $\sigma\in\mathbb{R}$ a shift such that $J_{\sigma}(u^{\ast})$ is invertible (note that in particular we have $\sigma\neq-\lambda^{*}$). Then, the spectra of $J(u^{\ast})$ and $\phi^{\prime}(u^{\ast})$ are connected in the following sense: $\mu\neq\lambda^{*}$ is an eigenvalue of $J(u^{\ast})$ if and only if $|\lambda^{*}+\sigma|/(\mu+\sigma)$ is an eigenvalue of $\phi^{\prime}(u^{\ast})$ and the eigenvalue $\lambda^{*}$ of $J(u^{\ast})$ corresponds to the eigenvalue [math] of $\phi^{\prime}(u^{\ast})$.

Proof.

Due to the assumption on the shift, $J_{\sigma}(u^{\ast})$ is invertible and $J_{\sigma}(u^{\ast})^{-1}\mathcal{I}$ has the eigenvalues $(\lambda_{1}+\sigma)^{-1},(\lambda_{2}+\sigma)^{-1},\dots$, including $(\lambda^{*}+\sigma)^{-1}$ with eigenfunction $u^{*}$. Note that the eigenvalues are shifted but the primal- and adjoint-eigenfunctions are the same as for $J(u^{\ast})$. We first consider the eigenvalue $(\lambda^{*}+\sigma)^{-1}$ of $J_{\sigma}(u^{\ast})^{-1}\mathcal{I}$ with eigenfunction $u^{*}$. Here we note that

[TABLE]

i.e., [math] is the corresponding eigenvalue of $\phi^{\prime}(u^{*})$. Now let $(\mu+\sigma)^{-1}$ be an eigenvalue of $J_{\sigma}(u^{\ast})^{-1}\mathcal{I}$ with $\mu\neq\lambda^{*}$. Since $J_{\sigma}(u^{\ast})^{-1}\mathcal{I}$ is a linear, bounded operator over $\mathbb{R}$, its (real) spectrum coincides with the (real) spectrum of the corresponding adjoint operator (cf. [HR11]) such that there exists a corresponding adjoint-eigenfunction $w$. By definition this means that

[TABLE]

for all $v\in V$. Using that $w$ is also an adjoint-eigenfunction of $J(u^{\ast})$ and the orthogonality of Lemma 3.2, we find that

[TABLE]

Thus, $w$ is an adjoint-eigenfunction of $\phi^{\prime}(u^{\ast})$ to the eigenvalue $|\lambda^{*}+\sigma|/(\mu+\sigma)$.

For the reverse direction let $(z,\mu)$ be an eigenpair of $\phi^{\prime}(u^{*})$, i.e., for all $v\in V$ we have

[TABLE]

For $\mu=0$ this directly leads to $z=u^{*}$, since $\mu=0$ yields the relation $J_{\sigma}(u^{*})^{-1}\mathcal{I}z=(J_{\sigma}(u^{*})^{-1}\mathcal{I}z,u^{*})_{H}\,u^{*}$. An application of $J_{\sigma}(u^{*})$ then implies that $u^{*}$ and $z$ coincide up to a multiplicative constant. Since both functions are normalized, we get $u^{*}=z$. In the case $\mu\neq 0$ we employ the adjoint-eigenfunctions $v_{1},v_{2},\dots$ of $J(u^{\ast})$ as test functions. This yields

[TABLE]

Note that the second term on the right-hand side vanishes due to Lemma 3.2. The definition of being an adjoint-eigenfunction then leads to

[TABLE]

for $j=1,2,\dots$ . Obviously, $(z,v_{j})_{H}$ cannot vanish for all test functions. This implies that $\mu=|\lambda^{*}+\sigma|/(\lambda_{j}+\sigma)$ for a certain index $j$. In other words, $\mu$ equals an eigenvalue of $J_{\sigma}(u^{*})^{-1}\mathcal{I}$ up to the factor $|\lambda^{*}+\sigma|$. ∎

With the obtained knowledge about the spectrum of $\phi^{\prime}(u^{\ast})$ we can directly apply Proposition 3.1 to deduce an abstract local convergence result with a rate depending on the spectrum of $J(u^{*})$ relative to the implemented shift as in the matrix case.

Theorem 3.4** (abstract local convergence).**

Consider the assumptions of Lemma 3.3 and let $\sigma\in\mathbb{R}$ be a shift sufficiently close to $-\lambda^{*}$. By $\mu$ we denote the (real) eigenvalue of $J(u^{*})$ which is closest to $-\sigma$ but different from $\lambda^{*}$. If the shift is selected such that

[TABLE]

then the iterations of the $J$-method (6) are locally convergent to $u^{\ast}$ in the $V$-norm. Furthermore, for every $\varepsilon>0$ there exists a neighborhood $S_{\varepsilon}$ of $u^{\ast}$ and a constant $C_{\varepsilon}>0$ such that

[TABLE]

and $n\geq 1$.

Below, we apply this abstract result to the GPEVP of Example 2.1.

4. Quantified local convergence of the shifted $J$-method for the GPEVP

As already mentioned in the introduction, we want to apply the damped $J$-method to the GPEVP (1), cf. Example 2.1. Thus, the task is to find an eigenfunction $u^{*}\in H^{1}_{0}(\mathcal{D})$ with $\|u^{*}\|^{2}_{L^{2}(\mathcal{D})}=1$ and corresponding eigenvalue $\lambda^{\ast}\in\mathbb{R}$ such that

[TABLE]

Recall that the potential satisfies $W\in L^{\infty}(\mathcal{D},\mathbb{R})$, whereas $\kappa\geq 0$ regulates the nonlinearity. Furthermore, $\mathcal{L}_{z}$ is the angular momentum operator with angular velocity $\Omega\in\mathbb{R}$. With these properties it is easily seen that the GPEVP can only have real eigenvalues. In the following we will also assume that

[TABLE]

This condition can be interpreted as that trapping frequencies are larger than the angular frequency. Physically speaking, this ensures that centrifugal forces do not become too strong compared to the strength of the trapping potential $W$.

4.1. The Gross-Pitaevskii energy

We consider homogeneous Dirichlet boundary conditions and the short-hand notation $L^{2}(\mathcal{D}):=L^{2}(\mathcal{D},\mathbb{C})$, $H^{1}_{0}(\mathcal{D}):=H^{1}_{0}(\mathcal{D},\mathbb{C})$, and $H^{-1}(\mathcal{D}):=H^{-1}(\mathcal{D},\mathbb{C})$ for the function spaces over $\mathbb{R}$, cf. the details in Example 2.1. We set

[TABLE]

This means that $\|\cdot\|$ equals the $L^{2}(\mathcal{D})$-norm and $(\,\cdot\,,\,\cdot\,)_{H}:=\Re(\,\cdot\,,\,\cdot\,)_{L^{2}(\mathcal{D})}$, where $(u,v)_{L^{2}(\mathcal{D})}:=\int_{\mathcal{D}}u\,\bar{v}\hskip 2.0ptdx$. Similarly, we equip $V$ with the inner product $(\,\cdot\,,\,\cdot\,)_{V}:=\Re(\,\cdot\,,\,\cdot\,)_{H^{1}(\mathcal{D})}$, where $(u,v)_{H^{1}(\mathcal{D})}:=\int_{\mathcal{D}}\nabla u\cdot\overline{\nabla v}\hskip 2.0ptdx$. As before, the weak formulation is given by $\mathcal{A}(u^{*})=A(u^{*},u^{*})=\lambda^{*}\mathcal{I}u^{*}$, where the nonlinear operator $A$ is defined in (3). The eigenvalue problem is equivalent to finding the critical points (on the manifold associated with the $L^{2}$-normalization constraint) of the energy functional $E\colon V\rightarrow\mathbb{R}$ given by

[TABLE]

cf. [DK10]. Recalling that $W_{\mbox{\rm\tiny R}}(x)=W(x)-\tfrac{1}{4}\,\Omega^{2}\,|x|^{2}$ and considering assumption (11), we see that the energy functional is bounded from below by a positive constant $c=c(\Omega,\mathcal{D})$, i.e.,

[TABLE]

Hence, $E$ is weakly lower semi-continuous and bounded from below, which yields the existence of a minimizer that is typically called a ground state.

If $\Omega=0$, i.e., in the absence of a rotating potential, then the GPEVP has infinitely many eigenvalues $0<\lambda_{1}^{\ast}<\lambda_{2}^{\ast}\leq\lambda_{3}^{\ast}\leq\cdots<\infty$, where the ground state eigenvalue $\lambda_{1}^{\ast}$ is simple, cf. [CCM10, HP20]. If $\Omega\not=0$, the smallest eigenvalue is typically no longer simple [BWM05]. This case refers to the physical phenomenon of a broken symmetry, where the ground state can have different shapes which differ in their number and location of vortices.

4.2. The $J$-operator for the GPEVP

In order to formulate the $J$-method for the GPEVP, we need to compute the linearization $J(u)$ according to (4). Hence, by calculating the Fréchet derivative of $\mathcal{A}(u)=A(u,u)$ using (3), we obtain

[TABLE]

where

[TABLE]

Note that the operator $J(u)$ induces an $\mathbb{R}$-bilinearform $\langle J(u)\,\cdot,\,\cdot\,\rangle$, i.e., we still consider $V$ as an $\mathbb{R}$-vector space and have bilinearity only for multiplicative constants in $\mathbb{R}$. Recall that we consider the real-Fréchet derivative, which also fits in the standard framework used in quantum mechanics. Moreover, the complex-Fréchet derivative does not exist for the present example.

We will show coercivity of $J(u)$ up to a shift in Lemma 5.1 below. Before that, it is worth to mention that the eigenvalue problem can be equivalently expressed in terms of $\mathcal{J}(u)$, which also yields the typical structure with standard inner products. We have the following result.

Proposition 4.1**.**

Consider the GPEVP with full operator $\mathcal{J}$ given by (14) and its real part $J$ given by (13). Then $\lambda^{\ast}\in\mathbb{R}$ is an eigenvalue with $L^{2}$-normalized eigenfunction $u^{\ast}\in V$, i.e.,

[TABLE]

if and only if

[TABLE]

Proof.

If $(\lambda^{\ast},u^{\ast})\in\mathbb{R}\times V$ solves $\langle\mathcal{J}(u^{\ast}),v\rangle=\lambda^{\ast}(u^{\ast},v)_{L^{2}(\mathcal{D})}$, then we can take the real part on both sides and obtain $\langle J(u^{\ast}),v\rangle=\lambda^{\ast}\Re(u^{\ast},v)_{L^{2}(\mathcal{D})}$. Vice versa, if $(\lambda^{\ast},u^{\ast})\in\mathbb{R}\times V$ solves $\Re\langle\mathcal{J}(u^{\ast}),v\rangle=\lambda^{\ast}\,\Re(u^{\ast},v)_{L^{2}(\mathcal{D})}$, then we can use test functions of the form $\mathrm{i}v\in V$ to obtain $\Im\langle\mathcal{J}(u^{\ast}),v\rangle=\lambda^{\ast}\,\Im(u^{\ast},v)_{L^{2}(\mathcal{D})}$. Multiplying the second equation with the complex number $\mathrm{i}$ and adding it to the first equation readily yields $\langle\mathcal{J}(u^{\ast}),v\rangle=\lambda^{\ast}(u^{\ast},v)_{L^{2}(\mathcal{D})}$. ∎

The proposition shows that we can either interpret the eigenvalue problem with real parts only (i.e., with $J$) or with the full operator $\mathcal{J}$. With analogous arguments, we can also prove that $J(u^{\ast})^{-1}\mathcal{I}=\mathcal{J}(u^{\ast})^{-1}\mathcal{I}^{\mathbb{C}}$, where

[TABLE]

This justifies that we can interpret the iterations of the damped $J$-method (7) equivalently with $J$ or $\mathcal{J}$. In particular, we have the following conclusion.

Conclusion 4.2**.**

Consider the GPEVP with full operator $\mathcal{J}$ given by (14) and its real part $J$ defined in (13). Given $u^{n}\in V$ with $\|u^{n}\|=1$ and a step size $\tau_{n}>0$, we assume that $J_{\sigma}(u^{n})^{-1}$ is invertible. Then the iterations of the $J$-method (7) can be equivalently characterized by the $\mathcal{J}$-iteration

[TABLE]

with $\mathcal{J}_{\sigma}(u):=\mathcal{J}(u)+\sigma\mathcal{I}^{\mathbb{C}}$ and $(\gamma_{n}^{\mathbb{C}})^{-1}:=(\mathcal{J}_{\sigma}(u^{n})^{-1}\mathcal{I}^{\mathbb{C}}\hskip 1.0ptu^{n},u^{n})_{L^{2}(\mathcal{D})}$. The assumed existence of $J_{\sigma}(u^{n})^{-1}\mathcal{I}$ implies the existence (and uniqueness) of $\mathcal{J}_{\sigma}(u^{n})^{-1}\mathcal{I}^{\mathbb{C}}$.

4.3. Local convergence

Finally, we apply Theorem 3.4 to the GPEVP.

Theorem 4.3** (quantified convergence for the GPEVP).**

Consider the GPEVP as described in Example 2.1 and let $u^{n}$ denote the iterates generated by the shifted $J$-method (without damping), i.e.,

[TABLE]

By $u^{\ast}\in V=H^{1}_{0}(\mathcal{D})$ we denote an $L^{2}$-normalized eigenfunction to (1) with eigenvalue $\lambda^{\ast}$. Assume that $\lambda^{\ast}$ is a simple eigenvalue of $J(u^{\ast})$ and that the shift $\sigma\not=-\lambda^{\ast}$ is such that $J_{\sigma}(u^{\ast})$ has a bounded inverse. Further, let the shift $\sigma$ be sufficiently close to $-\lambda^{\ast}$ and let $\mu$ be the eigenvalue of $J_{\sigma}(u^{\ast})$ so that $|\mu+\sigma|$ is minimal and

[TABLE]

Then there exists a neighborhood $S\subset V$ of $u^{\ast}$ such that for all starting values $u^{0}\in S$ we have

[TABLE]

Furthermore, for every $\varepsilon>0$ there is a neighborhood $S_{\varepsilon}\subset V$ of $u^{\ast}$ and a constant $C(\varepsilon)>0$ such that for all $u^{0}\in S_{\varepsilon}$ (and $n\geq 1$) it holds

[TABLE]

Hence, if the shift $\sigma$ is close to $-\lambda^{\ast}$ and $u^{n}$ close to $u^{*}$, then we have locally a linear convergence with rate

[TABLE]

for any $\varepsilon>0$.

Proof.

By definition of the problem (cf. Example 2.1) all eigenvalues are real. Furthermore, the spectrum is countable and does not have an accumulation point in $\mathbb{C}$. This is a direct consequence from the observation that $J_{\sigma}(u^{\ast})^{-1}\mathcal{I}\colon L^{2}(\mathcal{D})\rightarrow L^{2}(\mathcal{D})$ is a compact operator for all shifts $-\sigma$ that are not in the spectrum of $J(u^{\ast})$ (i.e., we have compact resolvents). Hence, we can apply Theorem 3.4 which readily proves the result. ∎

In the finite dimensional case based on a finite difference discretization, a corresponding result was presented in [JKM14]. Motivated by the above convergence result, a practical realization of the iterations can be based on the more natural formulation of the $J$-method given by (16). This is also the version for which we discuss the implementation in Appendix B.

5. Global convergence of the damped $J$-method for the GPEVP

In this section we come back to the question of invertibility of the operator $J$ (which hence also implies invertibility of $\mathcal{J}$). This will then lead to a globally convergent method.

5.1. Coercivity of the shifted $J$-operator

We first show that the operator $J$ is – up to a shift – coercive.

Lemma 5.1**.**

Given $u\in V$ and assumption (11), the operator $J(u)$ corresponding to the Gross-Pitaevskii operator satisfies a Gårding inequality. More precisely, for any $\sigma\geq\frac{\kappa}{3}\|u\|_{L^{4}(\mathcal{D})}^{4}/\|u\|^{4}$ the bilinear form

[TABLE]

is coercive and thus, the operator $J_{\sigma}(u)=J(u)+\sigma\mathcal{I}\colon V\to V^{*}$ is invertible.

Proof.

Consider $v\in V$. We start with considering the rotational gradient, for which we observe with Young’s inequality that

[TABLE]

Hence, with $W_{\mbox{\rm\tiny R}}=W(x)-\tfrac{1}{4}\,\Omega^{2}|x|^{2}$ and assumption (11) we have

[TABLE]

This leads to

[TABLE]

To estimate the negative part, we apply once more Young’s inequality with some parameter $\mu>0$ and the Cauchy Schwarz inequality to get

[TABLE]

Using this estimate in (5.1) with the particular choice $\mu=3\,\|u\|^{2}/\|u\|^{4}_{L^{4}(\mathcal{D})}$, we obtain

[TABLE]

Thus, we conclude

[TABLE]

Recall the definition of the energy $E$ in Section 4.1. Motivated by the previous lemma, we also define the shifted energy $E_{\sigma}(u):=E(u)+\frac{1}{2}\sigma\|u\|^{2}$ such that $E_{0}(u)=E(u)$. For $u\in V$ with normalization constraint $\|u\|=1$, a sufficient shift in the sense of Lemma 5.1 is thus given by

[TABLE]

Note that for a normalized function $u$ we can express the Rayleigh quotient in terms of the energy by

[TABLE]

In particular, this formula relates the eigenvalues with the energies of the eigenfunctions.

5.2. Feasibility of the $J$-method

For the feasibility of the damped $J$-method, we need to guarantee a priori that $J_{\sigma}(u^{n})$ stays invertible throughout the iteration process. We fix the shift $\sigma$ in the beginning of the iteration, e.g., by $\sigma:=\frac{4}{3}E(u^{0})$. Now, the aim is to show that the energy of the iterates does not increase such that $\sigma\geq\frac{4}{3}E(u^{n})$ for all $n\geq 0$.

Lemma 5.2**.**

Recall the shifted energy $E_{\sigma}(u):=E(u)+\frac{1}{2}\sigma\|u\|^{2}$ with $E(u)$ being defined in (12). Further recall the definition of the $J$-method (7) with iteration steps $\tilde{u}^{n+1}=(1-\tau)u^{n}+\tau\,\gamma_{n}\,J_{\sigma}(u^{n})^{-1}\mathcal{I}u^{n}$. Then, for any step size $\tau\leq\frac{1}{2}$ we have the following guaranteed estimate for the difference of the shifted energies

[TABLE]

Proof.

We start by establishing a couple of identities for different evaluations of $J$. Here we exploit the $L^{2}$-orthogonality (8), i.e., $(u^{n},\tilde{u}^{n+1}-u^{n})=0$. Together with $\|u^{n}\|=1$ this implies $(u^{n},\tilde{u}^{n+1})=1$. Using these facts, we observe that

[TABLE]

Next, using again the $L^{2}$-orthogonality together with the definition of $\tilde{u}^{n+1}$ in (7), we have

[TABLE]

With this equality, we conclude

[TABLE]

Here we need to have a closer look at the term $\langle J_{\sigma}(u^{n})(\tilde{u}^{n+1}-u^{n}),u^{n}\rangle$. By the definition of $J(u^{n})$ it is easily seen that

[TABLE]

Plugging this into (22) yields for the shifted energy

[TABLE]

Since $\Re\,z=\Re\,\overline{z}$ for all $z\in\mathbb{C}$, we conclude that

[TABLE]

Using that

[TABLE]

and $\Re(u^{n}(\overline{\tilde{u}^{n+1}-u^{n}}))=\Re(\tilde{u}^{n+1}\overline{u^{n}})-|u^{n}|^{2}$ we see that

[TABLE]

Finally, an application of the triangle and Young’s inequality yields the estimate

[TABLE]

which then implies the desired inequality if $\tau\leq\frac{1}{2}$. ∎

With this result we are now in the position to prove uniform boundedness of the energy of the iterates and even a reduction of the energy.

Theorem 5.3** (energy reduction).**

Given $u^{0}\in V$ with $\|u^{0}\|=1$, we let $\sigma\geq\frac{4}{3}E(u^{0})$ and $W(x)\geq\Omega^{2}|x|^{2}$. Then, there exists a $\tau^{\ast}>0$ (that only depends on $u^{0}$ and its energy) such that for all $0<\tau_{n}\leq\tau^{\ast}$ the sequence obtained by the damped $J$-method (7) is well-posed and strictly energy diminishing for all $n$, i.e., it holds

[TABLE]

where $E(u^{n+1})<E(u^{n})$ if $u^{n}$ is not already a critical point of $E$ and thus an eigenstate.

Proof.

We proceed inductively. From estimate (17) in the proof of Lemma 5.1 we know that $J_{\sigma}(u^{0})$ is coercive with constant $1/2$, i.e., $\|\nabla v\|^{2}=\|v\|_{V}^{2}\leq 2\,\langle J_{\sigma}(u^{0})v,v\rangle$. Together with (21) and the $L^{2}$-orthogonality (8), this implies

[TABLE]

Hence, we get

[TABLE]

Assume that $\tau_{0}\leq\tau^{\ast}\leq 2$, where $\tau^{\ast}$ is selected sufficiently small compared to $C_{0}$. Then we have that $\|\nabla\tilde{u}^{1}-\nabla u^{0}\|<1$ and the Sobolev embedding of $L^{4}(\mathcal{D})\hookrightarrow H^{1}_{0}(\mathcal{D})$ with constant $C_{S}$ implies that

[TABLE]

Consequently, we can use (19) and (17) to observe that the energy difference fulfills

[TABLE]

Hence, if

[TABLE]

then we have

[TABLE]

where we have used that the iterations increase the mass intermediately, i.e., $\|\tilde{u}^{1}\|\geq 1$ as shown in (9). Note that since $u^{0}$ and $u^{1}$ are normalized in $L^{2}(\mathcal{D})$, we can drop the shift, leading to $E(u^{1})\leq E(u^{0})$. Hence, we have $\sigma\geq\frac{4}{3}E(u_{0})\geq\frac{4}{3}E(u_{1})$ and Lemma 5.1 guarantees that $J_{\sigma}(u^{1})$ is still coercive. Inductively, we can repeat the arguments for $u^{n}$ with the same generic constant $C_{0}$ to show that for $\tau^{n}\leq\tau^{\ast}$ we have

[TABLE]

Since the energy is diminished in every iteration, the coercivity of $\langle J_{\sigma}(u^{n})\hskip 2.0pt\cdot\hskip 2.0pt,\hskip 2.0pt\cdot\hskip 2.0pt\rangle$ is maintained and all the iterations are well-defined. Finally, we note that because of (26) we have $E_{\sigma}(u^{n})=E_{\sigma}(\tilde{u}^{n+1})$ if and only if $u^{n}=\tilde{u}^{n+1}$. However, this can only happen if

[TABLE]

i.e., if $u^{n}$ is already an eigenfunction with eigenvalue $\lambda=\gamma_{n}$. ∎

It is interesting to note that the $L^{2}$-norm of $\tilde{u}^{n}$ cannot diverge. We see this in the following conclusion.

Conclusion 5.4**.**

In the setting of Theorem 5.3 it holds that $\|\tilde{u}^{n}\|\rightarrow 1$ for $n\rightarrow\infty$.

Proof.

In the proof of Theorem 5.3 we have seen that

[TABLE]

Since $E_{\sigma}(u^{n})$ is monotonically decreasing and bounded from below, we have $E_{\sigma}(u^{n})-E_{\sigma}(u^{n+1})\rightarrow 0$. This together with the Poincaré-Friedrichs inequality implies that

[TABLE]

for $n\rightarrow\infty$. ∎

5.3. Convergence and optimal damping

In this subsection we prove the convergence of the $J$-method for a suitable choice of damping parameters. We can make practical use of this result by selecting $\tau_{n}$ in each iteration step by the minimizer of a simple one-dimensional minimization problem.

Theorem 5.5** (global convergence).**

Suppose that the assumptions of Theorem 5.3 are fulfilled. Additionally assume that $\tau_{n}$ is selected such that it does not degenerate, i.e., is uniformly bounded away from zero. Then there exists a limit energy $E^{\ast}:=\lim_{n\rightarrow\infty}E(u^{n})$ and, up to a subsequence, we have that the iterates $u^{n}$ of the damped $J$-method converge strongly in $H^{1}(\mathcal{D})$ to a limit $u^{\ast}\in V$. The limit is an $L^{2}$-normalized eigenfunction with some eigenvalue $\lambda^{\ast}>0$, i.e.,

[TABLE]

and we have $E(u^{\ast})=E^{\ast}$. If $u^{\ast}$ is the only eigenfunction on the energy level $E^{\ast}$, then we have convergence of the full sequence $u^{n}$.

Proof.

The proof is similar to the arguments presented in [HP20, Th. 4.9]. First, Theorem 5.3 guarantees the existence of the limit $E^{\ast}:=\lim_{n\rightarrow\infty}E(u^{n})$. Hence, $u^{n}$ is uniformly bounded in $V$ and we can extract a subsequence (still denoted by $u^{n}$) that converges weakly in $H^{1}(\mathcal{D})$ and strongly in $L^{p}(\mathcal{D})$ (for $p<6$) to a limit $u^{\ast}\in V$ with $\|u^{\ast}\|=1$. Since $J_{\sigma}(u^{\ast})$ is a real-linear operator that depends continuously on the data and which induces the coercive bilinear form $\langle(J(u)+\sigma\mathcal{I})\hskip 3.0pt\cdot\hskip 2.0pt,\hskip 2.0pt\cdot\hskip 2.0pt\rangle$, we have that

[TABLE]

Together with the strong convergence $u^{n}\rightarrow u^{\ast}$ in $L^{4}(\mathcal{D})$, we conclude that

[TABLE]

This shows that

[TABLE]

Furthermore, we have seen in the proof of Theorem 5.3 (respectively Conclusion 5.4) that the strong energy reduction implies that for $n\rightarrow 0$

[TABLE]

and consequently we have with $\tilde{u}^{n+1}=(1-\tau_{n})u^{n}+\tau_{n}\gamma^{n}J_{\sigma}(u^{n})^{-1}\mathcal{I}u^{n}$ and the boundedness of $\tau_{n}$ that

[TABLE]

strongly in $H^{1}(\mathcal{D})$. Since we already know that $u^{n}$ converges weakly in $H^{1}(\mathcal{D})$ to $u^{\ast}$, we can now conclude that this is even a strong convergence and we have

[TABLE]

This shows that $u^{\ast}$ is an eigenfunction with eigenvalue $\gamma^{\ast}$. The strong $H^{1}$-convergence also implies convergence of the energies, i.e., $E^{\ast}=\lim_{n\rightarrow\infty}E(u^{n})=E(u^{\ast})$. ∎

For all sufficiently small $\tau_{n}$, Theorem 5.3 proves the energy reduction and Theorem 5.5 global convergence. However, since we do not know a priori what a sufficiently small value for $\tau_{n}$ is, we can combine the damped $J$-method with a line search algorithm that optimizes $\tau_{n}$ in each iteration step such that the energy reduction is (quasi) optimal. Theorems 5.3 and 5.5 show that such an optimal $\tau_{n}$ exists and that it does not degenerate to zero. We stress that finding such a $\tau_{n}$ does not require any additional inversions, which makes the procedure very cheap, cf. Appendix A for details.

Conclusion 5.6** ($J$-method with optimal damping).**

Consider a shift $\sigma$ such that the assumptions of Theorem 5.3 are fulfilled. Given $u^{n}\in V$ with $\|u^{n}\|=1$ the next iteration is obtained by selecting the optimal damping parameter with

[TABLE]

and defining $u^{n+1}$ as in (7). The approximations are energy diminishing and converge (up to a subsequence) strongly in $V$ to an $L^{2}$-normalized eigenfunction of the GPEVP.

Finally, if there is no rotation, we can even achieve guaranteed global convergence to the ground state provided that the selected initial value is non-negative.

Proposition 5.7** (global convergence to ground state).**

Assume the setting of Theorem 5.5. Furthermore, we consider that there is no rotation, i.e., $\Omega=0$, and a non-negative starting value $u^{0}\in V$, i.e., $u^{0}\geq 0$. If $\tau_{n}\leq 1$ and if the shift parameter $\sigma>0$ is selected sufficiently large, then the iterates $u^{n}$ of the damped $J$-method converge strongly in $H^{1}(\mathcal{D})$ to the unique (positive) ground state $u^{\ast}\geq 0$.

Proof.

If $\Omega=0$, then the problem can be fully formulated over $\mathbb{R}$ and admits a unique positive ground state $u^{\ast}\in V$, cf. [CCM10]. The only other ground state is $-u^{\ast}$. Furthermore, all excited states (i.e., all other eigenfunctions) must necessarily change their sign on $\mathcal{D}$, cf. [HP20, Lem. 5.4]. Hence, if we can verify that the iterates of the damped $J$-method (7) preserve positivity, then Theorem 5.5 guarantees that the global $H^{1}$-limit must be the desired ground state $u^{\ast}$.

Recall the damped $J$-iteration from (7), which (over $\mathbb{R}$) reduces to

[TABLE]

where $S_{\sigma}$ is the linear self-adjoint operator given by

[TABLE]

and $G$ characterizes the non-symmetric rank-$1$ remainder, i.e.,

[TABLE]

Analogously to the Sherman–Morrison formula for matrices [SM50], we can see that

[TABLE]

Consequently we can write the effect of the inverse as

[TABLE]

Since $S_{\sigma}$ is a self-adjoint elliptic operator, it preserves positivity, i.e., we have $S_{\sigma}^{-1}\mathcal{I}v\geq 0$ if $v\geq 0$. This immediately follows by writing the action of $S_{\sigma}^{-1}$ as an energy minimization problem. Starting (inductively) from a function $u^{n}\geq 0$ we conclude that

[TABLE]

If we can ensure that $1-(u^{n},S_{\sigma}^{-1}\mathcal{I}f)_{L^{2}(\mathcal{D})}>0$, then we have $J_{\sigma}(u^{n})^{-1}\mathcal{I}u^{n}\geq 0$. Let us hence consider $(u^{n},S_{\sigma}^{-1}\mathcal{I}f)_{L^{2}(\mathcal{D})}$, for which we obtain

[TABLE]

Since $S_{\sigma}$ is self-adjoint, we have

[TABLE]

where $\lambda_{\text{min}}(S_{0})>0$ is the minimal eigenvalue of $S_{0}$. Consequently, if the shift is such that

[TABLE]

then we have positivity of $J_{\sigma}(u^{n})^{-1}\mathcal{I}u^{n}$. Note that by the energy reduction property, we can bound $\|u^{n}\|^{3}_{L^{6}(\mathcal{D})}$ uniformly for all $n$ by a constant that only depends on the initial energy $E(u^{0})$. Together with the obvious positivity of $(1-\tau_{n})u^{n}$ for $\tau_{n}\leq 1$, we conclude the existence of a sufficiently large shift $\sigma$ so that $u^{n+1}\geq 0$ for all $n\geq 0$ and hence global convergence to the ground state. ∎

6. Numerical experiments

This section concerns the numerical performance of the proposed $J$-method enhanced by shifting and/or damping as outlined in Sections 2.3 and 2.4. As a general model, we seek critical points of the Gross-Pitaevskii energy (12) in a bounded domain $\mathcal{D}=(-L,L)^{2}$ for some parameter $L>0$. The particular choice of $L$ as well as the other physical parameters $\Omega$, $W$, and $\kappa$ will be specified separately in the various model problems below.

For the spatial discretization we always use bilinear finite elements on a Cartesian mesh of width $h=2^{-8}L$. We will not investigate discretization errors with respect to the underlying PDE. For approximation properties of discrete eigenfunctions we refer to the analytical results presented in [CCHM18, CCM10, CGZ10, HMP14] and to [CDM*+*14, HSW19, XX16] for a posteriori estimators and adaptivity. Our focus is the performance of the iterative eigenvalue solver promoted in this paper. As a measure of accuracy we will use the $L^{2}(\mathcal{D})$-norm of the residual $\mathcal{A}(u^{n})u^{n}-\lambda^{n}\mathcal{I}u^{n}$ given an approximate finite element eigenpair $(\lambda^{n},u^{n})$. We will stop the solver whenever the residual falls below the tolerance $\mbox{\rm\small TOL}=10^{-8}$.

For a better assessment of the performance of the $J$-method, we compare it with the projected $a_{z}$-Sobolev gradient flow introduced in [HP20]: given $u^{0}\in V$ with $\|u^{0}\|=1$, define for $n=1,2,\ldots$,

[TABLE]

We will refer to this approach as the $\mathcal{A}$-method (without shift). Note that every step of the $\mathcal{A}$-method can be interpreted as an energy minimization problem such that the iteration is positivity preserving for the case without rotation. Thus guarantees global convergence to the ground state for every nonnegative starting value. The corresponding shifted version, which we will also consider in the experiments, considers $\mathcal{A}_{\sigma}(u^{n}):=\mathcal{A}(u^{n})+\sigma\mathcal{I}$ in place of $\mathcal{A}(u^{n})$ in the definition of $\hat{u}^{n+1}$ and a corresponding adjustment of the normalization factor $\gamma_{n}$. According to the numerical experiments of [HP20], this method is representative for the larger class of gradient flows in terms of accuracy-cost ratios. The cost per iteration step for both $\mathcal{A}$- and $J$-method are proportional and of the same order. Tentatively, the $\mathcal{A}$-method is cheaper by a fixed factor, since (due to the rank-$1$ matrix that appears in the $J$-version) an additional linear system per step has to be solved when the Sherman-Morrison formula is used, cf. Appendix B. To what extent the computational overhead of the $J$-method can be reduced by suitable preconditioned iterative solvers is beyond the scope of the paper. Notwithstanding the above, the experiments below clearly show that the $J$-method easily compensates its possible computational overhead per step by a notedly smaller iteration count.

6.1. Ground state in a harmonic potential

In the first model problem, we consider a harmonic trapping potential with trapping frequencies $1/2$, i.e.,

[TABLE]

The angular momentum $\Omega$ is set to zero and the repulsion parameter to $\kappa=1000$. The size of the domain is chosen as $L=8$. This is larger than the Thomas-Fermi radius of the problem which can be estimated as $R^{\mbox{\scriptsize TF}}=\sqrt{2}\hskip 2.0pt(\kappa/\pi)^{1/4}\approx 5.97$, cf. [Bao14]. We are interested in computing the ground state, i.e., the global minimizer $u^{\text{gs}}$ of the energy (12). Note that, up to sign, $u^{\text{gs}}$ is the unique eigenfunction that corresponds to the smallest eigenvalue $\lambda^{\text{gs}}$ which is well-separated from the remaining spectrum. As an initial value for all variants of eigenvalue solvers we use the bi-quadratic bubble

[TABLE]

interpolated in the finite element space and normalized in $L^{2}(\mathcal{D})$. For this simple model problem, there are certainly more sophisticated initial guesses such as the Thomas-Fermi approximation. Our uneducated initial guess marks an additional challenge. As ground state energy we computed $E^{\text{gs}}:=E(u^{\text{gs}})\leq 6.019$. The corresponding eigenvalue approximation is $\lambda^{\text{gs}}\approx 17.93$. Clearly, the accuracy of these numbers is limited by the choice of the discretization parameter $h$. Mesh adaptivity as used in [HSW19] or a higher-order method would certainly help to improve on these numbers.

Figure 1 shows the evolution of the residuals during the iteration of several variants of the $J$ and $\mathcal{A}$-method. Unless specified differently, the general $J$- and $\mathcal{A}$-methods refer to a combination of damping and shifting. More precisely, we use damping for globalization of convergence until the residual falls below $10^{-3}$. Then we freeze the time step $\tau=1$ and switch to a Rayleigh shift strategy to possibly accelerate convergence, i.e., in each step we choose the shift $\sigma=-\lambda^{n}$ to be the current eigenvalue approximation. According to our numerical experience the coexistence of damping and shifting is hard to control. The transition from damping to shifting is clearly seen in the convergence plot of Figure 1. We observe linear but fairly slow convergence in the damped phase. As soon as we switch to shifting, the convergence is beyond linear. The $\mathcal{A}$-method performs similarly in the damping phase but diverges as soon as the shift is turned on. It is this phenomenon already observed in [JKM14] in a less extreme characteristic (see also the second model problem below) that motivated the derivation of a shift-sensitive $J$-method. Note that the improved rate can be explained by its connection to Newton’s method, cf. [JKM14, Sect. 3.4]. Convergence proofs with higher rate, however, are only given in the discrete setting for linear and particular nonlinear eigenvalue problems [Osb64, MV04].

In Figure 1 we also show results for the variants of the $J$- and $\mathcal{A}$-methods where either $\tau$, $\sigma$, or both are fixed. Their performance is in between the aforementioned combined approaches.

6.2. Exponentially localized ground state in a disorder potential

In the second model problem the non-negative external potential $W$ reflects a high degree of disorder and the repulsion parameter $\kappa$ is small. In this situation, the low-energy eigenstates essentially localize in the sense of an exponential decay of their moduli.

The numerical approximation of localized Schrödinger eigenstates in the linear case, i.e., for $\kappa=0$, has recently caused a large interest in the fields of computational physics and scientific computing [FM12, ADJ*+*16, Ste17, ADF*+*19, XZO19]. In particular, the results of [AHP20a] provide a mathematical justification of the observed localization. In the nonlinear case the phenomenon is still observable but locality deteriorates with increasing interaction [APV18, AP19, AHP20b]. Here, we choose $\kappa=1$ which leads to a fairly localized ground state as it can be seen in Figure 2 (right). Its computation turns out to be much more challenging than in the case of a harmonic trapping potential in the sense that convergence rates are slower and iteration counts larger. This is related to a possible clustering of the lowermost eigenvalues. In particular, we expect $\lambda_{1}/\lambda_{2}\approx 1$ in this example such that shifting can provide a considerable speed up. Due to the small repulsion parameter $\kappa$, however, we expect a significant gap within the first few eigenvalues as in the linear case.

We have tested the same solvers as in the previous subsection. The $J$-method involving shifting performs best by far. Surprisingly, the variant without adaptive time step in the damping phase even performed better. This is no contradiction with the theory as we are showing residuals rather than energies. Moreover, a locally optimal energy decrease does not necessarily lead to better global performance. Still the difference is not too big. As a general recommendation from our numerical experience we would favor to use an adaptive time step because it was more robust. Another difference with regard to the harmonic potential is that, this time, the $\mathcal{A}$-method reacts upon shifting in a positive way. For a few steps the convergence is indeed accelerated. However, thereafter the method turns back to a linear regime of convergence which cannot compete with the shifted $J$-method. Nevertheless, the $\mathcal{A}$-method does not fail completely as seen in the previous example. This may be caused by the small value of $\kappa$, compared with the considered potential.

6.3. Vortex lattices in a fast rotating trap

We close the numerical illustration of the $J$-method with a qualitative study of vortex lattice states in the present of fast rotating potentials. We choose a harmonic potential $W$ as in (28) and set $\Omega=0.99$, $\kappa=1000$, and the size of the computational domain to $L=10$.

We have computed four different eigenfunctions using the $J$-method. This was only possible using the shift-sensitivity of the $J$-method. We shall briefly describe the computational parameters. We use the bi-quadratic bubble (29) as the initial value for all computations. To compute the (tentative) ground state $u_{1}$ (see Figure 4, upper left) we used the combined strategy as before. However, we switched from damping to shifting only once the residual falls below $10^{-6}$. Switching earlier led to states of higher energy. E.g., switching at a tolerance of $10^{-3}$ lead to the eigenfunction $u_{2}$ depicted in the upper right of Figure 4. It is interesting to observe that while $E_{1}:=E(u_{1})<E(u_{2}):=E_{2}$ the corresponding eigenvalues are ordered the other way around.

Two further excited states are found by limiting the adaptive shift to the interval $[15.0,15.6]$ for $u_{3}$ (see lower left of Figure 4) and to the interval $[15.2,15.45]$ for the state $u_{4}$ that does not show any vortices (see lower right of Figure 4). In both cases we used the lower end of the interval as the shift in the damping phase and we switched to adaptive shifting at residual tolerance $10^{-4}$ for $u_{3}$ and $10^{-3}$ for $u_{4}$. Note that $u_{1}$ seems to be the global energy minimizer of the (discretized) problem but the exited states $u_{2},u_{3},u_{4}$ do not necessarily represent the next higher energy levels $2$ to $4$ but some levels of higher energy.

From this rather complicated derivation one can see that it is by no means trivial to compute these excited states. We shall also say that it is not always easy to control the shifting. If one shifts too early in the sense that the approximation is not yet close to the target eigenfunction (e.g. in terms of number of vortices) the procedure may fail completely. Despite this difficulty which is intrinsic to the nonlinear eigenvalue problem at hand, the $J$-method along with shifting and damping enables the selective approximation of excited states as well as the amplification of convergence beyond linear rates in the spirit of the Rayleigh quotient iteration even in this challenging regime of vortex pattern formation.

7. Conclusion

In this paper we have generalized the $J$-method proposed in [JKM14] to the abstract Hilbert space setting. This gives rise to a variational formulation that is straightforwardly accessible by Galerkin-type discretizations, e.g., based on finite elements. Moreover, we have transferred the proof of local convergence of the $J$-method from the discrete setting (cf. [JKM14]) to the abstract setting and recovered a quantitative convergence rate that depends on the spectral shift. Since this fast convergence is indeed a local feature, we have proposed a damped $J$-method. For the GPEVP, the damping step can be seen as a discretization of a generalized gradient flow and guarantees reduction of the energy associated to the Gross-Pitaevskii operator. This energy reduction is the key to the global convergence of the damped method.

We have proposed a combined strategy of damping and shifting, depending on the residual error. The damping part guides the iterates to a sufficiently small neighborhood of an eigenfunction. Therein, the shifting significantly improves the linear rate of convergence. With a Rayleigh-type shifting strategy remarkable speed-ups beyond linear convergence are observed. In numerical experiments we have demonstrated the excellent performance of the arising method and its suitability for both the computation of ground states and the selective computation of excited states. We believe that the proposed strategy can be also an efficient tool for treating other types of eigenvalue problems with nonlinearities, in particular those that can be rephrased as finding the critical points of constraint energy minimization problems.

Appendix A Energy-diminishing step size control

We consider the damped $J$-method (7) in the case of the Gross-Pitaevskii equation. In order to implement an efficient step size control, we consider the function

[TABLE]

that we want to minimize for $\tau\in(0,2)$. Based on $u^{n}$, we compute

[TABLE]

Note that this implies

[TABLE]

With this, we precompute various terms, which are given by

[TABLE]

as well as

[TABLE]

[TABLE]

Note that the terms $\alpha_{i}$, $\beta_{i}$, and $\zeta_{i}$ have to be computed only once per time step (e.g. within one loop over the grid elements). Finally, we define the function

[TABLE]

and observe that $f(\tau)$ is given by

[TABLE]

After precomputing $\alpha_{i}$, $\beta_{i}$, and $\zeta_{i}$, the function $f(\tau)$ can be cheaply evaluated. The minimization step, i.e., $\tau_{n}:=\mbox{arg min}\hskip 2.0pt\{f(\tau)|\hskip 3.0pt\tau\in(0,2)\}$ can be easily implemented using, e.g., a golden section search. Note that the energy of $u^{n+1}$ is now given by $f(\tau_{n})$.

Appendix B Matrix representation of the $\mathcal{J}$-operator

The $\mathcal{J}$-operator in case of the GPEVP was derived in Section 4.2 and we consider the iteration given by (16). For the implementation we need to discuss the handling of the nonlinear terms. Using the identity $(u\overline{v}+2v\overline{u})u\,\overline{w}=2\Re\left(u\overline{v}\right)u\,\overline{w}+|u|^{2}v\,\overline{w}$, we end up with the integrals

[TABLE]

Note that we only need to consider real test functions $w$ and decompose $u,v$ into its real and imaginary part, i.e., $u=u_{R}+\mathrm{i}u_{I}$, $v=v_{R}+\mathrm{i}v_{I}$. This is also done in the finite element discretization, i.e., we work with real vectors of double dimension. For the first integral we note that

[TABLE]

Introducing $M_{v}$ as the mass matrix weighted by $v$, this leads to the matrix representation

[TABLE]

For the second term we have |u|^{2}v\,\overline{w}=|u|^{2}\big{(}v_{R}w+\mathrm{i}v_{I}w\big{)}. Thus, the corresponding finite element matrix is block diagonal and reads $\text{diag}(M_{|u|^{2}},M_{|u|^{2}})$. Finally, we have

[TABLE]

A simple rearrangement shows that this corresponds to the rank-one matrix

[TABLE]

In total, a finite element discretization of the $J$-method calls for a solution of a linear system which is decomposed of several sparse matrices and the latter rank-$1$ update. This can be easily inverted using the Sherman-Morrison formula, cf. [SM50, Woo50].