Quantization and Bifurcation beyond Square-Integrable Wavefunctions

Probability interpretation is the cornerstone of standard quantum mechanics. To ensure the validity of the probability interpretation, wavefunctions have to satisfy the square-integrable (SI) condition, which gives rise to the well-known phenomenon of energy quantization in confined quantum systems. On the other hand, nonsquare-integrable (NSI) solutions to the Schrödinger equation are usually ruled out and have long been believed to be irrelevant to energy quantization. This paper proposes a quantum-trajectory approach to energy quantization by releasing the SI condition and considering both SI and NSI solutions to the Schrödinger equation. Contrary to our common belief, we find that both SI and NSI wavefunctions contribute to energy quantization. SI wavefunctions help to locate the bifurcation points at which energy has a step jump, while NSI wavefunctions form the flat parts of the stair-like distribution of the quantized energies. The consideration of NSI wavefunctions furthermore reveals a new quantum phenomenon regarding the synchronicity between the energy quantization process and the center-saddle bifurcation process.


Introduction
In the statistical formulation of quantum mechanics, a wavefunction ψ has to be square integrable (SI) to ensure the qualification of ψ * ψ as a probability density function. SI solutions to the Schrödinger equation can be used to determine the energy levels in a confined system. Nonsquare-integrable (NSI) solutions to the Schrödinger equation otherwise are ruled out and their role has been unknown till now. To investigate the role of NSI wavefunctions, we need a formulation of quantum mechanics, which does not require the SI condition. Among the nine different formulations of quantum mechanics [1], there is a formulation known as the quantum Hamilton-Jacobi (H-J) formalism [2,3], which meets our purpose. The quantum H-J formalism has been developed since the inception of quantum mechanics along the line of Jordan [4], Dirac [5] and Schwinger [6]. The main advantage of the classical H-J formalism is to give the frequencies of a periodic motion directly without solving the equations of motion. Analogous to its classical counterpart, the advantage of quantum H-J formalism is recognized as a method of finding energy eigenvalues directly without solving the related Schrödinger equation.
Based on the quantum H-J equation, Leacock and Padgett [2,3] proposed an ingenious method to evaluate energy eigenvalues by contour integral. This approach to energy eigenvalues E n is entirely independent of whether the related wavefunction is SI or not, and allows us to examine the participation of the NSI wavefunctions in the process of energy quantization. Apart from providing energy eigenvalues, quantum H-J formalism like its classical counterpart produces quantum Hamilton dynamics [7], from which complex quantum trajectories can be solved to describe the quantum motion associated with a given wavefunction. Probability interpretation isolates SI wavefunctions from NSI on the complex plane. As the third goal of the paper, we will contrast quantum trajectories of SI wavefunctions with those of NSI wavefunctions to manifest the invisible spin degree of freedom as a rotational motion on the complex plane.
The remainder of this paper is organized as follows: Section 2 presents quantum H-J formalism and the related method of determining energy eigenvalues. In Section 3, time average operation for NSI wavefunctions along a complex quantum trajectory is developed from the quantum H-J formalism to replace the ensemble average. The proposed time average operation is then used in Section 4 to derive the universal quantization laws regarding the kinetic energy and the quantum potential. Section 5 demonstrates the participation of NSI wavefunctions in the energy quantization process for a harmonic oscillator. Section 6 proposes a quantum dynamic description of energy quantization, in terms of which a new phenomenon regarding the synchronicity between quantization and bifurcation is revealed. Finally, both SI and NSI solutions to the Schrödinger equation are considered in Section 7 and their relations to spin degree of freedom are explained.

Quantum Hamilton-Jacobi Formalism
While the quantum H-J theory is general, here we consider its application to bound states, which have quantized energy levels and closed quantum trajectories. The quantum H-J approach to determining energy eigenvalues can be conceived of as an extension of the Wilson-Sommerfeld quantization rule [26]. In this approach, the quantum energy levels are given exactly by setting the quantum action variable equal to an integer multiple of Planck constant: where p(x) is called quantum momentum function (QMF) and the contour c is defined on the complex plane with the integer n being the number of poles of p(x) enclosed by c. The QMF p(x) is related to the quantum action function S and the wavefunction ψ as: with S satisfying the quantum H-J equation: ∂S ∂t + H(t, x, p)| p=∂S/∂x = ∂S ∂t and with ψ satisfying the Schrödinger equation: It appears that the quantum H-J Equation (3) and the Schrödinger Equation (4) are equivalent expressions via the relation S = −i ln ψ.
Leacock and Padgett [2,3] proposed an ingenious method to evaluate J(E) without actually solving p(x) from the quantum H-J Equation (3). They showed that for a given potential V(x), J(E) can be computed simply by a suitable deformation of the complex contour c and the change of variables in Equation (1). Once J(E) is found, the energy eigenvalues E n can be determined by solving E in terms of the integer n via the relation J(E) = n .
The quantum H-J approach to determining energy eigenvalues has two significant implications. Firstly, this approach suggests that the energy eigenvalue E n stems from the quantization of the action variable J, rather than from the quantization of the total energy E itself. Precisely speaking, the energy eigenvalue E n is the specific energy E at which the action variable J(E) happens to be an integer multiple of , i.e., J(E n ) = n . Inspired by this implication, the first goal of this paper is to reveal the internal mechanism causing the quantization of the action variable J and find out its relation to the energy quantization.
Secondly, the quantum H-J approach implies that the SI condition is not required throughout the process of determining energy eigenvalues, which means that whether wavefunctions are SI or not is unconcerned upon evaluating eigen energies. Based on this observation, our second goal here is to expose how SI and NSI wavefunctions cooperate to form the observed energy levels within a confining potential. For a given wavefunction ψ(t, x) either SI or NSI, the associated quantum dynamics can be described by the quantum Hamilton equations with the quantum Hamiltonian H given by Equation (3): where Q(x) is the complex quantum potential defined by: Quantum potential Q(t, x) is intrinsic to the quantum state ψ(t, x) and is independent of the externally applied potential V(x). The quantum Hamilton Equations (5) are distinct from the classical ones in two aspects: the complex nature and the state-dependent nature. The complex nature is a consequence of the fact that the canonical variables (x, p) solved from Equation (5) are, in general, complex variables. The state-dependent nature means that Equation (5) governs the quantum motion specifically in the quantum state described by ψ. The Hamilton Equations (5), which is usually regarded as the complex-extension of Bohmian mechanics, can be derived independently by the optimal stochastic control theory [27].
For a given wavefunction ψ(t, x), the complex contour c traced by x(t) can be solved from Equation (5a), along which the contour integral in Equation (1) then can be evaluated. The second Hamilton Equation (5b) is an alternative expression of the Schrödinger Equation (4) as can be shown by substituting p(t, x) from Equation (2) and Q(t, x) from Equation (6).

Time Average along a Complex Quantum Trajectory
The necessity of considering time average along a complex quantum trajectory comes from the fact that the action variable J introduced in Equation (1) is equal to the time-average kinetic energy, as will be shown below. For a particle confined by a time-independent potential V(x), we have wavefunction ψ(t, x) = e −iEt/ ψ E (x) and quantum action function S(t, x) = −Et − i ln ψ E (x), with which the quantum H-J Equation (3) can be recast into the following form: This is the energy conservation law in the quantum H-J formalism, indicating that the conserved total energy E comprises three terms: the kinetic energy E k = p 2 /2m, the applied potential V(x), and the quantum potential Q(x). When expressed in terms of the wavefunction ψ E (x), Equation (7) and Equation (4) become the time-independent Schrödinger equation: The energy conservation law (7) is valid for any solution ψ E (x) to the Schrödinger Equation (8), either SI or NSI. The Schrodinger Equation (8) has a continuum of solutions, unless it is supplemented with appropriate boundary conditions. Without loss of generality, we consider V(x) in the form of a potential well with the property V(x) → ∞ , as x → ±∞ . Due to the presence of the infinite potential, the probability of finding the particle at infinity is zero, i.e., ψ E (x) → 0, as x → ±∞. (9) This boundary condition gets rid of most of the solutions to Equation (8) and selects out only a discrete set of ψ E and E. Consequently, it is the boundary conditions in standard quantum mechanics that actually enforce the quantization. The boundary condition (9) originates from the fundamental requirement that wavefunctions must be SI, i.e., which allows the normalization of the total probability to unity. If the SI condition (10) or the boundary condition (9) is released, the total energy E will be still conserved, but no longer quantized, because the participation of NSI wavefunctions ψ E (x) in Equation (7) will result in an arbitrary total energy E other than E n . However, even if the total energy E is allowed to be varied continuously, there exist intrinsic quantization laws from which the energy eigenvalue E n can be recovered. In other words, probability interpretation with the accompanying SI condition is not the only way to arrive at the quantization. This issue was first addressed by Leacock and Padgett [2,3] and demonstrated in detail by Bhalla [28,29].
Although NSI wavefunctions ψ fail to serve as probability density functions in the assemble average Ω(x, p) ψ = ψ|Ω(x,p)|ψ , complex quantum trajectories for NSI wavefunction still exist, along which time average of Ω(x, p) can be defined to substitute for the assemble average Ω(x, p) ψ . The complex quantum trajectory describing the particle's motion in a confined system can be solved from Equations (5a) and (2), which together with ψ(t, x) = e −iEt/ ψ E (x) gives the governing equation as: where ψ E (x) is a general solution to Equation (8) with given energy E. The resulting complex trajectory x(t) serves as a physical realization of the complex contour c appearing in Equation (1), and allows the contour integral to be evaluated along the particle's path of motion. By treating the quantum Hamiltonian H(x, p) defined in Equation (7) as a Lyapunov function, the energy conservation law dH/dt = 0 implies that the autonomous nonlinear system (11) is Lyapunov stable (neutrally stable) with equilibrium points in the form of centers, irrespective of whether ψ E is SI or not. The trajectory solved from Equation (11) coincides with the Lyapunov contour lines defined by H(x, p) = E = constant, which are concentric curves surrounding equilibrium points.
The time average of Ω(x, p) along the particle's trajectory x(t) is defined as: where T is the period of oscillation of x(t). The quantum action variable J defined in Equation (1) is a ready example of taking time average along a complex contour. Letting c be a closed trajectory solved from Equation (11), we can rewrite the contour integral (1) in terms of the time-average kinetic energy as: where ω = 2π/T is the angular frequency of the periodic motion. Therefore, the Wilson-Sommerfeld quantization law J = n is simply an alternative expression of the energy quantization law E k T = n( ω/2). In general, the time average of an arbitrary function Ω(x, p) can be expressed in terms of a contour integral by using Equations (11) and (12): where c is the closed contour traced by x(t) on the complex plane, and the symbol "prime" denotes the differentiation with respect to x. Since QMF p(x) can be expressed as a function of x, we simply write Ω(x, p) as Ω(x) in the integrand. According to the residue theorem, the value of Ω(x, p) T is determined only by the poles of the integrand enclosed by the contour c and is independent of the actual form of c. We will see below that the discrete change of the number of poles in the integrant leads to the quantization of Ω(x, p) T .

General Quantization Laws without SI Condition
Let ψ E (x) be a general solution to the Schrödinger Equation (8) with a given energy E. We can treat the time average Ω(x, p) T as a function of the total energy E by noting that Ω(x, p) T is computed by Equation (14) with wavefunction ψ E (x), which in turn depends on the energy E. The time average Ω(x, p) T is said to be quantized, if its value manifests a stair-like distribution as the total energy E increases monotonically. We will derive several energy quantization laws originating from such a stair-like behavior of Ω(x, p) T , which are universal for all confined quantum systems. The energy quantization defined here denotes the discrete change of the considered energy, which is different from the definition in standard quantum mechanics, where energy quantization denotes the discrete energies satisfying the SI condition (10).
Firstly, we consider the quantization of the time-average kinetic energy. By substituting Ω(x, p) = p 2 /2m into Equation (14), we obtain: To evaluate the above contour integral, we recall a formula from the residue theorem: where Z f and P f are, respectively, the numbers of zero and pole of Ω(x) enclosed by the contour c. Using this formula in Equation (15) yields: where the integer n ψ = Z ψ − P ψ is the difference between the numbers of zero and pole of ψ E (x). It appears that the time average of the particle's kinetic energy in a confined potential is an integer multiple of ω/2. This is an universal quantization law independent of the confining potential V(x). Using Equation (17) in Equation (13), we recover the Wilson-Sommerfeld quantization law J = n ψ . The other quantized energy is the quantum potential Q. The evaluation of Equation (14) with Ω(x, p) = Q(x) gives: Applying Formula (16) once again, we arrive at the second energy quantization law: where integer n p = Z p − P p is the difference between the numbers of zero and pole of p(x). Like the quantization of E k T , Equation (19) reveals that the value of Q T is an integer multiple of ω/2, irrespective of the confining potential V(x). The Kinetic energy E k and the quantum potential energy Q, individually, are quantized quantities, and their combination leads to another quantization law. This can be verified from the combination of Equations (15) and (18): where in the integrand can be simplified further as: With the above simplification and the Formula (16), Equation (20) yields a new quantization law: where integer n ψ = Z ψ − P ψ is the difference between the numbers of zero and pole of ψ E (x). The three integers, n ψ , n p and n ψ , are solely determined by the wavefunction ψ E , which in turn is solved from Equation (8) with a prescribed energy E. As E increases, the three integers can only change discretely in response to the continuous change of E. Let E 0 < · · · < E n−1 < E n < · · · be the sequence of specific energies at which the integer n ψ experiences a step jump, n − 1 → n . With increasing E, the value of E k + Q T then assumes a stair-like distribution described by: The values of E k T and Q T have a similar distribution. It is noted that the wavefunction ψ E (x) solved from Equation (8) with an energy E in the interval E n−1 < E ≤ E n is generally NSI. Our next task is to clarify the roles of these NSI wavefunctions in the quantization process of E k T and Q T .

Energy Quantization beyond SI Wavefunctions
To elucidate how SI and NSI wavefunctions cooperate to form the observed quantization levels, we consider the typical quantum motion under a quadratic confining potential V(x) = x 2 /2. The related Schrödinger equation in dimensionless form is: where the total energy E is allowed to be any positive real number. A general solution to the Schrödinger Equation (23), which takes into account NSI wavefunctions, can be expressed in terms of the Whittaker function W(k, m, z) as: where F(α, β, z) is the hypergeometric function, and Γ(α) is the Gamma function. Detailed discussions on the above-mentioned special functions can be found in standard textbooks of physical mathematics [30]. For a given energy E, the obtained solution ψ E (x) is generally NSI, except for the energy eigenvalues E n = n + 1/2, n = 0, 1, 2, · · · , at which Equation (24b) becomes: Depending on whether n is even or odd, simplification of ψ n (x) is given respectively by: where we note 1/Γ(−m) = 0 in Equation (25) for negative integer −m. Combining the above two equations yields the eigenfunctions ψ n (x) = C 1 e −x 2 /2 H n (x) for the quantum harmonic oscillator. The eigenfunctions ψ n (x) are the only solutions to the Schrödinger Equation (23), satisfying the boundary condition (9) and the SI condition (10).
All the existing discussions on energy quantization in the harmonic oscillator focus on the SI eigenfunctions and their linear combinations. Here we are interested in the energy quantization related to the NSI wavefunctions described by Equation (24) with E = n + 1/2. According to Equations (17) and (21), the quantization of E k T and E k + Q T is determined by the numbers of zero and pole of ψ E (x) and ψ E (x). Examining the expression for ψ E (x) given by Equation (24b), we find that ψ E (x) and ψ E (x) do not have any pole over the entire complex plane, because the hypergeometric function F(α, β, z) and its derivative are analytic functions for any z ∈ C. Accordingly, we have P ψ = P ψ = 0, and: Hence the three quantum numbers, n ψ , n p and n ψ , can be determined by the two independent integers: Z ψ and Z ψ , the numbers of zero of ψ E (x) and ψ E (x), respectively. Regarding the computation of Z ψ , we can find the zero of ψ E by solving the roots of the Whittaker function according to Equation (24a): For a given energy E, the resulting root is denoted by x s (E), and Z ψ is the number of x s (E) satisfying Equation (28). The blue line in Figure 1 illustrates the variation of Z ψ with respect to the energy E. Similarly, Z ψ can be found by solving the roots of ψ E (x) = 0: Entropy 2018, 20, 327 9 of 22 The resulting root is denoted by ( ) and the number of ( ) satisfying Equation (29) for a given energy gives the value of . The red line in Figure 1 illustrates the variation of with respect to the energy .  The resulting root is denoted by x eq (E) and the number of x eq (E) satisfying Equation (29) for a given energy E gives the value of Z ψ . The red line in Figure 1 illustrates the variation of Z ψ with respect to the energy E.
As can be seen from Figure 1, when the total energy E increases monotonically, Z ψ and Z ψ exhibit a stair-like distribution in the form of: and Z ψ = 1, Z ψ = 0, as 0 < E ≤ 1/2. Based on the above distributions of Z ψ and Z ψ , the quantization laws derived in Equations (17), (19) and (21) now become: when the total energy E falls in the interval n − 1/2 < E ≤ n + 1/2. All the energies in Equation (31) have been expressed in terms of the multiples of ω. Consequently, as we increase the total energy E monotonically, E k T and E k + Q T increase in a stair-like manner with the step levels given by Equation (31), as shown in Figure 2. Up to this stage, the two components E k and Q in the energy conservation law (7) have been found to be quantized, while the third component, i.e., the externally applied potential V(x), is not a quantized quantity, which otherwise changes continuously with E via the relation: conservation law (7) have been found to be quantized, while the third component, i.e., the externally applied potential ( ), is not a quantized quantity, which otherwise changes continuously with via the relation: . Figure 2. The step changes of 〈 〉 and 〈 + 〉 occur at the SI wavefunctions with = + 1/2, as the total energy changes continuously in a harmonic oscillator. The flat parts of 〈 〉 and 〈 + 〉 are constituted by the NSI with ≠ + 1/2.
The most noticeable point is that the step change of 〈 〉 and 〈 + 〉 occurs at the specific energies = + 1/2, which coincide with the energy eigenvalues of the harmonic oscillator. In other words, the role of the SI condition amounts to determining the discrete energy at which the numbers of zero of ( ) and ( ) exhibit a step jump, while the role of the NSI wavefunctions ( ) with ≠ is to form the flat parts of the stair-like distribution as shown in Figure 1 and Figure 2, where the numbers of zero of ( ) and ( ), or equivalently the time-average energies 〈 〉 and 〈 + 〉 , keep unchanged.

Quantum Bifurcation beyond SI Wavefunctions
As the total energy increases, the wavefunction ( ) transits repeatedly from NSI states to a SI state, once coincides with an energy eigenvalue . In this section, we will show that the encounter with an energy eigenvalue not only causes a step jump of 〈 〉 and 〈 + 〉 , but also causes a nonlinear phenomenon -quantum bifurcation, where the number of equilibrium points of the quantum dynamics experiences an instantaneous change.
With ( ) given by Equation (24a), the quantum dynamics (11) assumes the following dimensionless form: The step changes of E k T and E k + Q T occur at the SI wavefunctions ψ E with E = n + 1/2, as the total energy E changes continuously in a harmonic oscillator. The flat parts of E k T and E k + Q T are constituted by the NSI ψ E with E = n + 1/2.
The most noticeable point is that the step change of E k T and E k + Q T occurs at the specific energies E n = n + 1/2, which coincide with the energy eigenvalues of the harmonic oscillator. In other words, the role of the SI condition amounts to determining the discrete energy E n at which the numbers of zero of ψ E (x) and ψ E (x) exhibit a step jump, while the role of the NSI wavefunctions ψ E (x) with E = E n is to form the flat parts of the stair-like distribution as shown in Figures 1 and 2, where the numbers of zero of ψ E (x) and ψ E (x), or equivalently the time-average energies E k T and E k + Q T , keep unchanged.

Quantum Bifurcation beyond SI Wavefunctions
As the total energy E increases, the wavefunction ψ E (x) transits repeatedly from NSI states to a SI state, once E coincides with an energy eigenvalue E n . In this section, we will show that the encounter with an energy eigenvalue not only causes a step jump of E k T and E k + Q T , but also causes a nonlinear phenomenon -quantum bifurcation, where the number of equilibrium points of the quantum dynamics experiences an instantaneous change.
With ψ E (x) given by Equation (24a), the quantum dynamics (11) assumes the following dimensionless form: where the total energy E is treated as a free parameter, whose critical values for the occurrence of bifurcation are to be identified. The quantum trajectories x(t) solved from Equation (33) provide us with a quantitative comparison between SI and NSI wavefunctions, which otherwise cannot be compared under the probability interpretation of ψ E (x). As can be seen from Equation (33), the equilibrium point x eq of the quantum dynamics is equal to the zero of ψ E (x), while the singular point x s is just the zero of ψ E (x). Hence the step changes of Z ψ and Z ψ shown in Figure 1 also imply the step changes of the numbers of the equilibrium points x eq and the singular points x s , respectively. In other words, we can say that the following two processes occur synchronously as E increases monotonically: one process is the quantization of E k T and E k + Q T regarding the step changes of Z ψ and Z ψ as discussed previously, and the other is the bifurcation of the quantum dynamics (33) regarding the step changes of the equilibrium points and singular points, as to be discussed below.
(1) SI wavefunctions: Firstly, we consider the special cases that the total energy E happens to be one of the eigen energies: E 0 = 1/2, E 1 = 3/2, and E 2 = 5/2. The related eigenfunctions and the eigen-dynamics derived from Equation (33) are given by: These three equations describe the velocity fields and their solutions give the eigen-trajectories for the first three SI states of the harmonic oscillator. It can be shown that the equilibrium points of Equation (34) are centers, while their singular points are saddles. For instance, . x = ix in Equation (34a) has an equilibrium point at the origin with solution given by x(t) = ce it , whose trajectories on the complex plane are concentric circles around the equilibrium point, showing that x = 0 is a center. To show singular points in Equation (34) are saddles, we consider the following complex-valued system with a singular point at the origin: .
where g(x) is analytic at x = 0. In a neighborhood of the origin, Equation (35) can be approximated by x = λ/x leads to the equivalent real-valued nonlinear system: Its solution is a set of hyperbolas expressed by x 2 showing that the singular point x = 0 in Equation (35) are saddles. The centers and saddles of Equation (34) generated by the SI wavefunctions ψ E (x) with E 0 = 1/2, E 1 = 3/2, and E 2 = 5/2 are illustrated in Figure 3, which displays the distribution and movement of the centers and saddles of the quantum dynamics (33) on the horizontal x axis, as the total energy E changes continuously along the vertical axis. Detailed discussions on the quantum trajectories of the SI wavefunctions ψ n (x) for a harmonic oscillator were reported in the literature [31,32]. Here our concern is the quantum trajectories of the NSI wavefunctions ψ E (x) with E = n + 1/2. Quantum trajectories in the first three quantization intervals of E will be examined below, from which a global picture of center-saddle bifurcation can be drawn.
(2) NSI wavefunctions ψ E (x) with 0 < E ≤ 1/2: The wavefunction ψ E (x) in this range of energy is NSI, except for E = 1/2. It seems to be a reasonable conjecture that NSI wavefunctions naturally give rise to unbound quantum trajectories; however, this is not the case. Figure 4 illustrates the quantum trajectories solved from Equation (33) for the NSI wavefunctions ψ E (x) with E = 0.1, E = 0.49, and E = 0.51. In spite of being generated by NSI wavefunctions, the resulting quantum trajectories are bound with slight deviations from the eigen-trajectories of E = 1/2, which are concentric circles around the equilibrium point at the origin, as described by Equation (34a). It appears that SI eigenfunctions are not isolated from the neighboring NSI wavefunctions, because their quantum trajectories can be deformed continuously into each other.  2) NSI wavefunctions ( ) with 0 < ≤ 1/2: The wavefunction ( ) in this range of energy is NSI, except for = 1/2. It seems to be a reasonable conjecture that NSI wavefunctions naturally give rise to unbound quantum trajectories; however, this is not the case. Figure 4 illustrates the quantum trajectories solved from Equation (33) for the NSI wavefunctions ( ) with = 0.1, = 0.49, and = 0.51. In spite of being generated by NSI wavefunctions, the resulting quantum trajectories are bound with slight deviations from the eigen-trajectories of = 1/2, which are concentric circles around the equilibrium point at the origin, as described by Equation (34a). It appears that SI eigenfunctions are not isolated from the neighboring NSI wavefunctions, because their quantum trajectories can be deformed continuously into each other.

3) NSI wavefunctions
( ) with 1/2 < ≤ 3/2: According to Figure 1, the number of equilibrium points of the quantum dynamics (33), jumps from one to two as across the energy eigenvalue = 1/2. This bifurcation phenomenon is illustrated in Figure 5. There is only one equilibrium point at the origin in the energy interval 0 < ≤ 0.5, while beyond the bifurcation point = 1/2, two equilibrium points come out from the origin. Particular attention is paid to the quantum trajectories of = 0.51 depicted in Figure 4d. At first glance, it looks like that the quantum trajectories of = 0.51 have a single equilibrium point at the origin. However, the enlargement of Figure 4d near the origin as illustrated in Figure 5a indicates that the single equilibrium point at the origin for = 1/2 splits into two equilibrium points as increases to 0.51. When increases to 3/2, the two equilibrium points (two centers) move further to = ±1, as described by the quantum dynamics (34b) and illustrated in Figure 5d. Coincident with the splitting of the equilibrium point at the bifurcation point = 1/2, a singular point emerges from the origin in the form of a saddle point. The resulting saddle point pattern in the vicinity of the origin is clearly manifested in Figure 5. It turns out that at the bifurcation energy = 1/2, two kinds of bifurcation occur simultaneously: one bifurcation regards the change of the number of equilibrium points from a single center at the origin into a pair of centers moving apart along the positive and negative real axis as increases from 1/2 to 3/2, and the other bifurcation regards the change from a center into a saddle at the origin.   (3) NSI wavefunctions ψ E (x) with 1/2 < E ≤ 3/2: According to Figure 1, the number of equilibrium points Z ψ of the quantum dynamics (33), jumps from one to two as E across the energy eigenvalue E = 1/2. This bifurcation phenomenon is illustrated in Figure 5. There is only one equilibrium point at the origin in the energy interval 0 < E ≤ 0.5, while beyond the bifurcation point E = 1/2, two equilibrium points come out from the origin. Particular attention is paid to the quantum trajectories of E = 0.51 depicted in Figure 4d. At first glance, it looks like that the quantum trajectories of E = 0.51 have a single equilibrium point at the origin. However, the enlargement of Figure 4d near the origin as illustrated in Figure 5a indicates that the single equilibrium point at the origin for E = 1/2 splits into two equilibrium points as E increases to 0.51. When E increases to 3/2, the two equilibrium points (two centers) move further to x eq = ±1, as described by the quantum dynamics (34b) and illustrated in Figure 5d.  Figure 4d near the origin to illustrate the split of the single equilibrium point at the origin into a pair of equilibrium points as increases from 0.5 to 0.51. In this energy interval, there are two equilibrium points and one singular point at the origin, corresponding to the energy levels = 2 and = 1 as shown in Figure 1.

4) NSI wavefunctions ( ) with 3/2 < ≤ 5/2:
At the energy eigenvalue = 3/2, the value of experiences the second step jump, and a new bifurcation is expected to form here. This prediction is confirmed in Figure 6a, where the enlargement of the velocity field near the origin shows that the saddle-point singularity at the origin for = 3/2 now transforms into a center for = 1.51. As increases further to = 1.6 and = 2 , flow circulation around the origin as a center becomes more apparent. Counting the new equilibrium point emerging from the origin and the already existing pair of centers, the number of equilibrium points increases from two to three as across = 3/2, and remains three in the interval 3/2 < ≤ 5/2. Coincident with the emergence of a new center from the origin at = 3/2, the singular saddle point previously residing at the origin now splits into a pair of saddles with their  Figure 4d near the origin to illustrate the split of the single equilibrium point at the origin into a pair of equilibrium points as E increases from 0.5 to 0.51. In this energy interval, there are two equilibrium points and one singular point at the origin, corresponding to the energy levels Z ψ = 2 and Z ψ = 1 as shown in Figure 1.
Coincident with the splitting of the equilibrium point at the bifurcation point E = 1/2, a singular point emerges from the origin in the form of a saddle point. The resulting saddle point pattern in the vicinity of the origin is clearly manifested in Figure 5. It turns out that at the bifurcation energy E = 1/2, two kinds of bifurcation occur simultaneously: one bifurcation regards the change of the number of equilibrium points from a single center at the origin into a pair of centers moving apart along the positive and negative real axis as E increases from 1/2 to 3/2, and the other bifurcation regards the change from a center into a saddle at the origin.
(4) NSI wavefunctions ψ E (x) with 3/2 < E ≤ 5/2: At the energy eigenvalue E = 3/2, the value of Z ψ experiences the second step jump, and a new bifurcation is expected to form here. This prediction is confirmed in Figure 6a, where the enlargement of the velocity field near the origin shows that the saddle-point singularity at the origin for E = 3/2 now transforms into a center for E = 1.51. As E increases further to E = 1.6 and E = 2, flow circulation around the origin as a center becomes more apparent. Counting the new equilibrium point emerging from the origin and the already existing pair of centers, the number of equilibrium points increases from two to three as E across E = 3/2, and remains three in the interval 3/2 < E ≤ 5/2. Coincident with the emergence of a new center from the origin at E = 3/2, the singular saddle point previously residing at the origin now splits into a pair of saddles with their separation increasing with E. The two saddles move to x s = ± √ 1/2 when E increases to 5/2, as described by Equation (34c) and illustrated in Figure 6d. separation increasing with . The two saddles move to = ± 1/2 when increases to 5/2, as described by Equation (34c) and illustrated in Figure 6d.  Part (a,b) are the enlargements of the velocity field near the origin to illustrate the emergence of a new equilibrium point (a center). In this energy interval, there are three equilibrium points and two singular points, corresponding to the energy levels = 3 and = 2 as shown in Figure 1. Part (d) plots the eigen trajectories for = 2.5 to show the three equilibrium points at = 0, ± 5/2 and two singular points at = ± 1/2.

5) Center-Saddle Bifurcation
When we proceed further, the bifurcations of the equilibrium centers and the singular saddles of the quantum dynamics (33) occur alternatively as increases. To gain a global picture of the bifurcation pattern, we solve the equilibrium points ( ) and the singular points ( ) from Equation (29) and Equation (28), respectively, and then plot them as functions of . The resulting plots generate two sequences of pitchfork bifurcation diagram as shown in Figures 7 and 8 for ( ) and ( ), respectively. It can be seen that the bifurcations of ( ) and ( ) occur alternatively at the critical energies = + 1/2 in such a way that the branches of ( ) bifurcate  Figure 6. Velocity fields and quantum trajectories in the energy interval 1.5 < E ≤ 2.5 show the quantum bifurcation that the number of equilibrium points jumps from two to three as energy across E = 1.5. Part (a,b) are the enlargements of the velocity field near the origin to illustrate the emergence of a new equilibrium point (a center). In this energy interval, there are three equilibrium points and two singular points, corresponding to the energy levels Z ψ = 3 and Z ψ = 2 as shown in Figure 1. Part (d) plots the eigen trajectories for E = 2.5 to show the three equilibrium points at x eq = 0, ± √ 5/2 and two singular points at x s = ± √ 1/2.

(5) Center-Saddle Bifurcation
When we proceed further, the bifurcations of the equilibrium centers x eq and the singular saddles x s of the quantum dynamics (33) occur alternatively as E increases. To gain a global picture of the bifurcation pattern, we solve the equilibrium points x eq (E) and the singular points x s (E) from Equations (29) and (28), respectively, and then plot them as functions of E. The resulting plots generate two sequences of pitchfork bifurcation diagram as shown in Figures 7 and 8 for x eq (E) and x s (E), respectively. It can be seen that the bifurcations of x eq (E) and x s (E) occur alternatively at the critical energies E n = n + 1/2 in such a way that the branches of x eq (E) bifurcate sequentially at E = 0 + (1/2), 2 + (1/2), 4 + (1/2), · · · , while the branches of x s (E) bifurcate sequentially at E = 1 + (1/2), 3 + (1/2), 5 + (1/2), · · · .  ( ) with respect to the total energy . The number of ( ) at each , denoted by the blue dots, is equal to the energy level as plotted in Figure 1. The branches of the ( ) curves start sequentially at = 0, 3/2, 7/2, ⋯, and bifurcate sequentially at = 1/2, 5/2, 9/2, ⋯. Except for the bifurcation points (the blue dots), the entire sequential bifurcation diagram is formed by the NSI wavefunctions ( ) with ≠ + 1/2. Figure 8. A sequence of pitchfork bifurcation curves shows the variation of singular points ( ) with respect to the total energy . The number of ( ) at each , denoted by the red dots, is equal to the energy level as plotted in Figure 1. The branches of the ( ) curves start sequentially at = 1/2, 5/2, 9/2, ⋯, and bifurcate sequentially at = 3/2, 7/2, 11/2, ⋯. Except for the bifurcation Figure 7. A sequence of pitchfork bifurcation curves shows the variation of equilibrium points x eq (E) with respect to the total energy E. The number of x eq (E) at each E, denoted by the blue dots, is equal to the energy level Z ψ as plotted in Figure 1. The branches of the x eq (E) curves start sequentially at E = 0, 3/2, 7/2, · · · , and bifurcate sequentially at E = 1/2, 5/2, 9/2, · · · . Except for the bifurcation points (the blue dots), the entire sequential bifurcation diagram is formed by the NSI wavefunctions Furthermore, it is worth noting that except for the bifurcation points (the blue dots in Figure 7 and the red dots in Figure 8), the sequential bifurcation diagram is constructed entirely by the NSI wavefunctions ψ E (x) with E = n + 1/2. Without the participation of the NSI wavefunctions, adjacent eigenfunctions lose their interconnection and a continuous description of the bifurcation sequence becomes impossible. The other perspective of center-saddle bifurcation can be gained from Figure 3, where we can see that centers and saddles appear alternatively at the origin as the total energy E increases monotonically along the vertical axis.  Figure 7. A sequence of pitchfork bifurcation curves shows the variation of equilibrium points ( ) with respect to the total energy . The number of ( ) at each , denoted by the blue dots, is equal to the energy level as plotted in Figure 1. The branches of the ( ) curves start sequentially at = 0, 3/2, 7/2, ⋯, and bifurcate sequentially at = 1/2, 5/2, 9/2, ⋯. Except for the bifurcation points (the blue dots), the entire sequential bifurcation diagram is formed by the NSI wavefunctions ( ) with ≠ + 1/2. Figure 8. A sequence of pitchfork bifurcation curves shows the variation of singular points ( ) with respect to the total energy . The number of ( ) at each , denoted by the red dots, is equal to the energy level as plotted in Figure 1. The branches of the ( ) curves start sequentially at = 1/2, 5/2, 9/2, ⋯, and bifurcate sequentially at = 3/2, 7/2, 11/2, ⋯. Except for the bifurcation points (the red dots), the entire sequential bifurcation diagram is formed by the NSI wavefunctions ( ) with ≠ + 1/2.
Furthermore, it is worth noting that except for the bifurcation points (the blue dots in Figure 7 and the red dots in Figure 8), the sequential bifurcation diagram is constructed entirely by the NSI wavefunctions ( ) with ≠ + 1/2 . Without the participation of the NSI wavefunctions, Figure 8. A sequence of pitchfork bifurcation curves shows the variation of singular points x s (E) with respect to the total energy E. The number of x s (E) at each E, denoted by the red dots, is equal to the energy level Z ψ as plotted in Figure 1. The branches of the x s (E) curves start sequentially at E = 1/2, 5/2, 9/2, · · · , and bifurcate sequentially at E = 3/2, 7/2, 11/2, · · · . Except for the bifurcation points (the red dots), the entire sequential bifurcation diagram is formed by the NSI wavefunctions ψ E (x) with E = n + 1/2.

(6) Synchronicity between quantization and bifurcation
We recall that the number of x eq (E) at each E is just the number of zero of ψ E (x), which gives the quantization level of E k + Q T . This relation indicates that the bifurcation of the equilibrium center point x eq (E) and the quantization of E k + Q T occur synchronously. Similarly, because the number of x s (E) at each E is the number of zero of ψ E (x), which gives the quantization level of E k T , the bifurcation of the singular saddle point x s (E) is thus synchronous with the quantization of E k T . Table 1 lists the numbers of saddles and centers, and the energy levels of E k + Q T and E k T for the first several energy intervals. As can be seen, the numbers of saddles and centers change synchronously with the change of E k + Q T and E k T . It is noted that as energy E varies continuously during the quantization and bifurcation processes, the instantaneous changes of energy levels and equilibrium points are triggered by the SI condition E = n + 1/2.

Spin Degree of Freedom beyond SI Wavefunctions
The role of the Schrödinger equation has long been considered as describing spinless particles only, because the Schrödinger charge current for the s-states of a hydrogen-like atom vanishes and produces no intrinsic angular momentum. However, based on the observation that in the absence of a magnetic field, the Pauli equation reduces to the Schrödinger equation, it was pointed out [33] that the Schrödinger equation must be regarded as describing an electron in an eigenstate of spin and not, as universally supposed, an electron without spin. According to the dBB trajectory approach [34,35], spin is interpreted as a dynamical property of electron motion and is attributed to a circulating movement of a point, i.e., to a pure orbital motion, but not to an extended spinning object. To be consistent with the Dirac theory and with the condition of Lorentz invariance, Holland [34] proposed that the Schrödinger charge current must be supplemented by a spin magnetization current, which is generated by a circulating flow of energy in the wave field of the electron [36,37].
The NSI solutions to the Schrödinger equation considered in the present paper might give an alternative explanation for the origin of particle's spin motion. The Schrödinger Equation (8) with given energy E actually has two independent solutions. It is surprising to find that the quantum trajectories generated by the two independent solutions are indistinguishable, except for their directions of rotation. Inspecting the quantum trajectories shown in Figures 4-6, it appears that all the trajectories, either generated by SI or NSI wavefunctions, rotate counterclockwise (CCW). In fact, all the trajectories produced by the general wavefunctions given by Equation (24) rotate in the same direction, because Equation (24) only gives one of the independent solutions. A complete general solution to the Schrödinger Equation (23) comprises two independent parts: where ψ CCW (x) is the solution considered previously and ψ CW (x) is the other independent solution, whose quantum trajectories rotate clockwise (CW). The wavefunction ψ CW (x) represents the second half of solutions to the Schrödinger Equation (23), which is NSI for any energy E and we usually take the neglect of it as granted.
The consideration of the wavefunctions ψ CW (x) helps to identify the additional degree of freedom independent of the particle's orbital motion. To highlight the difference between ψ CCW (x) and ψ CW (x), their velocity fields computed by Equation (11) with E = 1/2 are illustrated in Figures 9a,b, respectively. The velocity field of ψ CCW (x) is identical to Figure 4c, which shows circular flows surrounding the origin counterclockwise. By contrast, the velocity field of ψ CW (x) depicted in Figure 9b appears to be clockwise circulation around the origin. The quantum trajectories generated by ψ CW (x) are almost indistinguishable from those generated by ψ CCW (x), and the only difference between them is the directions of rotation. In addition to the orbital motion, the new degree of freedom manifested in the combination of ψ CCW (x) and ψ CW (x) is the dual directions of rotation, which is otherwise invisible along a single trajectory generated by either ψ CCW (x) or ψ CW (x). Due to their same spatial motion with dual directions, ψ CCW (x) and ψ CW (x) can be reasonably recognized as the same spatial solution to the Schrödinger equation but with opposing directions of spin. The evaluations of ( ) at = and = for = 1/2 are plotted in Figure 9d. Despite of the adverse nature between SI wavefunction ( ) and NSI wavefunction ( ), their total potentials exhibit a high degree of resemblance, which explains the observed similarity in the velocity fields of Figure 9a,b. The comparisons regarding the magnitudes of and the total potential ( ) for = 3/2 and = 5/2 are shown in Figure 10, where we observe that except for the region neighboring the origin, ( ) is close to ( ) and they become identical as | | → ∞. This asymptotic identity leads to the dual velocity fields derived in Equation (40). The reason underlying the opposite rotation of ψ CCW (x) and ψ CW (x) can be explained by using the asymptotic expansion property for the Whittaker function: With this property, the asymptotic expansions of ψ CCW (x) and ψ CW (x) take the following forms: The asymptotic quantum dynamics of ψ CCW (x) and ψ CW (x) then can be derived as: Irrespective to the energy E, both of the asymptotic quantum dynamics converge to the ground-state quantum dynamics (34a) with the only difference in their rotation direction. This similarity in the velocity field of ψ CCW (x) and ψ CW (x) has been ignored in the literature. Probability interpretation is only concerned with the square-integrable condition of ψ CCW (x) and ψ CW (x), as shown in Figure 9c, which obviously fails to explain the similarity in the quantum dynamics of ψ CCW (x) and ψ CW (x).
The factor dominating the similarity in the quantum dynamics of ψ CCW (x) and ψ CW (x) can be explained by Equation (5b): The total potential V T = V + Q, comprising the applied potential V and the quantum potential Q, dominates the time evolution of the QMF p. For the case of a harmonic oscillator, V T turns out to be (in dimensionless form): The evaluations of V T (ψ) at ψ = ψ CCW and ψ = ψ CW for E = 1/2 are plotted in Figure 9d. Despite of the adverse nature between SI wavefunction ψ CCW (x) and NSI wavefunction ψ CW (x), their total potentials exhibit a high degree of resemblance, which explains the observed similarity in the velocity fields of Figure 9a,b. The comparisons regarding the magnitudes of ψ and the total potential V T (ψ) for E = 3/2 and E = 5/2 are shown in Figure 10, where we observe that except for the region neighboring the origin, V T (ψ CCW ) is close to V T (ψ CW ) and they become identical as |x| → ∞ . This asymptotic identity leads to the dual velocity fields derived in Equation (40).  Schrödinger equation is a second-order differential equation with respect to its spatial coordinates, and its complete solution should be composed of two independent solutions. The common belief that Schrödinger equation is unable to describe spin motion seems to stem from our disregard of one of the independent solution. While probability interpretation of wavefunctions excludes the NSI solution ( ) from the general solution ( ), the spin degree of freedom is removed at the same time. Under the framework of quantum H-J formalism, we have seen that by incorporating ( ) with ( ) to form a general wavefunction as expressed by Equation (37), both spatial and spin motion can be described by the Schrödinger equation, and even for onedimensional quantum motion, the spin degree of freedom can be manifested as the dual rotations on the complex plane. Figure 10. The comparisons between CCW solution ψ CCW (x) and CW solution ψ CW (x) regarding the magnitudes of ψ and the total potential V T (ψ) for E = 3/2 and E = 5/2, respectively. Schrödinger equation is a second-order differential equation with respect to its spatial coordinates, and its complete solution should be composed of two independent solutions. The common belief that Schrödinger equation is unable to describe spin motion seems to stem from our disregard of one of the independent solution. While probability interpretation of wavefunctions excludes the NSI solution ψ CW (x) from the general solution ψ E (x), the spin degree of freedom is removed at the same time. Under the framework of quantum H-J formalism, we have seen that by incorporating ψ CW (x) with ψ CCW (x) to form a general wavefunction as expressed by Equation (37), both spatial and spin motion can be described by the Schrödinger equation, and even for one-dimensional quantum motion, the spin degree of freedom can be manifested as the dual rotations on the complex x plane.

Conclusions
Standard approach to energy quantization in confined quantum systems is to seek for the allowable energies E n for which the time-independent Schrödinger equation has square-integrable solutions. The obtained energy eigenvalue E n is recognized as the quantization level of the system's total energy. In this paper we have given a renewed interpretation for E n and considered NSI solutions ψ E (x) to the Schrödinger equation by releasing the SI requirement. The release of this requirement leads to several new findings as summarized in the following points: • Universal quantization laws: The total energy E = E k + Q(x) + V(x) derived from the time-independent Schrödinger equation is shown to be conserved, but not quantized. Regardless of the confining potential V(x), quantization always occurs in the kinetic energy E k T and the quantum potential Q T , whose values can only change by an integer multiple of ω/2. • Renewed meaning of the energy eigenvalues: The energy eigenvalues E n derived conventionally from the SI condition are shown to be the special energies at which the quantization levels of E k T and E k + Q T experience a step jump.

•
The origin of energy quantization: Energy quantization in a confined system originates from the discrete change of the numbers of zero of ψ E (x) and ψ E (x), whose values determine the quantization levels of E k T and E k + Q T .

•
Concurrence of Quantization and bifurcation: Bifurcations of equilibrium center points and singular saddle points of the quantum dynamics are shown to be synchronous, respectively, with the quantization process of E k T and E k + Q T , as the total energy E increases monotonically. • Undivided SI and NSI wavefunctions: probability interpretation isolates SI wavefunctions from NSI wavefunctions; however, under the quantum H-J formalism, SI and NSI wavefunctions are indivisible with continuously connected velocity field and quantum trajectories.

•
The role of NSI wavefunctions in energy quantization: Both SI and NSI wavefunctions contribute to the energy quantization. SI wavefunctions help to locate the bifurcation points at which E k T and E k + Q T have a step jump, while NSI wavefunctions form the flat parts of the stair-like distribution of the quantized energies.

•
The role of NSI wavefunctions in spin: The second-order Schrödinger equation generally contains two independent solutions with opposite rotation on the complex plane. The inclusion of both solutions allows the Schrödinger equation to describe the spatial motion as well as the spin motion.
At the present stage of this research, we can only say that if the SI requirement is released tentatively, above information can be gained from NSI wavefunctions. We need further experimental results to support the existence of NSI wavefunctions in confined quantum systems. The result of this paper gives a clue to the design of such an experiment. Based on the SI condition, what we consider to be quantized is the particle's total energy. On the other hand, if the SI condition is released, what to be quantized is the particle's time-average kinetic energy. Experiments on particle's motion in confining potentials can be performed to verify which prediction is correct.