Quantum harmonic oscillators with nonlinear effective masses in the weak density approximation

We study the eigen-energy and eigen-function of a quantum particle acquiring the probability density-dependent effective mass (DDEM) in harmonic oscillators. Instead of discrete eigen-energies, continuous energy spectra are revealed due to the introduction of a nonlinear effective mass. Analytically, we map this problem into an infinite discrete dynamical system and obtain the stationary solutions in the weak density approximation, along with the proof on the monotonicity in the perturbed eigen-energies. Numerical results not only give agreement to the asymptotic solutions stemmed from the expansion of Hermite-Gaussian functions, but also unveil a family of peakon-like solutions without linear counterparts. As nonlinear Schrödinger wave equation has served as an important model equation in various sub-fields in physics, our proposed generalized quantum harmonic oscillator opens an unexplored area for quantum particles with nonlinear effective masses.


Introduction
Quantum harmonic oscillator is the most important model system in quantum mechanics, which remarkably exhibits an exact, analytical solution with discrete (quantized) eigen-energies compared to the predictions of classical counterparts [1]. Instead with a given mass, m 0 , when particles (electrons or holes) move inside a periodic potential or interact with other identical particles, their motions differ from those in a vacuum, resulting in an effective mass [2]. With an effective mass, denoted as m * , the corresponding Schrödinger equation for a quantum particle in a one-dimensional harmonic oscillator, characterized by the spring constant k, has the form: position-dependence, chromatic dispersion may also have intensity-dependent dispersion (IDD) in the optical domains [11,12]. IDD, or in general the nonlinear corrections to the chromatic dispersion as a function of the wave intensity, has arisen in a variety of wave phenomena, such as shallow water waves [13,14], acoustic waves in micro-inhomogeneous media [15], ultrafast coherent pulses in quantum well waveguide structures [16], the saturation of atomic-level population [17], electromagnetically-induced transparency in a chain-Λ configuration [18], or nonlocal nonlinearity mediated by dipole-dipole interactions [19]. Inspired by IDD, in this work, we consider a quantum particle acquiring an probability density-dependent effective mass (DDEM), i.e., m * (|Ψ| 2 ), in a harmonic potential described by the following generalized Schrödinger equation: Here, the DDEM is approximated in the weak density condition by assuming b|Ψ| 2 < 1. Then, we have with the parameter b denoting the contribution from the nonlinear effective mass term. As one can see, when the nonlinear effective mass term is zero, i.e., b = 0, equation (2) is reduced to the well-known scenario for a quantum particle in a parabolic potential. Besides, optical waves IDD in the Gradient-index (GRIN) lenses [20] also share the same model equation described in equation (2) when we consider wave propagation along the z-coordinate (by replacing t by z), along with a confined transverse dimension denoted by x.
When k = 0, even though an approximated solution with a non-physical peak value 10 12 was illustrated in Ref. [11] for b < 0, it is proved rigorously that localized solitons exist only for b > 0 in Ref. [12]. However, with the supported harmonic potential, the scenarios can be totally different. In terms of solitons, the interplay between nonlinearity and harmonic potential has been studied for a long time [21]. Nevertheless, such localized solutions supported only with a nonlinear effective mass is never studied. Moreover, in the weak density approximation, equation (2) also shares the same mathematical form to the Salerno model in the continuum limit, which can be derived as a quantum modified discrete nonlinear Schrödinger equation, giving the time evolution of the field amplitude on the lattice [22][23][24]. Without considering any nonlinear potential but only with the nonlinear effective mass differs our results from the known ones.
However, when b ≠ 0, instead of the discrete energies, continuous energy spectra are revealed due to the introduction of a nonlinear effective mass. Analytical solutions for the corresponding eigen-energy and eigenfunction are derived by expanding the solutions in the weak density, i.e., b|Ψ| 2 . Numerical solutions obtained by directly solving equation (2) give good agreement to the analytical ones obtained from the expansion of Hermite-Gaussian functions. Moreover, we unveil a family of peakon-like solutions supported by DDEM, which has no counterpart in the linear limit. Our perturbed solutions and numerical results for this generalized quantum harmonic oscillator with nonlinear effective masses opens an unexplored area for quantum particles.
The paper is organized as follows: in Session II, we introduce the quantum harmonic oscillator into this generalized Schrödinger equation with nonlinear effective mass and reduce equation (2) into an infinite dynamical system. Then, by expanding b|Ψ| 2 and with the help of the eigen-solutions of quantum harmonic oscillator, we study the corresponding eigen-energy with the introduction of DDEM, as a function of the parameter b. The monotonicity of the perturbed eigen-energy is also proved. In section 3, explicitly, we derive the analytical solutions of eigen-energies and the corresponding wavefunctions for the ground and second-oder excited states in the asymptotical limit. The comparison between analytical solutions and numerical results is illustrated in section 4, demonstrating good agreement on the solutions with a smooth profile, stemmed from the expansion of Hermite-Gaussian wavefunctions. A new family of peakon-like solutions with a discontinuity in its first-order derivative is also unveiled, which has no linear counterparts. Finally, we summarize this work with some perspectives in Conclusion.

Quantum harmonic oscillator with DDEM
Without loss of generality, in the following, we set ÿ = 1, k > 0, m 0 = 1 for the simplicity in tackling equation (2). Here, by lookin for the stationary solutions Ψ(x, t) = ψ(x)e − iE t , we consider a family of differential equations parametrized by a continuous DDEM parameter b ≠ 0 of the form where E is the corresponding eigen-energy, x denotes a real variable for the coordinate, and ψ(x) is a square integrable function. This stationary Schrödinger wave equation can be seen as a generalized quantum harmonic oscillator. In addition to the stationary states considered here, this generalized nonlinear Schröoedinger equation is expected to provide an interesting platform, with non-stationary states, for further studies on the infinite dynamical system.
When b = 0, equation (3) becomes the well-known equation for the quantum harmonic oscillator, which supports eigen-function of the n-th order excited state in the position representation reads [25]: 2.1. Expanding the eigen-energies and eigen-functions with b|Ψ| 2 To investigate equation (2) with b ≠ 0 (but keeping k ≠ 1 first), we look for the solutions based on the expansion of the solution on the eigen-function f n (x). That is, .  where V m,n,p,q , and W m,n,p,q are defined as: As one can see from equation (7), we reduce the original partial differential equation into the infinite discrete dynamical system [26]. With the help of the recursive relation of Hermite polynomial H m (x), for example see Ref. [25], i.e., Now, we look for the stationary solution in the form: From now on, for simplicity, we assume k = 1. With the help of equation (10), next, we consider the perturbation on the energy deviated from the eigen-energy E n with the corresponding Hermite-Gaussian eigenmode f n (x).
Similar to the methodology used in dealing with the nonlinear mean field in the Gross-Pitaevskii equation (9), and arrive at a nonlinear eigen-energy equation: Once again, in equation (11), we can see that if P(E) → 0, then the resulting eigen-energy  = + E E n n 1 2 . By substituting f(x), obtained from equation (10), into equation (11), one can have the relation between E and P(E) near = + E n n 1 2 . In general, the perturbation approach illustrated above works for all the values of n. However, only when n is even, a neat formula can be conducted by taking the advantage of symmetric wavefunctions in ψ(x). For even numbers, 2n, the resulting eigen-energy E due to the introduction of the DDEM parameter b can be approximated as To compute the integrals shown in equation (12), one can utilize the Feldheim identity for the Hermite polynomials [25]: where n + p m, m + p n and m + n p; otherwise, the integral is zero. From equations (13) and (14), a direct calculation can yield

Monotonicity in the perturbed eigen-energy
), we prove that the two integrals inside the square brackets in equation (12) is monotonic, i.e., or equivalently n n n n n n n n 2 ,2 ,2 ,2 2 ,2 ,2 ,2 for n 0. In Appendix, the proof on the monotonicity for equation (16) and equation (17) is given in details. By using the upper and lower solution method developed in the variational calculus [27], we can further prove the existence of a positive solution (node-less state) through the corresponding Lagrangian for equation (2), i.e., ( ) Here, we also introduce the probability factor Q for this quantum harmonic oscillator with DDEM, by defining As the original generalized Schrödinger equation given in equation (2) preserves the U(1) symmetry, i.e., [ ] y q y  i exp , the conserved density for this model equation can be derived from Noether theorem [30]. It is marked that the local conservation law is also valid, i.e., [ ( . It is noted that equation (20) is only applicable when b ≠ 0. When |b| = 1, this probability factor Q can be approximated as which is reduced to the standard definition of probability for quantum wavefunctions. For b → 0, the corresponding Lagrangian density given in equation (18), as well as the conserved density given in equation (20), both go to |ψ| 2 . These two terms, Z(E) and Q(E), shown in equation (19), correspond to the Lagrangian of our generalized harmonic oscillator and the conserved quantity, respectively. As the DDEM parameter b → 0, the Lagrangian shown in equation (19) can be reduced to which is the Lagrangian for the linear equation, i.e., y y y - . By following the same concept in tackling weak non-linearity [31], the perturbation based on the expansion of the Hermite-Gaussian functions to deal with the DDEM ensures that when E → E n , one has Q(E) → 0.
3. Eigen-energies and eigen-functions obtained from perturbation 3.1. The ground state Now with the analytical formula give in equation (12), we explicitly give the perturbed eigen-energy E and eigenfunction for the ground state in our generalized quantum harmonic oscillator with a given DDEM parameter b. For the ground state, we can assume that B 0 ? B 2n , n = 1, 2, 3, L . Then, from equation (10), one has for n 1. Therefore, from equations (22) and (23), explicitly we have, noting that = With the coefficients above, the perturbed solution of ( ) y x b 0 can be conducted immediately as We notice that ( ) Again, with the orthonormality of f 2n (x), in the asymptotical limit, n → ∞ , the probability factor Q defined in equation (20) becomes: , by assuming B 2 ? B 2n , n = 0, 2, 3, L . Again, with equations (8) and (10), one can directly obtain: As a result, we we have . It is noted that here, the expansion starts from n = 2 as B 0 = 0. Again, we have ( ) 2 . Moreover, thee resulting probability factor Q(E) in the asymptotical limit, n → ∞ has the form:

Numerical results by direct simulations 4.1. The ground state
In order to verify the validity of our analytical solutions obtained by the expansion, we also perform the numerical calculations for equation (3) directly without applying any approximation. Here, the eigen-values/ eigen-functions are generated by substituting the eigen-function iteratively until a convergent eigen-value is reached. Moreover, the linear stability method is applied for the found eigen-function. All the found eigenfunctions are stable due to the harmonic potential. To maintain some level of formal rigor and mathematical correctness, we shall talk about finding solutions of differential equations [32]. To find the solutions of the eigenvalue problem with the nonlinear term, we connect with a quantum harmonic oscillator by solving equation (3) with Fourier spectral method. Using the matrix elements, we diagonalize the matrix numerically and perform the iteration to ensure that the truncated Fourier basis makes the eigen-value convergent. For low energy states, already the smallest basis of 512 elements gives more than sufficient accuracy.
In figure 1, we show the corresponding lowest eigen-mode of the generalized quantum harmonic oscillator described in equation (3), in the plot of probability factor versus eigen-energy Q-E. Starting from E 0 = 0.5, i.e., the eigen-energy of ground state in the standard quantum harmonic oscillator with b = 0; however, the eigenenergy is not a discrete value, but a continuous function due to the introduction of DDEM, i.e., b ≠ 0. Here, the initial guess solution has a single-hump profile, i.e., a Gaussian function stemmed from the zero-th order H n (x). With a positive value of b, such as b = 1 and b = 2, shown in Blue-and Red-colored curves in figure 1, the corresponding probability factor Q(E) presents an almost linear function of the eigen-energy E. Now, all the eigen-energy E 0 b are larger than that of E 0 . Compared to the analytical formula of ( ) have similar Gaussian shapes both for b = 1 and b = 2. But the resulting amplitude, as well as the width, is smaller with a larger value in the DDEM parameter, such as b = 2. The analytical solutions obtained by perturbed theory, depicted in dashed-curves in figure 2(a), also reflect this similarity.
However, when b is negative, there are two distinct regions in this Q-E curve, illustrated in the Green-and Yellow-colored backgrounds in figure 1. For the Green-colored region, the corresponding eigen-energy is smaller than E 0 = 0.5, but remains positive, i.e., < < E E 0 b 0 0 . The probability factor ( ) Q E b 0 is also linearly proportional to the eigen-energy E, as predicted by our analytical formula in equation (27). But, now the slope of Q-E curve changes its sign, as the b < 0. The resulting wavefunction ( ) y x b 0 , as shown in figure 2(b) for the marked eigen-energy E B = 0.25, still has a smooth profile. However, the corresponding width of wavefunction shrinks when E → 0. As a result, a singularity emerges at = E 0 s 0 for the ground state, in which no well-defined localized wavefunction can be supported. The singularity comes from the divergence of Q(E) near 1 + b|ψ| 2 = 0. Moreover, as one can see, our theoretical formula also breaks down when E approaches this singularity.
Unexpectedly, single-hump solutions can be supported even when < = E E 0 s 0 , as shown in the Yellowcolored region. As shown in figure 2(c) for the marked eigen-energy E C = − 0.25, instead of a smooth profile stemmed from the Gauss wavefunction, the resulting wavefunction of this family of peakon-like solutions has a discontinuity in its first-order derivative, similar to the peakon solution in the form of ( | |) x exp . Such peakonlike solutions are also already found in the IDD setting for optical waves, even without the introduction of harmonic oscillators [11,12]. As our perturbation theory starts from the eigen-basis of Hermite-Gaussian functions, it is not applicable to this family of peak-like solutions.

The excited states
In addition to the ground state, the founded second order excited states, both numerically and analytically, are also depicted in have three humps in their profiles and share the similar profile as the 2nd order Hermite-Gaussian function. By comparing the solid-and dashed-curves, corresponding to our numerical results and analytical solutions, respectively, one can see nearly perfect agreement for the solutions around the eigen-energy E 2 .
Moreover, a discontinuous profile emerges when b < 0 and < E E s . When the eigen-energy is smaller than the value at the first singularity E s 2 ,1 but larger than the value at the second singularity E s 2 ,2 , for example E F = 0.75, the peakon-like solution illustrated in Blue-color in figure 2(f), has a profile of ( | |) x exp in two of the humps in the sidebands. It is noted that the profile in the central hump remains a smooth function. Nevertheless, when the eigen-energy is smaller than the value at the second singularity E s 2 ,2 , for example E G = − 0.25, the corresponding eigen-function has discontinuities in all the three humps, as depicted in Red-color in figure 3(f).
In figure 3, we plot all the founded eigen-energies, up to n = 4, by depicting the solution family with the same number of humps in the eigen-functions ( ) y x n b in the same colors. One can see clearly that, all the Q-E curves start from the eigen-energies = + E n n 1 2 of a standard quantum harmonic oscillator, i.e., b = 0. Around these energy values, E n , our perturbation theory works perfectly, giving the linear dependence of Q(E) on the eigenenergy, along with the inversely proportional relation to the DDEM parameter b. In particular, as depicted in the dashed-curves, it can be seen clearly that our analytical solutions given in equation (33) illustrate good agreement to the numerical solutions near E 2 = 2.5.
However, when b turns negative and the supported eigen-energy E n b is away from the starting energy value is the same, i.e., equal to n + 1. Before Conclusion, we remark the stability of the founded eigen-solutions of our generalized quantum harmonic oscillator with a probability density-dependent effective mass (DDEM). As confined by the external harmonic oscillator, all the found eigen-solutions are stable numerically. The validity of our perturbation theory is limited to the eigen-energy around the known one = + E n 2 n 2 1 2 . It is expected that our analytical formula breaks down when E n b approaches E n s . As for the possible bifurcation maps, how to develop an analytical method to find the solutions for these peakon-like solutions, as well as around the singularities, is a challenge, which goes beyond the scope of the current work but deserves further studies. are depicted in the same colors, with the labelled starting eigen-energies E n , for n = 0, 1, 2, 3, 4. Analytical solutions based on the perturbation theory given in equation (33) are also depicted in the Black dashed-curves, which illustrate good agreement to the numerical solutions near E 2 = 2.5.

Conclusion
By introduction a probability density-dependent effective mass (DDEM) for a quantum particle in harmonic oscillators, we propose a generalized Schrödinger equation to embrace the nonlinear effective mass. With the help of orthonormal property of Hermite-Gaussian functions, we reduce this partial differential equation into an infinite discrete dynamical system and find the corresponding stationary solution by expanding b|Ψ| 2 . The monotonicity of perturbed solutions is also approved rigorously. The resulting eigen-energy spectra is no longer discretized, but continuous due to the introduction of a nonlinear effective mass. With the comparison to numerical results obtained by direct simulations, the validity of our analytical formula in the asymptotic limit, in terms of the probability factor as a function of the eigen-energy, Q(E), can be easily verified, in particular for the solutions stemmed from the expansion of Hermite-Gaussian functions. However, the nonlinear effective mass also introduces a new family of peakon-like solutions with a discontinuity in their first-order derivative, which definitely deserves further studies.
It is noted that what we illustrated in this work is based on tackling equation (2), under the weak density approximation. It is also possible to go beyond the weak density approximation by directly studying equation (1). However, according to quantum mechanics, the effective mass m * (x) does not commute with the should be taken into account when the quantum particle is considered. It has been well studied with the nonlinear Schrödinger wave equation, or the Gross-Pitaevskii equation in general, where the nonlinear terms come from Kerr-effect, or the mean-field interaction. With the eigen-energy and eigen-function illustrated in this work, our proposed generalized quantum harmonic oscillator opens an unexplored area for quantum particles with nonlinear effective masses. A number of promising applications and directions for further exploration may be identified when particles accessing nonlinear correction to their effective mass. Similar models related to our proposed generalized quantum harmonic oscillators, but in more complicated settings involve off-resonant self-induced transparency (SIT) solitons [33,34] spatially-periodic refractivity doped with two-level systems (TLS) [35,36], electromagnetically-induced transparency (EIT) via resonant dipole-dipole interactions [37,38], and the continuum limit of the Salerno model [22][23][24].

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Appendix
Here, we give the details to prove the inequality shown in equations (16) and (17).
First of all, from equation (8), one can see that