Existence of Optical Vortex Solitons in Photorefractive Media

TL;DR

This paper proves the existence of optical vortex solitons in photorefractive media, providing theoretical constraints and numerical methods to analyze their properties in nonlinear optical physics.

Contribution

It introduces an existence theory for vortex solitons in photorefractive media, using constrained minimization and finite element methods.

Findings

Existence of positive exponentially decaying vortex solutions.

Constraints on propagation constants for vortex stability.

Numerical computation methods for vortex profiles.

Abstract

Optical propagation and vortices in nonlinear media have been intensively studied in modern optical physics. In this paper, we establish constraints regarding the propagation constant and provide an existence theory and numerical computations for positive exponentially decaying solutions for a class of ring-profiled solitons in a type of nonlinear media known as a photorefractive nonlinearity. Our methods include constrained minimization and finite element formalism, and we study the vortex profile and its propagation by fixing the energy flux.

Equations

$$\begin{array}{c}i{\varphi }_z+\frac{1}{2k_1}{\nabla }_{\perp }^2\varphi -\Delta n_1(I)\varphi =0,\\ i{\psi }_z+\frac{1}{2k_2}{\nabla }_{\perp }^2\psi -\Delta n_2(I)\psi =0.\end{array}$$

$$\Delta n_1(I)=\frac{P}{1+I}\text{and}\Delta n_2(I)=\frac{Q}{1+I},$$

$$\varphi (r)=e^{i({\beta }_{1}z+m_1\theta )}u(r),\psi (r)=e^{i({\beta }_{2}z+m_2\theta )}v(r),$$

$$\begin{array}{c}(r{u}_r{)}_r-\frac{m_1^2}{r}u-\frac{ru}{1+(u^2+v^2)}={\beta }_{1}ru,0<r<R,\\ (r{v}_r{)}_r-\frac{m_2^2}{r}v-\frac{rv}{1+(u^2+v^2)}={\beta }_{2}rv,0<r<R\end{array}$$

$$u(0)=u(R)=v(0)=v(R)=0.$$

$$(r{u}_r{)}_r-\frac{m^2}{r}u-\frac{ru}{1+\alpha u^2}$$

$$u(0)=u(R)$$

$$u(r)>0\text{ for all }r\in (0,R).$$

$$I(u)=\frac{1}{2}{\int }_0^R\left\{r{u}_r^2+\frac{m^2}{r}{u}^2+\frac{r}{\alpha }\ln (1+\alpha u^2)\right\}dr$$

$$P(u)=2\pi {\int }_0^{R}r{u}^{2}dr=P_0,$$

$$I_0=\underset{u\in \mathcal{U}}{\inf }\left\{I(u)|P(u)=P_0\right\},$$

$$\mathcal{U}=\{u(r)\text{ absolutely continuous over }[0,R],u(0)=u(R)=0,\mathcal{E}(u)<\infty \}$$

$$\mathcal{E}(u)={\int }_0^R\left\{r{u}_r^2+\frac{m^2}{r}{u}^2\right\}dr.$$

$$\beta$$

$$\underset{0<r<R}{\max }\{u(r){\}}^2>-\frac{1}{\alpha }\left[{\left(\beta +\frac{m^2}{R^2}\right)}^{-1}+1\right]$$

$$\beta >-\frac{m^2}{R^2}-1,$$

$$u(r{)}^2\le C_{{ϵ}_0}{e}^{-{\sqrt{ϵ}}_{0}r}$$

$${ru{u}_r|}_0^R-{\int }_0^{R}r{u}_r^{2}dr$$

$${\left({\int }_0^{\delta }r{u}_r^{2}dr\right)}^{1/2}{\left({\int }_0^{\delta }\frac{u^2}{r}dr\right)}^{1/2}\ge {\int }_0^{\delta }u|u_r|dr>{\int }_0^{\delta }\frac{ϵ}{r}dr=\infty .$$

$$-{\int }_0^{R}r{u}_r^{2}dr\ge {\int }_0^R\left(\beta +\frac{m^2}{r^2}+\frac{1}{1+\alpha u^2}\right)r{u}^{2}dr>\left(\beta +\frac{m^2}{R^2}\right){\int }_0^{R}r{u}^{2}dr.$$

$$\frac{R^2}{r_0^2}{\int }_0^{R}r{u}_r^{2}dr\ge {\int }_0^{R}r{u}^{2}dr$$

$$-\frac{{r_0}^2}{R^2}{\int }_0^{R}r{u}^{2}dr$$

$$\eta u_{rr}(\eta )-\frac{m^2}{{\eta }^2}\eta u(\eta )-\frac{\eta u(\eta )}{1+\alpha u(\eta {)}^2}$$

$$u_{rr}(\eta )$$

$$0\ge \left(\beta +\frac{m^2}{{\eta }^2}+\frac{1}{1+\alpha u(\eta {)}^2}\right)u(\eta ).$$

$$u(\eta {)}^2\ge -\frac{1}{\alpha }\left(\frac{1}{\beta +\frac{m^2}{{\eta }^2}}+1\right)>-\frac{1}{\alpha }\left(\frac{1}{\beta +\frac{m^2}{R^2}}+1\right).$$

$$\beta >-\frac{m^2}{R^2}-1,$$

$$\Delta =\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial }{\partial r}\right).$$

$$\Delta u(r)$$

$$\Delta u(r{)}^2$$

Full text

1 Introduction

Vortices has long been intensely studied in may areas in physics and mathematics. In particular, quantum vortices, also known as Abrikosov vortex or fluxon, are vortex solutions of $2$- or $3$-dimensional wave equations first predicted by Lars Onsager in the context of superfluid. Mainly developed by Alexi Abrikosov, vortex solutions later arose from the Ginzburg-Landau equations in the context of superconductivity and served as a keystone in the explanation of type–II superconductors using the Ginzburg-Landau equations which describes the phase transition of superconductors using variation principle.

In the studies of optical beams, Chiao, Garmire and Townes [32] developed their monumental work in 1964 on the nonlinear propagation of light. They derived the nonlinear Schrödinger equation from simplifying the Ginzburg-Landau equations in $1+1$-dimensions and demonstrated the existence of self-reinforced wave packets, later known as “solitons” [29]. It is the phenomenon of nonlinear propagation of high-intensity light beams enabling the light beams to produce their own waveguide and propagate without spreading. They then proposed that that, in $2$- or $3$-dimensions, the the propagation wave vortices could also shown self-reinforced or soliton behaviors.

Intuitively, vortex solitons are phase singularities where energy flows around a point. Approaching the center of the vortex, the velocity reaches infinity, and yet due to the finite energy flux, the intensity of the center will vanish [2]. One may essentially picture this kind of ring-profiled vortex soliton as a ring of light around a black spot directed along its axis of propagation [37]. Such vortices have widespread applications in many areas of mathematics and physics, including condensed matter physics, particle interactions, cosmology, superconductivity, quantum information processing, and wireless communication [37, 21, 25, 1, 2, 39, 11, 19, 13], and have since been intensively studied in optics, both theoretically and experimentally [16, 26, 43, 17, 8, 27, 4], and observed in various nonlinear media [21, 25, 11, 19, 13].

In this paper, we establish a rigorous mathematical theory for a new type of nonlinearity encountered in an experiment reported by Fleischer et al. [20], known as a photorefractive nonlinearity [28, 7]. Photorefractive effect essentially describes the behavior of certain materials responding to light by altering their refractive index [15]. Similar applied analysis papers have pursued the existence of this nonlinearity in various media in an attempt to determine the geometrical properties and constraints of the potentially observable vortex solitons, and a number of bounds have been derived for the energy flux and propagation constants [36, 24, 23, 37, 5]. Proving the theoretical existence of the photorefractive nonlinearity will provide insights on the future direction of experimental research.

In modern optics, light is theoretically described by a complex-valued wave function governed by the nonlinear Schrödinger equations [33, 34, 40, 41, 3, 9, 10, 12]. In the experimental setup used by Fleischer et al. to create photonic lattice solitons via nonlinear optical induction [20], the wave function is written as two coupled nonlinear differential equations:

[TABLE]

In these equations, $\psi$ describes the slowly varying amplitude of the lattice wave and $\phi$ is the soliton-forming probe; both are complex-valued functions. $k_{1},k_{2}$ relate to the anisotropy of the indices of refraction. Here, $\nabla^{2}_{\perp}$ denotes the Laplace operator over the transverse plane of the coordinate. $\Delta n_{1},\Delta n_{2}$ are the nonlinear index changes induced by the total intensity $I=|\phi|^{2}+|\psi|^{2}$. Fleischer et al. chose a photorefractive screening technique in which the change in index is given by

[TABLE]

where $P,Q\neq 0$ denote coupled parameters [30].

As our focus here is a ring-profiled spatial soliton, we expect the field to be written in polar coordinates over $\mathbb{R}^{2}$ using an $m$-vortex ansatz [42, 18] as

[TABLE]

where $\beta_{1},\beta_{2}\in\mathbb{R}$ are the wave propagation constants, $\theta$ represent the phase, and $m_{1},m_{2}\in\mathbb{Z}$ are known as vortex numbers, or topological charges of vortices. $u,v$ are real-valued, radial profile functions which gives rise to the amplitude of the solitons. The presence of the the vortex at $r=0$ requires $u(0)=v(0)=0$, and as in [18], the self-concentrating effect of solitons suggest $u(r),v(r)$ may be assumed to vanish at a large distance $R>0$ [42]. Thus, we here present the boundary condition $v(R)=u(R)=0$ for some $R>0$.

Using the changes in index given by (1.4), we arrive at a boundary value problem for the coupled nonlinear differential equations

[TABLE]

satisfying the boundary conditions

[TABLE]

In this paper, we prove the existence of semitrivial solutions of the coupled system (1.8)–(1.9) by simplifying the coupled system into a single ordinary differential equation: these are cases where $v$ is a scalar multiple of $u$. Here, we use the substitution $v=\sqrt{\alpha-1}\ u$, $\alpha\geq 1$, for convenience. In this case, as we will later prove in Lemma A.1, if both $u,v$ are nonzero, $|m_{1}|=|m_{2}|$ and $\beta_{1}=\beta_{2}$. If either $u,v$ is equal to [math], the problem is already a single equation with $\alpha=1$. Hence, we can simplify (1.8) as a two-point boundary value problem for the nonlinear ordinary differential equation

[TABLE]

with the boundary condition

[TABLE]

with an undetermined parameter $\beta$ and prescribed $R$, for any given $m\in\mathbb{Z}$

Recall the physical meaning of $u$ as the soliton amplitude. We are interested in the nontrivial solutions which are classic, positive solutions of $u$.

Definition 1.1** (Positive Solutions).**

A positive solution of (1.10) and (1.11) means $u\in C^{2}[0,R]$ is a classic solution to (1.10) and (1.11) such that

[TABLE]

By classic solution, we mean $u$ satisfy (1.10) and (1.11) pointwise almost everywhere.

To tackle this problem, we shall use the methods of calculus of variation which will be discussed in section 3. This relies on a constrained minimization approach that views (1.10) as an eigenvalue problem and that the propagation constant $\beta$ arises as a Lagrange multiplier. In this sense, the differential equation can be viewed as the Euler–Lagrange equation of the action functional

[TABLE]

subject to the constraint functional

[TABLE]

where $P_{0}>0$ is fixed and $\beta$ appears via the Lagrange multiplier of the constrained optimization problem

[TABLE]

with $\mathcal{U}$ denoting an admissible class

[TABLE]

and the energy functional defined by

[TABLE]

The functional $P(u)$ has been used as the beam power [41, 38], energy flux [12], or stability integral [35]. In our case, we refer to it as the energy flux. Similar to the method described in [42], it suffices to show that and the optimal solution $u\in\mathcal{U}$ under action functional (1.13) is a positive solution satisfying (1.10) and (1.11) pointwise.

Our major results are as follows:

Theorem 1**.**

For the existence of positive solutions of (1.10)–(1.11), the following conditions are necessary:

[TABLE]

where $r_{0}\approx 2.404825$ is the first zero of the Bessel function $J_{0}$ **[42]**. 2. 2.

[TABLE]

with $\max_{0<r<R}\{u(r)\}$ denoting the peak of the profile curve. 3. 3.

If

[TABLE]

then the solution $u$ is bounded by the exponentially decaying estimate

[TABLE]

as $r\to R^{-}$. Here, $\epsilon_{0}$ is a constant bounded by $\epsilon_{0}>2\left(\beta+\dfrac{m^{2}}{R^{2}}+1\right)>0$, and $C_{\epsilon_{0}}$ depends on $\epsilon_{0}$.

Theorem 2**.**

There exists solution pair $(u,\beta)$ satisfying the differential equations (1.10)–(1.11) pointwise, where $u$ is a positive solution following Definition 1.1 and $\beta\in\mathbb{R}$. Such a solution can be obtained by solving the constrained minimization problem (1.15), from which $\beta$ arises as a Lagrange multiplier.

The remainder of this paper is structured as follows. In section 2, we provide proofs for Theorems 1. In section 3, we prove our Theorem 2, which is the main theorem demonstrating the existence of positive solutions using the constrained minimization approach. Section 4 presents a selection of results from numerical computations using the finite element method and provides a summary of the results. The proof of a lemma used in the paper is provided in the appendix.

2 Proof of Theorem 1

To prove (1.18), multiplying by $u$ and integrating by parts over $r$, we obtain

[TABLE]

We first want to show that the first term is [math]. Recall that we are interested in a positive solution as defined in Definition 1.1. Since $u\in C^{2}[0,R]$, as $r\to R$, $u_{r}$ is uniformly bounded which implies $ruu_{r}\to 0$ as $r\to R$. We will then prove that $ru(r)u_{r}(r)\to 0$ as $r\to 0$. Suppose that $\liminf_{r\to 0}\{ru|u_{r}|\}\neq 0$. Then, there exists some $\epsilon>0$ and $\delta\in(0,R]$ such that $ru|u_{r}|>\epsilon$ for any $r\in(0,\delta)$. Using the Cauchy–Schwarz inequality, we have that

[TABLE]

This is a contradiction with regard to the energy functional, which implies that

$\liminf_{r\to 0}\{ru|u_{r}|\}=0$. Hence, we can establish the following inequality:

[TABLE]

Recall Poincaré’s inequality[22]

[TABLE]

where $r_{0}\approx 2.40483$ is the first zero of the Bessel function $J_{0}$. Together with (2.3), we have that

[TABLE]

which results in $\beta<-\dfrac{m^{2}+r_{0}^{2}}{R^{2}}$, as claimed. ∎

To prove (1.19), let $u\in C^{2}[0,R]$ satisfy boundary condition (1.11). As we are interested in nontrivial solutions of $u$, we know that there exists a maximizer $\eta\in(0,R)$ such that $u(\eta)=\max\{u\}$, which also implies $u_{rr}(\eta)\leq 0$ and $u_{r}(\eta)=0$. At $r=\eta$, (1.10) becomes

[TABLE]

which simplifies to

[TABLE]

Taking the concavity at the peak $u_{rr}(\eta)\leq 0$, we arrive at the inequality

[TABLE]

This simplifies to the following expression for the maximum of $u$, as claimed:

[TABLE]

Furthermore, we see that if an desired inequality

[TABLE]

holds between the prescribed variables $\beta,m,R$, then $\max\{u\}>0$. This can serve as a sufficient but not necessary condition for the existence of nontrivial solution. ∎

To prove (1.20)–(1.21), let $\Delta$ denote the Laplacian operator in polar coordinates, defined by

[TABLE]

Rewriting (1.10) by replacing $(ru_{r})_{r}$ with $r\Delta u$ and simplifying gives

[TABLE]

Using the product rule $\Delta u^{2}\geq 2u\Delta u$, we then obtain

[TABLE]

Because $\lim_{r\to R^{-}}\dfrac{1}{1+\alpha u^{2}}=1$, for all $\epsilon>0$, there exists $R^{\prime}\in(0,R)$ such that

[TABLE]

for all $r\in(R^{\prime},R)$. Substituting (2.14) into (2.13), we then have

[TABLE]

Insert the condition $\beta>-\dfrac{m^{2}}{R^{2}}-1$ and select a small $\epsilon$ such that $\epsilon<\beta+\dfrac{m^{2}}{R^{2}}+1$. It can then be shown that there exists some $R^{\prime}$ such that, for all $r\in[R^{\prime},R]$,

[TABLE]

for some $\epsilon_{0}>0$. Consider the decreasing exponential function $\xi$ with respect to its exponential coefficient $\sqrt{\epsilon_{0}}$ and a constant $C>0$ defined as

[TABLE]

for which the Laplacian is computed to be

[TABLE]

Subtracting (2.16) from (2.17) gives

[TABLE]

Because $u^{2}(r)\geq u^{2}(R)=0$ and $Ce^{-\sqrt{\epsilon_{0}}r}>Ce^{-\sqrt{\epsilon_{0}}R}$ for all $r\in[R^{\prime},R]$, we find that

[TABLE]

This establishes the following inequality with respect to exponential decay:

[TABLE]

We further denote this as $u^{2}(r)<C_{\epsilon_{0}}e^{-\sqrt{\epsilon_{0}}r},\ r\in[R_{\epsilon_{0}},R]$ to emphasize the dependence of $C,R$ on $\epsilon_{0}$, similar to [23]. ∎

3 Proof of Theorem 2

To prove Theorem 2, we employ a variational principle and constrained minimization problem. Recall the definitions of the action functional (1.13), constraint functional (1.14), admissible class (1.16), and the minimization problem (1.15).

The differential equation (1.10) can be viewed as an eigenvalue problem where the propagation constant $\beta$ is undetermined and acts as a Lagrange multiplier. In order to prove the existence of a solution pair $(u,\beta)$, it suffices to show that a solution to the minimization problem (1.15) exists subjected to the prescribed value of energy flux $P_{0}$.

We first shows the existence of feasible solution to the constraint minimization problem (1.15) by showing the existence of $u\in\mathcal{U}$ such that $P(u)=P_{0}$ and $I(u)<\infty$. Since $x\geq\dfrac{1}{\alpha}\ln(1+\alpha x)$ for $x\geq 0,\alpha>0$, we see that

[TABLE]

Consider the following function $f\in\mathcal{U}$ defined as

[TABLE]

Since the function was designed so that

[TABLE]

along with the property $f(0)=f(R)=0$, we see that $f$ belongs to the feasible set of (1.15). By direct computations, we can compute for $I(f)$. Using the following computations

[TABLE]

we arrive at a finite value for $I(f)$

[TABLE]

Since $f$ is a feasible solution to (1.15) with finite objective value, the optimization problem is well defined. Note the lower bound of the action functional $I\geq 0$; the boundedness of $I$ gives us a weakly convergent subsequence $\{u_{k}\}_{k=1}^{\infty}$, which we can use it as minimizing sequence of (1.15) such that

[TABLE]

Using boundedness and the fact that $\{u_{k}\}_{k=1}^{\infty}$ is minimizing, it follows that there exist $C>0$ independent of $k$ such that

[TABLE]

This inequality is further developed into

[TABLE]

Because the distributional derivative of $u$ satisfies $||u|_{r}|\leq|u_{r}|$ [6] and the functionals $I$ and $P$ are both even, we have

[TABLE]

which implies that we may modified minimizing sequence $\{u_{k}\}_{k=1}^{\infty}$ so that it consists of nonnegative valued functions $u_{k}>0$.

By setting $r=\sqrt{x^{2}+y^{2}}$, we can view these functions as radially symmetric functions over the disk $B_{R}=\{(x,y)\in\mathbb{R}^{2}|x^{2}+y^{2}\leq R^{2}\}$ which all vanish at boundary $\partial B_{R}$. Together with (3.6), $\{u_{k}\}_{k=1}^{\infty}$ is shown to be bounded under the radially symmetrically reduced norm, with definition

[TABLE]

for the standard Sobolev space $W_{0}^{1,2}(B_{R})$. Notice that we are using $\{u_{k}\}_{k=1}^{\infty}$ to denote both the minimizing subsequence as well as the radial function view. Since $u_{k}$ are radial symmetric function, we can write it as $u=u(r)$ which satisfy $u(R)=0$. Given that the sequence is monotonically decreasing but bounded in terms of the action functional, it is sufficient to show that $u_{k}$ converges weakly to a minimizer $u$ on $W_{0}^{1,2}(B_{R})$. Here $u$ is used to denote the limit of sequence.

Using the compact embedding $W^{1,2}(B_{R})\subset\subset L^{p}(B_{R})$ for $p\geq 1$, the sequence $u_{k}$ converges strongly to $u$ in $L^{p}(0,R)$. Furthermore, for any $\epsilon\in(0,R)$, $\{u_{k}\}_{k=1}^{\infty}$ is a bounded sequence in the space $W^{1,2}(\epsilon,R)$. If we apply compact embedding $W^{1,2}(\epsilon,R)\subset\subset C[\epsilon,R]$, we see that $u_{k}\to u$ as $k\to\infty$ uniformly over $[\epsilon,R]$.

To show that $u(0)=0$, let $r_{1},r_{2}\in(0,R)$ such that $r_{1}<r_{2}$, we have Cauchy–Schwarz inequality

[TABLE]

where the constant $C$ is provided in (3.5). By letting $k\to\infty$, this becomes

[TABLE]

In view of (3.6) and applying Fatou’s Lemma, we have

[TABLE]

We see that $\dfrac{u^{2}}{r}\in L^{p}(0,R)$. This implies that the right hand side of (3.9) tends to zero as $r_{1},r_{2}\to 0$, which further implies that the limit exists for $u(r)^{2},r\to 0$ and it is

[TABLE]

Thus, the boundary condition for $u(0)=0$ is achieved. Recall (3.1), the $\dfrac{r}{\alpha}\ln(1+\alpha u^{2})$ term is bounded by

[TABLE]

Using the above results, we can conclude that $u$ obtained from minimizing sequence $\{u_{k}\}_{k=1}^{\infty}$ for the minimization problem (1.15) satisfies the desired solution, i.e. positive solution and satisfying the boundary conditions. The functional $I$ is weakly lower semicontinuous

[TABLE]

and the constrains can be achieved as

[TABLE]

Consequently, $\beta$ appears as an Lagrange multiplier.

Since our solution $u$ is nonnegative, suppose there exists $r_{0}\in(r,R)$ such that $u(r_{0})=0$, $r_{0}$ would be a minimizer with $u_{r}(r_{0})=0$. Recall the uniqueness theorem for the initial value problem of ordinary differential equations, we must have the trivial solution $u(r)=0$ for all $r\in(0,R)$ which contradict with the value for energy flux $P_{0}$. Therefore the positive condition $u(r)>0,r\in(0,R)$ is proved.

∎

4 Computations Using Finite Element Method

To approximate solutions for (1.15) and their corresponding eigenvalue $\beta$ given a fixed energy flux $P_{0}$, we employ variational principle and a type of finite element method known as Ritz-Galerkin method to approximate solutions to boundary valued problems in weak formulation.

Viewing the the admissible class $\mathcal{U}$ defined in (1.16) as a vector space, we can describe its elements $u$ as a linear combination of basis functions $\{\psi_{j}\in C^{1}[0,R]\}_{j=1}^{\infty}$ with coefficients $c^{i}\in\mathbb{R}$:

[TABLE]

The vector space is then equipped with the weighted inner product

[TABLE]

to be viewed as a Hilbert Space. To approximate the solution numerically, we will choose a set of Schauder basis[31] so that we can approximate the solution (4.1) using finite dimension $N\in\mathbb{N}$ as

[TABLE]

Since the trigonometric functions can serve as a set of Schauder basis, the specific choice of basis functions we will use for the weighted inner product is

[TABLE]

The orthogonal basis $\{\psi_{j}\}_{j=1}^{N}$ can be obtained by applying Gram–Schmidt process on (4.4) such that $\langle\psi_{i},\psi_{j}\rangle=\delta_{ij}$. To minimize error, we have implemented the modified Gram–Schmidt[14] in code. Note that the energy flux can be computed with the inner product as follows:

[TABLE]

Our code obtains an orthonormal basis by applying the modified Gram–Schmidt process on $\left\{\sin\left(j\pi r/R\right)\right\}_{j=1}^{N}$. As $u$ can be viewed as a minimum of the action functional (1.13), in view of (1.15) and (4.3), we arrive at the optimization problem

[TABLE]

where we minimize $c\in\mathbb{R}^{N}$. We use MATLAB’s Optimization Toolbox and the Chebfun package to numerically solve for (4.6). For convenience, all subsequent numerical computations have $\alpha=1$ as the curves with $\alpha>1$ generally have the same shape. The propagation constant then can be computed as follows: let $\lambda\in\mathbb{R}$ be the Lagrange multiplier, so that $\langle I^{\prime}(u),\tilde{u}\rangle=\lambda\langle P^{\prime}(u),\tilde{u}\rangle$. $I^{\prime}(u)$ and $P^{\prime}(u)$ are given by

[TABLE]

where $\tilde{u}$ is an arbitrary test function. The expression $\langle I^{\prime}(u),\tilde{u}\rangle=\lambda\langle P^{\prime}(u),\tilde{u}\rangle$ can then be expanded as

[TABLE]

Applying integration by parts to the first term, we have

[TABLE]

This gives the weak version of (1.10). We can therefore claim that $\beta=-4\pi\lambda$. We rewrite (4.9) as

[TABLE]

and finally arrive at an expression for $\beta$ that is suitable for numerical computation:

[TABLE]

To compute the error in the numerical computations, taking (1.10) and rearranging into the form

[TABLE]

as the left-hand side is [math], we can use this to calculate the computational error of our eigenvalue by integrating over the square of the right-hand side of (4.13).

[TABLE]

As an example of the numerical results, Fig. 1 shows the profile of the amplitude of the soliton if we fix $R=40$ and $P_{0}=200$ while varying the vortex number $m$. We see that increasing $m$ corresponds to an outward shift of the peak. Figure 2 demonstrates the relationship between the profile function and energy flux $P_{0}$. An increase in $P_{0}$ corresponds to an increase in the height of the peak of the profile curve. Furthermore, beyond some threshold, the profile curves start to become visibly asymmetric over $(0,R)$, which is likely to be the threshold governing the behavior as $R\to\infty$. The reason for this conjecture is that the asymmetry will cause $u$ to vanish near $R$; and thus, for an increase $\Delta R$ in $R$, the increase $\Delta I,\Delta P$ in the corresponding integrals is likely to be small due to the asymmetry over the domain. Physically, $R$ is chosen to be large enough such that $u$ is nearly vanished some distance before $R$ for a self-reinforced soliton [42]. In this view, only the asymmetric results have applications in physics. Recall Theorem 1; the additional condition $\beta>-\frac{m^{2}}{R^{2}}-1$, which implies exponential decay near $R$, has the potential of categorizing solutions for the $R\to\infty$ case. Unfortunately an increase in $R$ drastically increase the computation’s error, dimensions required, and time to be verified with confidence. The $R\to\infty$ case could be studies more in future research.

Some values for propagation constants $\beta$ are also computed. For example, for $R=20$, $m=1$, and $P_{0}=1,50,100,200,400,800$, the numerically computed eigenvalues are $\beta=-1.1419$, $-0.9328$, $-0.7818$, $-0.6185$, $-0.4792$, $-0.3679$ with corresponding errors \Delta\beta=2.9550$$1\text{\times}{10}^{-4}\text{\,}\mathrm{,} 5.7607$$1\text{\times}{10}^{-4}\text{\,}\mathrm{,} 4.4359$$1\text{\times}{10}^{-4}\text{\,}\mathrm{,} 3.6516$$1\text{\times}{10}^{-4}\text{\,}\mathrm{,} 5.6023$$1\text{\times}{10}^{-4}\text{\,}\mathrm{,} $0.0015$. The following table presents some of the computed eigenvalues along with their errors for $R=20$, $m=1$, and various $P_{0}$.

[TABLE]

Finally, Fig. 3 shows $\beta$ as a function of $P_{0}$ for $R=20$ and $m$ = $1$, $2$, $3$, $4$, $5$. The error bars show the computational error computed with (4.14) due to the finite number of dimensions ($N=20$). The corresponding upper bounds, proved with Lemma 1.18, are shown as dotted lines.

Appendix A Additional Lemmas

Lemma A.1**.**

If $u$ and $v$ are nonzero scalar multiples of each other that satisfy positive solutions to (1.8), then $|m_{1}|=|m_{2}|$ and $\beta_{1}=\beta_{2}$.

• Proof

Suppose that both $u$ and $v$ are nonzero and positive. By making one a scalar multiple of the other, we can set $v=u\sqrt{\alpha-1}$, $\alpha>1$, and then (1.8) reduces to

[TABLE]

Subtracting (A.1) from (A.2), we can simplify the difference as

[TABLE]

Given that $m_{1},m_{2},\beta_{1},\beta_{2}$ are constants while there exist more than one different $r^{2}u>0$ in the domain, the only way to make this equation hold for all $r\in(0,R)$ is to have [math] on both sides, implying that $|m_{1}|=|m_{2}|$ and $\beta_{1}=\beta_{2}$, as claimed.