Fourth Painlevé and Ermakov equations: quantum invariants and new exactly-solvable time-dependent Hamiltonians

During the last several decades, physicists have realized the importance of nonlinear equations in the study of physical systems, even in cases where linear equations govern the dynamical laws. For instance, in quantum mechanics, the nonlinear Riccati equation [1, 2] has played a fundamental role in the construction and study of new exactly-solvable models, for it relates the factorization method [3, 4, 5, 6, 8, 7], with the Darboux transformation [9]. The latter can be formulated in the so-called supersymmetric quantum mechanics (SUSYQM) [13, 12, 10, 11]. In this regard, systems with spectrum on-demand are obtained either by adding energy level not present in the original model, or by removing levels through the Darboux-Crum transformation [14]. This formalism represents an outstanding progress in the study of quantum systems; it has allowed extending the families of exactly-solvable models  [15, 16] beyond the conventional harmonic oscillator, hydrogen atom, the interaction between diatomic molecules, and some few others. Moreover, quantum mechanics in the non-Hermitian regime has been explored by generalizing the factorization method and allowing the Riccati equation to be a complex-valued function. In this form, the realitiy of the spectrum is preserved [17, 18, 19] in systems with either broken and unbroken $PT$-symmetry, extending the conventional systems with $PT$-symmetry and real spectrum [20, 21, 22]. From the latter, the nonlinear Ermakov equation [23, 24, 25] emerges naturally from the complexified Riccati one, the solutions of which ensures the regularity of the new complex-valued potentials. For more details on the applications of the Riccati and Ermakov equations in physics, see [2].

Among other nonlinear equations, we have the Painlevé transcendentals, a family of six nonlinear equations $P_{\rm I}$-$P_{\rm VI}$ with complex parameters, whose solutions are in general transcendental, that is, theyr are not be expressed in terms of classical functions [26, 27]. Nevertheless, for some specific values of the parameters, a seed function can be used to generate a complete hierarchy of solutions through the Bäcklund transformation [28], which can be thought as a nonlinear counterpart of the recurrence relations. In particular, the fourth Painlevé equation can be taken into a Riccati equation with the appropriate choice of the parameters. we thus solve a “simpler” nonlinear equation instead. Also, the fourth Painlevé equation has also brought new result in the trend of orthogonal polynomials, where new families were discovered though the hierarchies of rational solutions in terms of the generalized Okamoto, generalized Hermite and Yablonskii-Vorob’ev polynomials [29]. The Painlevé transcendental have also found interesting applications in the study of physical models in nonlinear optics [30], quantum gravity [31], and SUSYQM [33, 34, 32], to mention some. Interestingly, the fourth Painlevé equation arises quite naturally in third-order shape-invariant SUSYQM [33], where the parameters of the Painlevé transcendental define the eigenvalues of the new Hamiltonians. The respective intertwining operators serve at the same time as ladder operators, from where the eigenfunctions are determined. A striking feature of this approach is that, in general, the so-constructed intertwining operators are not in general factorizale in terms of first-order operators. Thus, the results obtained in this way generalize those of [35]. It is worth to notice that higher-order ladder operators have been also studied, in a different way, in the context of supersymmetric (SUSY) partners for the stationary oscillator in both the Hermitian [16, 36] and non-Hermitian regimes [18].

Although a vast literature on families of solvable stationary systems is available, the time-dependent counterparts have not been widely explored. The difficulty lies in the dynamical law, the Schrödinger equation, which is defined in terms of a partial differential equation that, in general, can not be reduced to an ordinary differential equation. Under some circumstances, we can extract information of the system through approximation techniques such as the sudden and the adiabatic approximations [37]. The latter restricts the range of applicability of the so-obtained solutions. Despite all these difficulties, time-dependent phenomena find exciting applications in physical systems such as electromagnetic traps of charged particles [38, 41, 39, 40], plasma physics [42], and in optical-analogs under the paraxial approximation [43, 44, 45].

In contradistinction to the stationary cases, the lack of an eigenvalue equation in time-dependent systems prevents us from implementing the Darboux-transformation directly, and some workarounds are in order. The latter has been addressed by Bargov-Samsonov [46, 47], where some intertwining operators allow us to relate an exactly solvable Schrödinger equation with another unknown one. In analogy to the stationary Darboux transformation, the solutions of the new models are inherited from the former one. Let us mention that orthogonality is no longer a property that can be taken for granted [48], since the method by itself does not provide essential information about the system such as the constants of motion, which have to be determined separately. Despite such a difficulty, several new families of exactly-solvable time-dependent potentials have been reported in the literature [48, 49, 50, 51].

Among nonstationary quantum systems, the parametric oscillator [52, 53, 54, 55] is perhaps the most well-known model that admits a set of exact solutions. Lewis and Riesenfeld [52] addressed the problem by noticing the existence of a nonstationary eigenvalue equation associated with the appropriate constant of motion (quantum invariant) of the system in which the time dependence appears in the coefficients of the related ordinary differential equation. The latter eigenvalue equation can indeed be factorized in such a way that the Darboux transformation¹¹1An eigenvalue equation for the time-dependent Hamiltonian is still attainable in the context of the adiabatic approximation [37]. is applied with ease [56, 57], resulting in a new quantum invariant rather than a Hamiltonian. Then, the appropriate ansatz allows to determine the respective Hamiltonian and time-dependent potentials with ease [56]. The solutions, and the complex-phases introduced by Lewis-Riesenfeld, are inherited from the former system, ensuring an orthogonal set of solutions for the new system.

In this work, we combine the solutions of the Ermakov and fourth Painlevé equation to address the construction of new time-dependent Hamiltonians. This is achieved by considering third-order intertwining operators in the spatial variable with time-dependent coefficients. Those operators generate a third-order shape-invariant condition with respect to an unknown quantum invariant, which is introduced as a deformation of the one associated with the parametric oscillator. In this form, by working with a quantum invariant rather than a Hamiltonian, we generalize the construction presented in [33, 34], and the time dependence is introduced into the intertwining operators through the solutions of the Ermakov equation, which ensure the regularity of the resulting quantum invariant and its eigenfunctions at each time. On the other hand, with aid of the appropriate reparametrization, the fourth Painlevé equation is achieved, the parameters and solutions of which determine the spectral information and the exact form of the quantum invariant. Then, we modify the transitionless tracking algorithm [58] to construct the time-dependent Hamiltonians from the quantum invariant, from where the respective solutions of the Schrödinger equation are determined with the addition of a time-dependent complex-phase.

The text is structured as follows. In Sec. 2, we introduce the basic notions of shape-invariance for time-dependent Hamiltonians and their respective quantum invariants. Then, a couple of differential ladder operators of third-order are introduced such that an initial, and unknown, quantum invariant satisfies a higher-order shape-invariant relationship. From the latter, the explicit form of the ladder operators and the quantum invariants are determined. In Sec. 3, with the aid of the third-order ladder operators, we determine the respective spectral information of the quantum invariant. In Sec. 4, the time-dependent Hamiltonians associated with the quantum invariants are identified, together with the solutions to the Schrödinger equation. In Sec. 5, we discuss the solutions of the Ermakov equation for some specific time-dependent frequency profiles. In particular, it is shown that the constant-frequency case leads to periodic potentials, whose solutions are in agreement with the Floquet theorem. Also, the appropriate limit to recover the well-known stationary results is presented. In turn, in Sec. 6, we consider some particular solutions of the Painlevé equation obtained through solutions of the Riccati equation, in the form of rational solutions, or by solutions of another nonlinear models such as the derivative nonlinear Schrödinger equation. For completeness, in App. A, we briefly revisit the parametric oscillator and its solutions through the approach of Lewis-Riesenfeld.

In quantum mechanics, the quantum invariants play a fundamental role in determining the exact solutions of the quantum models. For stationary systems (time-independent Hamiltonians), it is straightforward to realize that the Hamiltonian is a constant of motion, such that it leads to an eigenvalue equation in the form of a Sturm-Liouville problem. For time-dependent Hamiltonians, the determination of such quantum invariants becomes a challenging task in most of the cases. A prime example is given by the parametric oscillator, where the respective quantum invariants are determined with relative ease. As pointed out in [52], such an invariant admits a nonstationary eigenvalue equation, where its spectral information leads to the solutions of the Schrödinger equation. For the sake of self-consistency, in App. A we provide a brief discussion on the matter.

Thoughout this manuscript, we construct new time-dependent models and their respective solutions using a new approach based on the ladder operator structure associated with the quantum invariant. The method relies on the existence of an unknown quantum invariant $\hat{I}_{1}(t)$, which admits a set of ladder operators $\{\hat{A}(t),\hat{A}^{\dagger}(t)\}$ defined in coordinate representation as $N$-order differential operators in the spatial variable with time-dependent coefficients. In particular, to reduce the possible family of quantum invariants, we consider $\hat{I}_{1}(t)$ as a deformation of the invariant associated with the parametric oscillator, $\hat{I}_{0}(t)$. Thus, for a fixed order $N$, we determine the exact form of $\hat{I}_{1}(t)$. It is worth to mention that ladder operators of first and second-order have been reported for the parametric oscillator [54, 56] and the nonstationary singular oscillator [57, 59]. Nevertheless, in those cases, the quantum invariants are already known, and the ladder operators are constructed following the polynomial structure of the eigenfunctions. Throughout the rest of the manuscript, we focus in third-order differential ladder operators and the related quantum invariants. In this case, we determine several families of time-dependent Hamiltonians $\hat{H}_{1}(t)$, together with the respective solutions to the Schrödinger equation.

Before proceeding, we would like to stress the meaning of shape-invariance in supersymmetric quantum mehcnics (SUSYQM) for both stationary and nonstationary systems. It is said that two time-independent Hamiltonians $\hat{H}_{\pm}$ are shape-invariant [13] if their respective potentials $V_{+}(x;\{c_{n}\})$ and $V_{-}(x;\{d_{n}\})$, with $\{c_{n}\}$ and $\{d_{n}\}$ two sets of constant parameters, are related by the condition $V_{+}(x;\{c_{n}\})=V_{-}(\{d_{n}\})+S(\{c_{n}\})$, where $S(\{c_{n}\})$ is a function of the set of parameters $\{c_{n}\}$ and independent of $x$. In turn, for nonstationary systems, we define the shape-invariance considering two quantum invariants $\hat{I}_{\pm}(t)$ by means of the relationship $\hat{I}_{+}(t)=\hat{I}_{-}(t)+f_{0}$, with $f_{0}$ a real constant. The latter implies that the respective time-dependent Hamiltomnians $\hat{H}_{\pm}(t)$ are related as $\hat{H}_{+}(t)=\hat{H}_{-}(t)+f(t)$, with $f(t)$ an arbitrary real-valued function of time. As a consequence, the solutions of the respective Schrödinger equations differ only from a global time-dependent complex-phase $e^{i\int^{t}dt^{\prime}f(t^{\prime})}$, see [56] for details.

From the previous considerations, let us introduce two unknown quantum invariants $\hat{I}_{1,2}(t)$, where $\hat{I}_{1}(t)$ is the invariant operator under consideration, and $\hat{I}_{2}(t)$ serves as an auxiliary operators to determined $\hat{I}_{1}(t)$. From the latter, there are two time-dependent Hamiltonians $\hat{H}_{1,2}(t)$ such that

$$\frac{d\hat{I}_{j}(t)}{dt}=i\left[\hat{H}_{j}(t),\hat{I}_{j}(t)\right]+\frac{\partial\hat{I}_{j}(t)}{\partial t}=0\,,\quad j=1,2\,.$$

(1)

We focus on self-adjoint time-dependent Hamiltonians and quantum invariants. It is thus guaranteed that the nonstationary eigenvalue equations

$$\hat{I}_{j}(t)\phi^{(j)}_{n}(x,t)=\Lambda^{(j)}_{n}\phi^{(j)}_{n}(x,t)\,,\quad j=1,2,$$

(2)

lead to real and time-independent eigenvalues $\Lambda^{(j)}_{n}$, and orthogonal nonstationary eigenfunctions $\phi^{(j)}_{n}(x,t)$ that satisfy the finite-norm condition $|\langle\phi^{(j)}_{n}(t)|\phi^{(j)}_{n}(t)\rangle|<\infty$. The orthogonality condition can be used as a base to construct the vector spaces $\mathcal{H}_{j}(t)=Span\{\phi^{(j)}_{n}(x,t)\}_{n=0}^{\infty}$, with $j=1,2$.

As mentioned above, we construct the quantum invariant $\hat{I}_{1}(t)$ from the third-order SUSYQM shape-invariant condition, defined in term of the intertwining relationships

$$\hat{I}_{1}(t)\hat{A}^{\dagger}(t)=\hat{A}^{\dagger}(t)\left[\hat{I}_{1}(t)+2\lambda\right]\,,\quad\hat{I}_{1}(t)\hat{A}(t)=\hat{A}(t)\left[\hat{I}_{1}(t)-2\lambda\right]\,,$$

(3)

where $\hat{A}(t)$ and $\hat{A}^{\dagger}(t)$ are also the annihilation and creation operators, respectively, for $\phi_{n}^{(1)}(x,t)$. Indeed, the action of $\hat{A}^{\dagger}(t)$ ($\hat{A}(t)$) on $\phi_{n}^{(1)}(x,t)$ increases (decreases) the eigenvalue $\Lambda_{n}^{(1)}$ by $2\lambda$ units.

Now, to provide a specific form to the quantum invariants $\hat{I}_{1,2}(t)$, we consider them as deformations of the quantum invariant associated with the parametric oscillator $\hat{I}_{0}(t)$, see Eq. (A-5), that is,

$$\hat{I}_{j}(t)=\hat{I}_{0}(t)+R_{j}(\hat{x},t)\equiv-\sigma^{2}\frac{\partial^{2}}{\partial x^{2}}+ix\sigma\dot{\sigma}\frac{\partial}{\partial x}+R(x,t)+R_{j}(x,t)\,,\quad j=1,2,$$

(4)

where $\sigma\equiv\sigma(t)$ given in (A-7) is solution to the Ermakov equation associated with the parametric oscillator (A-6), and $\dot{\sigma}\equiv d\sigma/dt$. The function $R(x,t)$ is given as

$$R(x,t)=\left(\frac{\dot{\sigma}^{2}}{4}+\frac{1}{\sigma^{2}}\right)x^{2}+i\frac{\dot{\sigma}\sigma}{2}\,,$$

(5)

and $R_{1,2}(x,t)$ are real-valued functions, where $R_{1}(x,t)$ is determined from the shape-invariant condition (3), and $R_{2}(x,t)$ defines the auxiliary invariant $\hat{I}_{2}(t)$. As discussed in [33], higher-order intertwining operators can be factorized as products of first-order operators for a reducible factorization, or as combination of first and second-order operators for irreducible factorizations [60]. In the sequel, we focus on the reducible factorization of the intertwining relationships (3); however, for the sake of completeness, we discuss the more general irreducible factorization in this section. We thus decompose the set of intertwining operators $\{\hat{A}(t),\hat{A}^{\dagger}(t)\}$ as the product of first and second-order differential operators $\{\hat{Q}^{\dagger}(t),\hat{Q}(t)\}$ and $\{\hat{M}^{\dagger}(t),\hat{M}(t)\}$, respectively, as follows:

$$\hat{A}^{\dagger}(t)=\hat{Q}^{\dagger}(t)\hat{M}(t)\,,\quad\hat{A}=\hat{M}^{\dagger}(t)\hat{Q}(t)\,.$$

(6)

The new operators give rise the additional set of intertwining relationships of the form

$$\hat{I}_{1}(t)\hat{Q}^{\dagger}(t)=\hat{Q}^{\dagger}(t)\left[\hat{I}_{2}(t)+2\lambda\right]\,,\quad\hat{I}_{2}(t)\hat{Q}(t)=\hat{Q}(t)\left[\hat{I}_{1}(t)-2\lambda\right]\,,$$

(7)

$$\hat{I}_{2}(t)\hat{M}(t)=\hat{M}(t)\hat{I}_{1}(t)\,,\quad\hat{I}_{1}(t)\hat{M}^{\dagger}(t)=\hat{M}^{\dagger}(t)\hat{I}_{2}(t)\,,$$

(8)

where in the latter we have introduced the auxiliary quantum invariant $\hat{I}_{2}(t)$ as an intermediate. In contradistinction to (3), the Eqs. (7)-(8) by themselves do not define a shape-invariant relation. Nevertheless, their combined action take us back to the shape-invariant condition (3), see Fig. 1 for details.

Now, given that $\hat{Q}(t)$ and $\hat{Q}^{\dagger}(t)$ are considered as first-order differential operators, we use the general form introduced in [57], that is,

$$\hat{Q}^{\dagger}=\sigma\frac{\partial}{\partial x}+\mathfrak{w}(x,t)\,,\quad\hat{Q}=-\sigma\frac{\partial}{\partial x}+\mathfrak{w}^{*}(x,t)\,,$$

(9)

with $\mathfrak{w}(x,t)$ a complex-valued function, $\sigma\equiv\sigma(t)$ given in (A-7), and $f^{*}$ stands for the complex-conjugate of $f$. In turn, the second-order differential operators $\hat{M}(t)$ and $\hat{M}^{\dagger}(t)$ are constructed as a generalization of those reported in [60] by introducing the time-dependent coefficients of the form

		$\displaystyle\hat{M}^{\dagger}=\sigma^{2}\frac{\partial^{2}}{\partial x^{2}}-2g(x,t)\frac{\partial}{\partial x}+b(x,t)\,,$		(10)
		$\displaystyle\hat{M}=\sigma^{2}\frac{\partial^{2}}{\partial x^{2}}+2g^{*}(x,t)\frac{\partial}{\partial x}+b^{*}(x,t)-2\left[g^{\prime}(x,t)\right]^{*}\,,$

where $b(x,t)$ and $g(x,t)$ are complex-valued functions, and $g^{\prime}(x,t)$ stands for the partial derivative of $g(x,t)$ with respect to $x$.

The complex-valued functions $\mathfrak{w}(x,t)$, $g(x,t)$ and $b(x,t)$ are determined from the respective intertwining relationships. For instance, $\mathfrak{w}(x,t)$ is obtained after substituting (9) in (7), leading to

$$\mathfrak{w}(x,t)=-i\frac{\dot{\sigma}}{2}x+W(z(x,t))\,,\quad z(x,t):=\frac{x}{\sigma}\,.$$

(11)

Given that the solution to the Ermakov equation $\sigma(t)$ is a nodeless function, we can guarantee that the reparametrizated variable $z(x,t)$ is non-singular for $t\in\mathbb{R}$. In turn, $W(z(x,t))$ is a real-valued function that solves the Riccati equations

$$z^{2}+R_{1}(z)=\partial_{z}W+W^{2}\,,\quad z^{2}+R_{2}(z)=-\partial_{z}W+W^{2}-2\lambda\,,\quad\partial_{z}\equiv\frac{\partial}{\partial z}\,.$$

(12)

From (12) we also get $R_{2}(z)-R_{1}(z)=-2\partial_{z}W-2\lambda$, which resembles the conventional relationship between the potential and the super-potential of the conventional stationary SUSYQM construction [33, 34].

On the other hand, the functions $g(x,t)$ and $b(x,t)$ are computed after inserting (10) in (8), leading to

		$\displaystyle g(x,t)=i\frac{\sigma\dot{\sigma}}{2}x+\sigma G(z(x,t))\,,$		(13)
		$\displaystyle b(x,t)=i\dot{\sigma}xG(z(x,t))-i\frac{\dot{\sigma}\sigma}{2}-\frac{\dot{\sigma}^{2}}{4}x^{2}+B(z(x,t))\,,$

with the real-valued functions $G(z(x,t))$ and $B(z(x,t))$ determined from the nonlinear relationships

		$\displaystyle B=2G^{2}+\partial_{z}G-\left(z^{2}+R_{2}(z)\right)+\gamma\,,$		(14)
		$\displaystyle z^{2}+R_{1}(z)=-2\partial_{z}G+G^{2}+\frac{\partial^{2}_{z}G}{2G}-\frac{(\partial_{z}G)^{2}}{4G^{2}}-\frac{d}{4G^{2}}+\gamma\,,$
		$\displaystyle z^{2}+R_{2}(z)=2\partial_{z}G+G^{2}+\frac{\partial^{2}_{z}G}{2G}-\frac{(\partial_{z}G)^{2}}{4G^{2}}-\frac{d}{4G^{2}}+\gamma\,,$

with $\gamma$ and $d$ constants of integration with respect to $z(x,t)$, that is, those constants do not depend on $x$ or $t$. From (14) we obtain a complementary relationship of the form $R_{2}(z)-R_{1}(z)=4\partial_{z}G$ that, together with the one obtained from (12), gives

$$W(z)=-2G(z)-\lambda z\,.$$

(15)

A differential equation for $G(z)$ can be found after substituting (15) into any of the Riccati equations in (12) and comparing with  (14). The straightforward calculation leads to

$$\partial^{2}_{z}G=\frac{(\partial_{z}G)^{2}}{2G}+6G^{3}+8\lambda zG^{2}+2\left[\lambda^{2}z^{2}-(\gamma+\lambda)\right]G+\frac{d}{2G}\,,$$

(16)

where the following reparametrizations:

$$y=\sqrt{\lambda}z\,,\quad G=\frac{\sqrt{\lambda}}{2}w(y)\,,\quad\alpha=\frac{\gamma}{\lambda}+1\,,\quad\beta=\frac{2d}{\lambda^{2}}\,,$$

(17)

allow us to rewrite (16) as the fourth Painlevé differential equation [61]

$$w_{yy}=\frac{(w_{y})^{2}}{2w}+\frac{3}{2}w^{3}+4yw^{2}+2(y^{2}-\alpha)w+\frac{\beta}{w}\,.$$

(18)

Solutions for the fourth Painlevé equation have been extensively studied in the literature, in particular it is known that $w(y)$ can be determined in terms of elementary functions [26, 27].

Before finishing this section, we would like to recall that the factorization (6) allowed us to find the functions $R_{1,2}(x,t)$, which define uniquely the respective quantum invariants $\hat{I}_{1,2}(t)$ in terms of the solutions of the fourth Painlevé equation (18). A summary of the steps followed so far is presented in the diagram of Fig. 1, where the time-dependent Hamiltonians $\hat{H}_{1,2}(t)$ are discussed in Sec. 4.

As pointed out in the previous section, the shape-invariant condition (3) implies that $\hat{A}(t)$ and $\hat{A}^{\dagger}(t)$ are the ladder operators for the nonstationary eigenfunctions $\phi_{n}^{(1)}(x,t)$ of $\hat{I}_{1}(t)$. The latter indeed allows determining the spectral information (2), for $j=1$. To this end, we first determine the zero-mode eigenfunction, which is an element in the kernel of the annihilation operator $\mathcal{K}_{A}\equiv Ker(\hat{A}(t))=\{\phi^{(1)}\}$, with $\hat{A}(t)\phi^{(1)}=0$. However, in our case, the annihilation operator under consideration is a differential third-order one, and thus $\mathcal{K}_{A}$ is composed of three linearly independent zero-mode solutions, $\mathcal{K}_{A}=\{\phi^{(1)}_{0;1},\phi^{(1)}_{0;2},\phi^{(1)}_{0;3}\}$. Nevertheless, we must verify whether the elements in $\mathcal{K}_{A}$ fulfill the finite-norm condition. With the zero-modes already identified, the remaining eigenfunctions are computed from the iterated action of the creation operator $\hat{A}^{\dagger}(t)$ on the zero-mode eigenfunctions, and the respective eigenvalues increase by $2\lambda$ at each iteration.

For convenience, in this section we consider the case for which $\hat{M}(t)^{\dagger}$ factorizes as the product of two first-order operators, that is, a reducible case. Let us consider the factorization

$$\hat{M}^{\dagger}(t)\equiv\hat{M}_{1}^{\dagger}(t)\hat{M}_{2}^{\dagger}(t)\,,\quad\hat{M}(t)=\hat{M}_{2}(t)\hat{M}_{1}(t)\,,$$

(19)

where $\hat{M}_{1,2}(t)$ are first-order operators constructed in analogy to (9) as

$$\hat{M}_{1}^{\dagger}(t):=\left(\sigma\frac{\partial}{\partial x}-i\frac{\dot{\sigma}}{2}x+W_{1}(z)\right)\,,\quad\hat{M}_{2}^{\dagger}(t):=\left(\sigma\frac{\partial}{\partial x}-i\frac{\dot{\sigma}}{2}x+W_{2}(z)\right)\,.$$

(20)

The straightforward calculations show that the real-valued functions $W_{1}(z)$ and $W_{2}(z)$ are given by

$$W_{1}=-G+\left(\frac{G_{z}-\sqrt{-d}}{2G}\right)\,,\quad W_{2}=-G-\left(\frac{G_{z}-\sqrt{-d}}{2G}\right)\,.$$

(21)

From the latter result, it is clear that the factorization of $\hat{M}(t)$ requires $d<0$. Recall that $\hat{M}^{\dagger}(t)$ intertwines the quantum invariant $\hat{I}_{1}(t)$ with $\hat{I}_{2}(t)$. Thus, to inspect the respective intertwining relationships fulfilled by $\hat{M}_{1,2}(t)$ we introduce a new auxiliary quantum invariant

$$\hat{\mathfrak{I}}\tilde{\phi}_{n}(x,t)=\tilde{\Lambda}_{n}\tilde{\phi}_{n}(x,t)$$

(22)

with $\tilde{\mathcal{H}}=Span\{\tilde{\phi}_{n}\}_{n=0}^{\infty}$ the respective vector space composed with the finite-norm solutions. The spectral information of $\hat{\mathfrak{I}}(t)$ is not relevant, for it just serves as an aid to solve the eigenvalue problem associated with $\hat{I}_{1}(t)$. For this reason, the respective Hamiltonian associated with $\hat{\mathfrak{I}}(t)$ is not considered throughout the rest of the text.

The new auxiliary invariant satisfies the intertwining relationships

$$\hat{I}_{1}(t)\hat{M}_{1}^{\dagger}(t)=\hat{M}_{1}^{\dagger}(t)\hat{\mathfrak{I}}(t)\,,\quad\hat{\mathfrak{I}}(t)\hat{M}_{2}^{\dagger}(t)=\hat{M}_{2}^{\dagger}(t)\hat{I}_{2}(t)\,.$$

(23)

Given that both operators $\hat{M}_{1,2}(t)$ are of first-order, the relationships (23) are then equivalent to

		$\displaystyle\hat{I}_{1}(t)=\hat{M}_{1}^{\dagger}(t)\hat{M}_{1}(t)+\epsilon_{1}\,,\quad$		$\displaystyle\hat{\mathfrak{I}}(t)=\hat{M}_{1}(t)\hat{M}_{1}^{\dagger}(t)+\epsilon_{1}\,,$		(24)
		$\displaystyle\hat{\mathfrak{I}}(t)=\hat{M}^{\dagger}_{2}(t)\hat{M}_{2}(t)+\epsilon_{2}\,,$		$\displaystyle\hat{I}_{2}(t)=\hat{M}_{2}(t)\hat{M}_{2}^{\dagger}(t)+\epsilon_{2}\,,$

where the substitution of (20)-(21) into (24) leads to

$$\epsilon_{1}=\gamma-\sqrt{-d}\,,\quad\epsilon_{2}=\gamma+\sqrt{-d}\,.$$

(25)

From (7) and (23), we can see that the first-order operators define mappings among the vector spaces $\mathcal{H}_{1}(t)$, $\mathcal{H}_{2}(t)$ and $\tilde{\mathcal{H}}(t)$ in the following form:

		$\displaystyle\hat{M}_{2}^{\dagger}(t):\mathcal{H}_{2}(t)\rightarrow\tilde{\mathcal{H}}(t)\,,\quad$		$\displaystyle\hat{M}_{2}(t):\tilde{\mathcal{H}}(t)\rightarrow\mathcal{H}_{2}(t)\,,$		(26)
		$\displaystyle\hat{M}_{1}^{\dagger}(t):\tilde{\mathcal{H}}(t)\rightarrow\mathcal{H}_{1}(t)\,,\quad$		$\displaystyle M_{1}(t):\mathcal{H}_{1}(t)\rightarrow\tilde{\mathcal{H}}(t)\,,$
		$\displaystyle\hat{Q}^{\dagger}(t):\mathcal{H}_{2}(t)\rightarrow\mathcal{H}_{1}(t)\,,\quad$		$\displaystyle Q(t):\mathcal{H}_{1}(t)\rightarrow\mathcal{H}_{2}(t)\,.$

From the latter mappings, we construct the elements of $\mathcal{K}_{A}$ (zero-modes), that is, the eigenfunctions annihilated by $\hat{A}(t)$. We thus have

$$\hat{A}(t)\phi_{0;k}^{(1)}=\hat{M}^{\dagger}(t)\hat{Q}(t)\phi_{0;k}^{(1)}=\hat{M}_{1}^{\dagger}\hat{M}_{2}^{\dagger}(t)\hat{Q}(t)\phi_{0;k}^{(1)}=0\,,\quad k=1,2,3,$$

(27)

where it is worth discussing three different cases.

$\bullet$ $\hat{Q}(t)\phi_{0;1}^{(1)}=0$. Here, a first solution is determined, up to a normalization constant, by solving a trivial first-order differential equation.

$\bullet$ $\hat{M}_{2}^{\dagger}(t)\hat{Q}(t)\phi_{0;2}^{(1)}=0$ with $\hat{Q}(t)\phi_{0;2}^{(1)}\neq 0$. From the mappings defined by the intertwining relationship (7), it is clear that $\hat{Q}(t)\phi_{0;2}^{(1)}=\phi_{0}^{(2)}\in\mathcal{H}_{2}(t)$, with the latter being annihilated by $\hat{M}_{2}^{\dagger}(t)$. Thus, in order to determine the zero-mode $\phi_{0;2}^{(1)}$, we should solve the first-order differential equation $\hat{M}_{2}^{\dagger}(t)\phi_{0}^{(2)}=0$, and from it we determine $\phi_{0;2}^{(1)}$, together with the respective eigenvalue, after mapping $\phi_{0}^{(2)}$ through $\hat{Q}^{\dagger}(t)$ (see Fig. 2).

$\bullet$ $\hat{M}_{1}^{\dagger}(t)\hat{M}_{2}^{\dagger}(t)\hat{Q}(t)\phi_{0;3}^{(1)}=0$ with $\hat{M}_{2}^{\dagger}(t)\hat{Q}(t)\phi_{0;3}^{(1)}\neq 0$. In this case, the non-null element $\hat{M}_{2}^{\dagger}(t)\hat{Q}(t)\phi_{0;3}^{(1)}=\tilde{\phi}_{0}\in\tilde{\mathcal{H}}(t)$ is annihilated by $\hat{M}_{1}^{\dagger}(t)$. To extract the zero-mode $\phi_{0;3}^{(1)}$, we solve the first-order differential equation $\hat{M}_{1}^{\dagger}(t)\tilde{\phi}_{0}=0$. Then, we take $\tilde{\phi}_{0}$ to $\mathcal{H}_{1}(t)$ by consecutively performing the mappings $\hat{M}_{2}(t)$ and $\hat{Q}^{\dagger}(t)$ (see Fig. 2).

The complete procedure is summarized in the scheme depicted in Fig. 2. The straightforward calculations lead to

		$\displaystyle\phi_{0;1}^{(1)}(x,t):=\mathcal{N}_{0;1}^{(1)}\,\frac{e^{\frac{i}{4}\frac{\dot{\sigma}}{\sigma}x^{2}}}{\sqrt{\sigma}}e^{\int^{z}dz^{\prime}W(z^{\prime})}\,,$		(28)
		$\displaystyle\phi_{0;2}^{(1)}(x,t)=\mathcal{N}_{0;1}^{(1)}\left[W(z)-W_{2}(z)\right]\frac{e^{\frac{i}{4}\frac{\dot{\sigma}}{\sigma}x^{2}}}{\sqrt{\sigma}}e^{-\int^{z}dz^{\prime}W_{2}(z^{\prime})}\,,$
		$\displaystyle\phi_{0;3}^{(1)}(x,t)=\mathcal{N}_{0;3}^{(1)}\left[-2\sqrt{-d}+(W(z)-W_{2}(z))(W_{1}(z)+W_{2}(z))\right]\frac{e^{\frac{i}{4}\frac{\dot{\sigma}}{\sigma}x^{2}}}{\sqrt{\sigma}}e^{-\int^{z}dz^{\prime}W_{1}(z^{\prime})}\,,$

with the respective eigenvalues $\Lambda_{0;1}^{(1)}=0$, $\Lambda_{0;2}^{(1)}=\epsilon_{2}+2\lambda=\gamma+\sqrt{-d}+2\lambda$ and $\Lambda_{0;3}^{(1)}=\epsilon_{1}+2\lambda=\gamma-\sqrt{-d}+2\lambda$. The terms $\mathcal{N}_{0;j}^{(1)}$ stand for the normalization factors that might depend on time. From (28), the rest of the eigenfunctions are determined from the action $\left[\hat{A}^{\dagger}(t)\right]^{n}$, for $n=0,1,\cdots$, on each element $\phi_{0;j}^{(1)}$. By doing so, we generate at most three sequences of eigenfunctions, where the eigenvalues $\Lambda_{0;j}^{(1)}$, for $j=1,2,3$, increase by $2n\lambda$.

For the conventional stationary oscillator, it is well-known that the creation operator does not lead to finite-norm eigenfunctions. On the other hand, the one-step SUSY partner Hamiltonians admit a creation operator for which a finite-norm eigenfunction is achived. In the context of SUSYQM, such an eigenfunction is the so-called missing state [16]. Thus, it is natural to look for the solutions that are annihilated by the creation operator. In the case under consideration, we have constructed $\hat{A}^{\dagger}(t)$ as a third-order differential operator, which admits three linearly independent eigenfunctions, and at least one finite-norm solution is possible. The existence of the latter implies a truncation of the sequences generated from the zero-modes $\phi_{0;j}^{(1)}$. We thus define $\mathcal{K}_{A^{\dagger}}:=Ker(\hat{A}^{\dagger}(t))=\{\Phi_{0;1}^{(1)},\Phi_{0;2}^{(1)},\Phi_{0;3}^{(1)}\}$ as the set containing the finite-norm eigenfunctions of $\hat{A}^{\dagger}(t)$. If the set $\mathcal{K}_{A^{\dagger}}$ is empty, three infinite sequences are generated (see Sec. 6.2). In turn, if $\mathcal{K}_{A^{\dagger}}$ contains one single element, we generate at most two infinite sequences, together with one finite-dimensional sequence (see Sec. 6.3), which in particular could be a singlet (see Sec. 6.1).

Following the same steps as in (28), it is straightforward to show that the eigenfunctions of $A^{\dagger}$ are

		$\displaystyle\Phi_{0;1}^{(1)}(x,t)=\overline{\mathcal{N}}_{0;1}^{(1)}\,\frac{e^{\frac{i}{4}\frac{\dot{\sigma}}{\sigma}x^{2}}}{\sqrt{\sigma}}e^{\int^{z}dz^{\prime}W_{1}(z^{\prime})}\,,$		(29)
		$\displaystyle\Phi_{0;2}^{(1)}(x,t):=\overline{\mathcal{N}}_{0;2}^{(1)}\left[W_{1}(z)+W_{2}(z)\right]\frac{e^{\frac{i}{4}\frac{\dot{\sigma}}{\sigma}x^{2}}}{\sqrt{\sigma}}e^{\int^{z}dz^{\prime}W_{2}(z^{\prime})}\,$
		$\displaystyle\Phi_{0;3}^{(1)}(x,t):=\overline{\mathcal{N}}_{0;3}^{(1)}\left[\epsilon_{2}+2\lambda+(W_{1}(z)+W_{2}(z))(W_{2}(z)-W(z))\right]\frac{e^{\frac{i}{4}\frac{\dot{\sigma}}{\sigma}x^{2}}}{\sqrt{\sigma}}e^{-\int^{z}dz^{\prime}W(z^{\prime})}\,,$

where the respective eigenvalues are given by $\overline{\Lambda}_{0,1}^{(1)}=\epsilon_{1}=\gamma-\sqrt{-d}$, $\overline{\Lambda}_{0,2}^{(1)}=\epsilon_{2}=\gamma+\sqrt{-d}$ and $\overline{\Lambda}_{0,3}^{(1)}=-2\lambda$.

In general, we can not say which solutions in (28) and (29) fulfill the finite-norm condition, since it depends on the specific solutions of the fourth Painlevé equation. However, we may get more insight by considering the possible behavior of the asymptotics. To this end, let us suppose that the real-valued functions $W_{1}(z)$, $W_{2}(z)$ and $W(z)$ are smooth, and such that they converge to a finite value for $|z|\rightarrow\infty$. Then, finite-norm solutions are achieved depending on the convergence of the exponential functions in (28)-(29). For instance, if $e^{\int^{z}dz^{\prime}W(z^{\prime})}\rightarrow 0$ for $|z|\rightarrow\infty$, then $\phi_{0;1}^{(1)}(x,t)$ becomes a finite-norm solution, whereas $\Phi_{0;3}^{(1)}(x,t)$ do not. The same analysis can be extended to the rest of solutions in (28)-(29). In such a case, we may conclude that at most three out of the six solutions have a finite-norm. This indeed corresponds to the solutions discussed in Sec. 6.1 and Sec. 6.2.

So far, we have determined the families of exactly solvable quantum invariant $\hat{I}_{1}(t)$, related to the fourth Painleé transcendents, which fulfill a third-order SUSYQM shape-invariant condition. Nevertheless, the respective Hamiltonian of the system $\hat{H}_{1}(t)$ has not been identified yet. The latter is required to properly define the Schrödinger equation that characterize the quantum system under consideration. Such a task have been addressed in previous works using the factorization method for time-dependent Hamiltonians [56, 57] by imposing the appropriate ansatz. In this work, we consider an alternative approach based on the transitionless tracking algorithm [58, 62]. In this form, additional information is obtained about the classes of time-dependent Hamiltonians that can be constructed in term of the nonstationary eigenfunctions.

In App. A we have discussed the nonstationary eigenvalue equation of the quantum invariant associated with the parametric oscillator. Remarkably, the results from Lews-Riesenfeld [52] hold for any quantum invariant²²2Although, an orthogonal set of eigenfunctions can not be taken for granted for any general quantum invariant.. Therefore, for the quantum invariants $\hat{I}_{1,2}(t)$ constructed in Sec. 2 we can determine the respective time-dependent Hamiltonians $\hat{H}_{1,2}(t)$ such that the Schrödinger equation, in coordinate-free representation,

$$i\partial_{t}|\psi_{n}^{(j)}(t)\rangle=\hat{H}_{j}(t)|\psi_{n}^{(j)}(t)\rangle\,,\quad|\psi_{n}^{(j)}(t)\rangle=e^{i\theta_{n}^{(j)}(t)}|\phi_{n}^{(j)}(t)\rangle\,,\quad j=1,2,$$

(30)

is fulfilled, where the wavefunctions and nonstationary eigenfunctions are recovered from the coordinate-representation $\psi_{n}^{(j)}(x,t)=\langle x|\psi_{n}^{(j)}(t)\rangle$ and $\phi_{n}^{(j)}(x,t)=\langle x|\phi_{n}^{(j)}(t)\rangle$, respectively. Notice that the wavefunctions and eigenfunctions differ by just a time-dependent complex-phase, which is computed from (30) through the expectation value

$$\frac{d}{dt}\theta_{n}^{(j)}(t)=\langle\phi_{n}^{(j)}(t)|\left[i\partial_{t}-\hat{H}_{j}(t)\right]|\phi_{n}^{(j)}(t)\rangle\,,\quad j=1,2.$$

(31)

provided that the Hamiltonian is already known. However, in our case, both the Hamiltonian and the complex-phase are unknown, and a workaround should be implemented. To this end, we consider the time-evolution operator $\hat{U}_{j}(t;t_{0})$, that is, an operator that maps a solution defined at a time $t_{0}$ into one defined at a time $t$, $|\psi_{n}^{(j)}(t)\rangle=\hat{U}_{j}(t;t_{0})|\psi_{n}^{(j)}(t_{0})\rangle$. Given that both the Hamiltonians $\hat{H}_{j}(t)$ and the quantum invariants $\hat{I}_{j}(t)$ are self-adjoint, it follows that the time-evolution operator is unitary and it takes the diagonal form

$$\hat{U}_{j}(t;t_{0}):=\sum_{n=0}^{\infty}|\psi_{n}^{(j)}(t)\rangle\langle\psi_{n}^{(j)}(t_{0})|=\sum_{n=0}^{\infty}e^{i\left[\theta_{n}^{(j)}(t)-\theta_{n}^{(j)}(t_{0})\right]}|\phi_{n}^{(j)}(t)\rangle\langle\phi_{n}^{(j)}(t_{0})|\,,$$

(32)

where it is worth to recall that $\langle\phi_{n}^{(j)}(t^{\prime})|\phi_{m}^{(j)}(t)\rangle\neq\delta_{n,m}$, for $t^{\prime}\neq t$. Therefore, the set $\{|\psi_{n}^{(j)}(t)\rangle\}_{n=0}^{\infty}$ inherits the orthogonality from the set $\{|\phi_{n}^{(j)}(t)\rangle\}_{n=0}^{\infty}$. In this form, we can be built-up the vector spaces $\overline{\mathcal{H}}_{j}(t)=Span\{|\phi_{n}^{(j)}(t)\rangle\}_{n=0}^{\infty}$, which under the definition of inner-product are equivalent to the vector spaces $\mathcal{H}_{j}(t)$ introduced in Sec. 2.