Manifestation of instability in the quasiclassical limit of the spectrum of the quartic double well

TL;DR

This paper investigates how quantum spectra of a 1D quartic double well reflect classical instability as Planck's constant approaches zero, revealing signatures of the continuum limit and the persistence of tunneling near critical energies.

Contribution

It identifies quantum spectral signatures of classical instability in the quartic double well as 5 approaches zero, linking quantum behavior to classical chaos.

Findings

Signatures of classical instability appear in the quantum spectrum near critical energy.

Tunneling persists only at energies close to the critical energy.

Quantum spectrum approaches the continuum limit as 5 0.

Abstract

Within the Bohr's correspondence principle, the quantum theory should reproduce the classical world when $\hbar\to 0$. In practice, the discrete energies come close to each other before passing to continuum, causing that highly excited states become inaccessible. However it is possible to identify signatures that the spectrum is approaching the continuum limit. In this work, for $\hbar \to 0$, we focus our attention on the spectrum of the 1D quartic double well and find signatures of the classical instability in the quantum spectrum near the critical energy, and explore how the tunneling phenomenon remains only at energies close to this critical energy.

Tables

Table 1. Table 1: Number of states below the critical energy E c subscript 𝐸 𝑐 E_{c} for different values of ℏ Planck-constant-over-2-pi \hbar .

$\hslash$	1		$\frac{1}{10}$		$\frac{1}{100}$		$\frac{1}{200}$		$\frac{1}{500}$		$\frac{1}{1000}$		$\frac{1}{2000}$
No. states	10		94		950		1898		4746		9490		18980

Equations

$$H=\frac{p^2}{2m}+a{x}^2+b{x}^4,$$

$$J=\oint pdx=2{\int }_{x_1}^{x_2}pdx=2{\int }_{x_1}^{x_2}\sqrt{2}\sqrt{E-V(x)}dx,$$

$$\nu =\frac{1}{T}=\frac{\partial E}{\partial J},$$

$$T=\frac{\partial J}{\partial E}={\int }_{x_1}^{x_2}\frac{\sqrt{2}}{\sqrt{E-V(x)}}dx.$$

$$E_{k+1}^{\pm }-E_k^{\pm }=\hslash \omega =\frac{2\pi \hslash }{T}.$$

$$T=2\pi \hslash {\rho }^{\pm }(\bar{E}).$$

$$2T=2\pi \hslash \rho (\bar{E}),\text{(for E<Ec=0).}$$

$$T=2\pi \hslash \rho (\bar{E}),\text{(for E>Ec=0),}$$

$$\Psi (x)=\sum _{k=1}^N{c}_k{\psi }_k(x),$$

$${\psi }_k(\Omega ,x)=\frac{\sin (\pi (x-k\Omega )/\Omega )}{\pi (x-k\Omega )/\Omega },$$

$$V_{k{k}^{\prime }}$$

$$T_{kk}$$

$$T_{k\ne k^{\prime }}$$

$${\psi }_k(x)=\frac{(-1{)}^{N-k}}{(2h_N{)}^{1/2}}\frac{H_N(x)}{(x-x_k)}{e}^{-x^2/2},$$

$$V_{k{k}^{\prime }}$$

$$T_{kk}$$

$$T_{k\ne k^{\prime }}$$

$${\hat{a}}_{\Omega }=\frac{1}{\sqrt{2\hslash }}\left(\frac{\hat{x}}{\Omega }+i\Omega \hat{p}\right)$$

$${\hat{a}}_{\Omega }^{†}=\frac{1}{\sqrt{2\hslash }}\left(\frac{\hat{x}}{\Omega }-i\Omega \hat{p}\right),$$

$$\hat{H}=-\frac{\hslash }{4{\Omega }^2}({\hat{a}}_{\Omega }^{†}-{\hat{a}}_{\Omega }{)}^2+b{\left(\frac{\hslash }{2{\Omega }^2}\right)}^2({\hat{a}}_{\Omega }^{†}+{\hat{a}}_{\Omega }{)}^4+\frac{a\hslash }{2{\Omega }^2}({\hat{a}}_{\Omega }^{†}+{\hat{a}}_{\Omega }{)}^2.$$

$$H_{n,n}$$

$$H_{n-2,n}$$

$$H_{n-4,n}$$

$$J=\oint pdx=2\pi \hslash \left(n+\frac{\mu }{4}+\frac{d}{2}\right),$$

$$\mathcal{T}=\exp \left(-2\sqrt{\frac{2}{\hslash }}{\int }_{x_1}^{x_2}\sqrt{|E-V(x^{\prime })|}d{x}^{\prime }\right)$$

$$\lambda =\sqrt{-2a}=\sqrt{20}\approx 4.47214.$$

$$2\pi \hslash \rho (E)$$

$$T=\sqrt{\frac{2}{b}}{\int }_{x_1}^{x_2}\frac{dx}{\sqrt{(x-x_1)(x_2-x)(x+x_1)(x+x_2)}}.$$

$$x_{2,1}=\frac{\sqrt{-a\pm a\sqrt{1+4bE/a^2}}}{\sqrt{2b}}.$$

$$I(x,x_1,x_2)$$

$$\underset{x\to x_1}{\lim }\arcsin \left(\sqrt{\frac{(x_1+x_2)(x_2-x)}{2x_2(x_1-x)}}\right)=\underset{x\to x_1}{\lim }\arcsin \left(i\sqrt{\frac{x_2^2-x_1^2}{2x_2}}\frac{1}{\sqrt{x-x_1}}\right)$$

$$\approx \frac{i}{2}\left(\log \frac{2(x_2^2-x_1^2)}{x_2}-\log (x-x_1)\right),$$

$$\underset{x\to x_1}{\lim }I(x,x_1,x_2)=\frac{-2}{x_1+x_2}K\left(1-\frac{4x_1{x}_2}{(x_1+x_2{)}^2}\right),$$

Full text

classical correspondence, anharmonic oscillator, Lyapunov exponents

I Introduction

Since the beginning of quantum mechanics, potentials which admit analytical solution have played a key role to provide insights of the phenomena at the quantum scale. Right after them, we have non-analytically solvable problems which can be only accessed via approximation methods. One of the simplest non-analytically solvable problems is the quartic oscillator including both quadratic and quartic terms Lay (1997); Reid (1970); Bay and Lay (1997); Chen et al. (2003); Turbiner and del Valle (2021a); Turbiner (2010). In particular, for negative quadratic term added to a positive quartic term, we have the double well with a Mexican hat symmetric form which provides one of the simplest examples of symmetry breaking Joger et al. (2014); Delabaere and Pham (1997); Meier et al. (2018). The double well problem offers a good model to test theories of tunneling and to explore asymptotic approximations Fujimura and Sakai (2011); Joger et al. (2014); Meier et al. (2018). The ground and lowest states are well established and can be reached with a desired accuracy Fernández and Garcia (2017) by several approximation methods even though some of them show faster convergence than others (see discussion Okun and Burke (2021); Turbiner and del Valle (2021b) for instance). Less studied is the regime of highly excited states in the spectrum.

By its simplicity, the quartic potential also offers a good chance to explore quantum-classical correspondence Plata and Llorente (1992). For example, in Novaes (2003) a correspondence between the Husimi function and the classical trajectory for high quantum numbers was found and it was concluded that the zero point energy is relevant when dealing with low energy states. Similar correspondences have been observed for many-body models between the trajectory of the semiclassical approximation and the Husimi functions de Aguiar et al. (1991); Furuya et al. (1992); Nader et al. (2021); Arranz et al. (1996); Bakemeier et al. (2013). Apart from the Husimi and trajectories, classical correspondence can also be observed in the density of states Bastarrachea-Magnani et al. (2014); Stránský et al. (2014); Luna-Acosta et al. (2000); Puebla et al. (2016). However this can be achieved only if highly excited states are estimated to a certain accuracy.

Comparing quantum and classical analysis is specially valuable if some signatures of chaotic behavior can persist in the quantum analogous. The out-of-time ordered correlators (OTOC), which quantify the degree of non commutativity in time between two operators, have been widely used as a common tool to investigate chaos in quantum systems Maldacena and Stanford (2016); Wang et al. (2021); Roberts and Stanford (2015); Fan et al. (2017); Luitz and Bar Lev (2017); Torres-Herrera et al. (2018); Borgonovi et al. (2019); Yan et al. (2019). In general the exponential growth rate of the OTOC is related to positive Lyapunov exponents Chávez-Carlos et al. (2019), which is interpreted as high sensitivity to initial conditions i.e. connected to the notion of chaos. However, the OTOC can also grow exponentially fast in one-degree-of-freedom quantum systems that are not globally chaotic, but are critical Hummel et al. (2019a); Chávez-Carlos et al. (2019). This is the case of the quartic double well potential which exhibits periodical trajectories but it contains a critical energy at the local maximum.

In this manuscript we test different approximation methods in quantum mechanics in order to chose the most suitable option to access to highly excited states and compare them to a semiclassical method (EBK). Then we use the complete spectrum below the critical energy to observe classical correspondence in the density of states when $\hbar\to 0$. Finally we identify signatures of the positive Lyapunov exponent on the spectrum.

II Generalities

The Hamiltonian we are considering is the following

[TABLE]

where $p$ represents the momentum, $x$ the coordinate and $m$ the mass which, for simplicity, will be set from now on to the unity, $m=1$. In the quantum case, the variable $x$ and $p$ are promoted to operators $x\to\hat{x}$ and $p\to\hat{p}=-i\hbar\frac{d}{dx}\,$. The potential is plotted in Fig. 1 for particular constants $a=-10$ and $b=1$. It shows degenerated minima at $x=\pm\sqrt{\frac{-a}{2b}}=\pm\sqrt{5}$ and a local maximum at the center of coordinates $x=0$ which represents an unstable point of equilibrium at a critical energy $E_{c}=0$.

From a classical point of view, we can obtain the period of the orbits by considering action-angle coordinates

[TABLE]

where $x_{i}$ ($i=1,2$) are the classical turning points. The period in which the coordinates undergo an oscillation is given by

[TABLE]

and therefore the period for a particular value of the energy $E$ can be obtained by the following integral

[TABLE]

In Fig. 2 the classical period $T$ is plotted as a function of the energy. The divergence corresponds to the critical energy $E_{c}=0$, associated to the unstable equilibrium point at $x=0$, which physically implies that the particles at the double well potential take infinite time to reach the unstable equilibrium point.

In quantum mechanics, since the potential grows to infinite for $x\rightarrow\pm\infty$, there are infinite bound states, however the number of states with energy below the critical $E_{c}=0$ is finite. The Hamiltonian commutes with the parity operator, $\hat{P}\psi(x)=\psi(-x)$, consequently the energy eigenfunctions are symmetric ($P=1$) or anti-symmetric ($P=-1$ ), $\hat{P}\psi_{k}^{\pm}(x)=\pm\psi_{k}^{\pm}(x)$. Below the critical energy, the eigenstates organize in pairs ($P=\pm$) of quasidegenerate levels. The probability of tunneling between the two wells diminish as a function of the energy gap between these quasi-degenerate pairs Bhattacharyya and Pathak (1996); Yuste and Sánchez (1993), and is it weaker for states in the bottom of the wells. The energy gap between parity partners increases as their energy approaches the top of the double well, and so does the probability of tunneling. The energy gap between the lowest levels has been studied to a certain depth Zhou et al. (2003); Park et al. (1998) for different separations between the two potential minima. This type of solution to the Schrödinger equation leads to the concept of instantons Shuryak (1988); Zinn‐Justin (1984); Vaĭnshteĭn et al. (1982).

Within the approach of quantum mechanics, the frequency and therefore the period of the quantum dynamics inside either well at a certain energy region is given by the energy difference of consecutive levels with the same parity as

[TABLE]

By defining the energy density of states of a given parity as $\rho^{\pm}(\bar{E})=1/(E_{k+1}^{\pm}-E_{k}^{\pm})$, where $\bar{E}=(E_{k+1}^{\pm}+E_{k}^{\pm})/2$, the period as a function of the energy is given by

[TABLE]

The total density of states, taking into account the quasi double degeneracy of energy-level pairs with different parity is $\rho(\bar{E})=2\rho^{\pm}(\bar{E})$, therefore

[TABLE]

For states with energy above the double well, $E>E_{c}=0$, the double quasi-degeneracy is broken and the period and density of states are related as follows

[TABLE]

where $\rho(\bar{E})$ is the total density of states given by $\rho(\bar{E})=1/(E_{k+1}-E_{k})$ with $E_{k}$ the energy of levels irrespective of their parity.

Following the Bohr’s correspondence principle, the results from quantum mechanics, as $\hbar\to 0$, should reproduce those from classical mechanics where discrete energy levels or tunneling do not exist. However as $\hbar$ decreases, the number of bound states below the critical energy increases as $1/\hbar$, even though the potential remains fixed. In this case, large matrices are needed in order to obtain the spectrum below the critical energy. Accelerating the convergence of the energy levels is a priority to access to highly excited states close to the critical energy. In the following section we briefly describe approximation methods we used to obtain the spectrum from the ground state up to the critical energy $E_{c}=0$.

III Methods

Within the framework of variational methods, the simplest choice for building the trial function is to consider a linear superposition in a certain complete and orthogonal basis

[TABLE]

where $N$ is the size of the basis set and $\psi_{i}$ the basis functions. In the following subsections we briefly describe different basis.

III.1 Sinc Method

Within the Sinc method, the wavefunction is a linear superposition of sinc functions

[TABLE]

where $\Omega$ is the spacing between maximas of neighbouring sinc functions and $k$ is an integer index which controls the location of such maximas and takes values from $-k_{max}$ to $k_{max}$.

The matrix elements of the Hamiltonian evaluated in the set of sinc functions are

[TABLE]

The spacing parameter $\Omega$ is treated as variational and its optimal value is obtained by solving $\frac{d}{d\Omega}{\rm Tr}[\hat{H}]=0$. The size of the matrix representations is $N\times N$ where $N=2k_{max}+1\,$ is the number of sinc functions. For a detailed description we refer to Amore (2006).

III.2 Lagrange Mesh Method

In the LMM the wave function is composed by a superposition of Lagrange functions which satisfy the Lagrange condition $\psi_{k}(x_{k^{\prime}})=\lambda_{k}^{-1/2}\delta_{kk^{\prime}}\,,$ thus facilitating the integration of the matrix elements using the Gauss quadrature. In particular for the Hermite mesh, the Lagrange functions are

[TABLE]

where $H_{N}$ are the Hermite polynomials, $x_{i}$ the zeroes of $H_{N}$, $N$ the size of the mesh and $h_{N}=2^{N}N!\sqrt{\pi}$. The matrix elements for the potential and the kinetic terms are

[TABLE]

For details we refer to Amore (2006). For numerical calculations one can use the Mathematica package del Valle (2022).

III.3 Hermite basis

In this method the wave function is considered as a linear superposition composed by eigenfunctions of the harmonic oscillator i.e. Hermite polynomials weighted by a Gaussian function Verguilla-Berdecia (1993). The simplest way to construct the matrix representation in this basis, is rewriting the Hamiltonian (1) in terms of the usual ladder operators $\hat{a}=\frac{1}{\sqrt{2\hbar}}(\hat{x}+i\hat{p})$ and $\hat{a}^{\dagger}=\frac{1}{\sqrt{2\hbar}}(\hat{x}-i\hat{p})$ acting on eigenstates of the harmonic oscillator $|n\rangle$. However, the convergence rate of the energy in the truncated basis can be accelerated by introducing a scaling $x^{\prime}\to\Omega x$ Amore et al. (2005), which in turns transform the annihilation and creation operators as

[TABLE]

and

[TABLE]

respectively. The scaling parameter is optimized by $\frac{d}{d\Omega}{\rm Tr}\hat{H}=0$. In terms of scaled operators, the Hamiltonian reads

[TABLE]

Thus the matrix representation of the Hamiltonian is tridiagonal with non-zero elements

[TABLE]

and their corresponding transposed elements. The basis is truncated at $n_{max}$ and the size of the matrix is $N\times N$ where $N=n_{max}+1$.

III.4 Einstein-Brillouin-Keller

This method, named after Einstein-Brillouin-Keller Keller (1958) (EBK), belongs to the semiclassical methods since it is an improvement to the Bohr-Sommerfeld quantization. The quantization rule is the following

[TABLE]

where $J$ is the angle-action variable, $n$ is a non-negative integer (the quantum number), $\mu$ and $d$ are the Maslov index: $\mu$ corresponds to the number of turning points in the classical trajectory ($\mu=2$ in this case) and $d$ corresponds to the number of reflections with a hard wall (absent in this case $d=0$). Since EBK is a semiclassical method, it is expected to correctly reproduce the spectrum for large quantum number $n$ Hruska et al. (1997) and therefore to be a good approximation to access to highly excited states when $\hbar\to 0$.

IV Results

We performed numerical diagonalization of the Hamiltonian (1) for different values of $\hbar$ in order to approach the classical limit and identify some signatures of classical instability in the spectrum of the quartic oscillator. Our results are organized in four subsections. For simplicity we focus our attention to a quartic potential with constant values $a=-10$ and $b=1$.

IV.1 Convergence

We tested the convergence of the energy for different levels and different values of $\hbar$ as a function of the basis size $N$ for the following basis: (i) Sinc, (ii) LMM basis, (iii) Hermite functions and (iv) the Hermite functions with scaled coordinate. We measure the rate of convergence by computing the energy difference $\Delta E_{n}{(N)}=|E_{n}(N)-E^{ref}_{n}|$ where $n$ is the quantum number, $N$ the size of the basis and $E^{ref}_{n}$ a reference energy. We use as reference the results obtained by the corresponding method with the basis size $N=2000$. Since the three methods involve numerical diagonalization, we use the routine diasym of LAPACK and perform diagonalization in fortran90. Obtaining the matrix representation does not imply a major difference in computing time for any method. As for the optimization of the scaling parameters $\frac{d}{d\Omega}Tr\hat{H}=0$ it is necessary to solve a third (Hermite) or sixth-degree (Sinc) equation in $\Omega$. This process does not represent a significant addition on computational time even for large basis size. However LMM requires to obtain $N$ roots of the Hermite polynomials to build the mesh thus increasing considerably the computing time for large values of $N$.

In Fig. 3 we plot $\Delta E_{n}(N)$ as a function of the basis size considering the four basis. One can appreciate that the fastest convergence is always obtained with the Sinc basis, while the lowest is for LMM. One can notice that from $\hbar=\frac{1}{10}$ the convergence rate of Lagrange basis is very slow. This is mainly because the mesh points must be redistributed by introducing a scaling parameter (see Baye (2015)) and treating it as variational parameter. In LMM there is no closed expression to estimate its optimal value, instead numerical minimization must be performed to obtain its optimal value. Since the optimal value of the scaling parameter highly depends on the basis size $N$ and the quantum number $n$, the analysis is impractical for highly excited states. On the other hand, the optimization of the scaling parameter $\Omega$, corresponding to the spacing between maximas within the Sinc and squeezing within the Hermite basis, accelerates the rate of the convergence by using the closed expression $\frac{d}{d\Omega}{\rm Tr}{\hat{H}}=0$. Finally we compared the accuracy of the Sinc method to the semiclassical method EBK. It can be seen from Fig. 4 that EBK provides correctly at least six decimal digits in the entire region below the critical energy $E_{c}=0$ for $\hbar=\frac{1}{2000}$.

IV.2 Quantum-classical correspondence

As $\hbar\to 0$ it can be expected that the energy difference between discrete states vanishes. Therefore more bound states are contained below the critical energy even though the potential remains fixed. The number of states which lie below the critical energy can be counted using the EBK results and they can also be corroborated with the Sinc method. In Table 1, the number of bound states with energy below the critical for different values of $\hbar$ is presented.

Whereas the EBK does not have a limitation for the number of bound states, the quantum methods are limited by the size of the matrix to perform diagonalization. The minimum value reached in this work is $\hbar=\frac{1}{2000}$ where 18980 bound states lie below the critical energy.

In Fig. 5 we plot the continuous classical period Eq.(2) along with the discrete quantum density Eq.(5), for different values of $\hbar$. We notice that the discrete points from the quantum mechanics are distributed along the classical continuous period and the emergence of the divergence corresponding to the critical energy becomes evident as $\hbar\to 0$. Before the critical energy $E<E_{c}=0$, according to Figure 4, the difference between the energies obtained by diagonalization in the Sinc basis and the EBK is of the order $10^{-6}$ energy units, because of that in Fig.5 the EBK values (green dots) lie behind the results coming from the Sinc method(purple dots) and only a few of them can be appreciated in the region $E<0$.

IV.3 Tunneling decay

One of the most intriguing properties of the quantum theory is the possible quantum tunneling in the region of the double well which is classically forbidden. The energy gap between quasi-degenerated pairs of different parity is a consequence of the tunneling phenomenon which appears in the energy region inside the double well. It is well known that the phenomenon of tunneling is less likely to be observed if the wells are too separated from each other Zhou et al. (2003). In this case we kept the potential fixed but the tunneling effect in quantum mechanics should vanish when approaching the classical limit as $\hbar\to 0$. Our high accurate results, provided by the Sinc method, give us access to the energy gap between two consecutive energy levels corresponding to different parities. On the other hand we can also use our high accurate energies to estimate the transmission coefficient by using the WKB approximation

[TABLE]

where $x_{1}$ and $x_{2}$ are the classical turning points.

In Figure 6 we plot the energy gap of parity-partners energy levels $\Delta E$ at panel (a), and the transmission coefficient $\mathcal{T}$ at panel (b) as a function of their energy average $\bar{E}$ for different values of $\hbar$. It can be noticed that the energy gap and the transmission coefficient decreases exponentially as the energy becomes more negative and the levels go deeper into the double well. This can be understood because the energy barrier is narrower close to the local maximum at $x=0$ corresponding to the critical energy $E=0$. As $\hbar$ decreases, the energy gap and transmission coefficient become detectable only for energies close to the critical energy $E_{c}=0$. Indicating, thus, that tunneling phenomenon can only be observed close to the critical energy when approaching the classical limit. Given the linear behavior of the transmission coefficient and the energy gap of parity partner in the log-log scale, at Figure 6 panel (c), we are able to establish that the tunneling increases as a function of the energy gap as ${\mathcal{T}}\propto(\Delta E)^{\alpha}$ and decay as ${\mathcal{T}}\propto\hbar^{-\alpha}$, where $\alpha=5/2$.

IV.4 Manifestation of the positive Lyapunov exponent

From a classical treatment of the Hamiltonian (1) it can be easily shown that the positive Lyapunov exponent which characterizes the critical behavior of the systems at the stationary point $(x_{s},p_{s})=(0,0)$ $(\dot{x}=0,\dot{p}=0)$ is (see Appendix B and Pilatowsky-Cameo et al. (2020))

[TABLE]

On the other hand, it can be shown that the density of states, close to the critical energy, behaves as

[TABLE]

where $\lambda$ is the Lyapunov exponent. In Appendix A we discuss the derivation of (IV.4) in detail. Thus, the relation between the density of states and the logarithm of the energy is linear and the Lyapunov exponent $\lambda$ governs the slope in agreement with the results from Ref.Hummel et al. (2019b).

In Fig. 7 the density of states, for energies lower but very close to $E_{c}=0$, is plotted against the logarithm of the energy for the smallest value of $\hbar=1/2000$. We can observe that the data fits to a linear function with negative slope $\beta=-0.4192$. The slope given by the fit is similar to that in Eq. (IV.4) $\beta=-\frac{2}{\lambda}=-0.447214$ in terms of the Lyapunov exponent. This is an indication that the positive Lyapunov exponent governs the behavior of the quantum density of states close to the critical energy. We think that the deviation of the slope would vanishes for smaller values of the Planck constant $\hbar$ than those reached in this work.

V Conclusions

The Sinc method accelerate the rate of convergence of the energy levels, allowing us to explore highly excited states close to the critical energy when approaching the classical limit. The optimization of spacing parameter $\Omega$ plays an essential role to increase the rate of convergence. The highly accurate results by the Sinc method allow to explore how the tunneling phenomenon decay when approaching the classical limit and how it can be detected only close to the critical energy. On the other hand, the semiclasical method EBK provides at least six correct decimal digits for the entire energy region of interest in this work, thus facilitating the analysis of the density of states.

We observed that the highly excited states converge to the critical energy forming a logarithmic divergence in the density of states. The quantum results go to the classical results as $\hbar\to 0$ and they successfully reproduce the divergence on the classical period corresponding to the stationary point for the minimum value considered in this work $\hbar=1/2000$. Furthermore, we found a manifestation of the positive Lyapunov exponent on the density of states of the spectrum of the quartic harmonic oscillator. Therefore we can conclude that the Lyapunov exponent governs the logarithm divergent behavior of the density of states for energies close to the critical $E_{c}=0$.

VI Acknowledgments

DJN thanks the kind hospitality of the Brown University Department of Chemistry. This project was partially supported by Fulbright COMEXUS Project NO. P000003405.

Appendix A The classical period on the vicinity of the critical energy

The classical period given by (2) can be conveniently written in terms of the turning points as

[TABLE]

where $x_{1},x_{2}$ are given by

[TABLE]

The indefinite integral has the following form

[TABLE]

where $F(\phi|u)$ is a incomplete elliptical integral of first kind. It is clear that (15) vanishes when evaluating in $x_{2}$. Before evaluating (15) in $x_{1}$, we should explore the behavior of the arcsin function in the limit $x\to x_{1}$

[TABLE]

[TABLE]

taking only the leading term of the series expansion. By substituting this expression in the indefinite integral (15) we obtain

[TABLE]

where $K(u)$ is the complete elliptic integral of first kind. Substituting the limit in (13) we obtain the classical period in terms of the turning points

[TABLE]

[TABLE]

So far, $x_{1}$ and $x_{2}$ have been treated as arbitrary only depending on the energy (14). However we are interested in the turning points (14) when the energy approaches the critical $E\to 0$. In this limit case, we can use the approximation

[TABLE]

to obtain the leading term of the turning points $x_{2}\rightarrow\sqrt{\frac{-a}{b}}$ and $\lim_{E\rightarrow 0^{-}}x_{1}=\sqrt{\frac{|E|}{-a}}.$ For this turning points close to the critical energy, the period is

[TABLE]

Finally, considering the asymptotic behavior of the complete elliptical integral

[TABLE]

we find that the classical period behaves as

[TABLE]

which in turns can be writen in terms of the Lyapunov exponent $\lambda=\sqrt{-2a}$. Finally, by recalling that for $E<0$ the total (irrespective of the parity) density of states is $2\pi\hbar\rho(E)=2T$, we obtain Eq.(IV.4) of the main text.

Appendix B The Lyapunov exponent

In this appendix we discuss in detail how to obtain the value of the Lyapunov exponent corresponding to the stationary point ($x=0$,$p=0$) of the phase-space corresponding to the the quartic oscillator (1). The canonical Hamilton equations read

[TABLE]

Therefore we can write

[TABLE]

with

[TABLE]

A stationary stationary point satisfies $F(x_{0},p_{0})=(0,0)$, which is equivalent to $x(t)=x(t_{0})=x_{0}$ and $p(t)=p(t_{0})=p_{0}$. In order to find the Lyapunov exponent associated to $x_{0}$, let us build the Jacobian matrix $A=D_{x}F(x,p)$ and evaluate it at the stationary point Pilatowsky-Cameo et al. (2020)

[TABLE]

The eigenvalues of (19) are $\lambda=\pm\sqrt{-2a}\,$. The Lyapunov exponent corresponds to the maximum real part of the eigenvalues. Since $a=-10$ in our case of interest, the Lyapunov exponent is the positive eigenvalue.