Alternative quantisation condition for wavepacket dynamics in a hyperbolic double well

We propose an analytical approach for computing the eigenspectrum and corresponding eigenstates of a hyperbolic double well potential of arbitrary height or width, which goes beyond the usual techniques applied to quasi-exactly solvable models. We map the time-independent Schr\"odinger equation onto the Heun confluent differential equation, which is solved by using an infinite power series. The coefficients of this series are polynomials in the quantisation parameter, whose roots correspond to the system's eigenenergies. This leads to a quantisation condition that allows us to determine a whole spectrum, instead of individual eigenenergies. This method is then employed to perform an in depth analysis of electronic wave-packet dynamics, with emphasis on intra-well tunneling and the interference-induced quantum bridges reported in a previous publication [H. Chomet et al, New J. Phys. 21, 123004 (2019)]. Considering initial wave packets of different widths and peak locations, we compute autocorrelation functions and Wigner quasiprobability distributions. Our results exhibit an excellent agreement with numerical computations, and allow us to disentangle the different eigenfrequencies that govern the phase-space dynamics.


Introduction
Analytical modeling is widely used in many areas of physics. It provides key insight and interpretational power, which may be unavailable in purely numerical approaches. Although numerical models are versatile and extremely useful for quantitative comparisons, the physics involved may be difficult to extract. For that reason, analytic models are employed to establish paradigms, or distill the essential features of a physical system. For instance, the harmonic oscillator is widely used in many areas of physics, such as quantum optics, solid state physics, or molecular physics to describe modes of the electromagnetic field, lattice or molecular vibrations (see for example [1]). Within strong-field laser-matter interaction, the Gordon-Volkov solution [2,3] has been widely used approximate the electron dynamics by field-dressed plane waves. This solution is exact and constitutes an important ingredient in constructing the strong-field approximation, which is a semi-analytic approach that can be linked to interfering electron orbits [4]. An orbit-based interpretation was vital to the description of strong-field phenomena as the laser-induced rescattering or recombination of an electron with its parent ion [5,6,7]. This interpretation led to the inception of eigenvalues for the potentials proposed earlier by [18,19,20].
The approach developed here is then applied to molecular tunneling. This is motivated by recent studies of strong-field enhanced ionisation in stretched molecules, in which momentum gates in phase space have been identified using Wigner quasiprobability distributions. These gates allow a direct intra-molecular population flow and were attributed to the system's non-adiabatic response to a a strong laser field [21]. Recently, however, we have shown that momentum gates also exist for static fields, or even in the field-free case [22]. The key physical mechanism facilitating such gates is quantum interference, which provides a bridge for the electronic wavepacket to reach the other centre and ultimately the continuum. These quantum bridges perform a clockwise rotation in phase space, whose frequency depends on the initial wavepacket and the internuclear separation. However, it is yet not understood what properties of the system determine these frequencies. The analytical model developed in this work is ideally placed for an in-depth study of how the initial electronic wavepacket influences this motion, and how it is related to the system's eigenfrequencies. Specifically, a hyperbolic double well potential has several desirable properties for the molecular toy model. First, the limit V (x) → 0 as x → ±∞ allows for the existence of continuum of states for positive electron energies. This is in contrast to the models of double-well potential by e.g. [23] or [24] for which V (x) → ∞ as x → ±∞. Second, it allows to faithfully model the binding potential in the region of interest (i.e. close to the central barrier [22]). Third, the location of the (symmetric) wells and peak value of V (x) may be independently tuned. Finally, although other hyperbolic double-well models such as those developed by [19] or [20] can reliably model the central potential barrier, the one we are using leads to an impenetrable barrier by classical means. This is important to rule out other populationtransfer mechanisms.
This article is organised as follows. In Sec. 2, the method developed by us is outlined, including how the Schrödinger equation can be reduced to Heun's equation (Sec. 2.1), the quantisation condition we propose (Sec. 2.2), determining the number of bound states (Sec. 2.3) and how to construct appropriate wave packets and ascertain their time evolution (Sec. 2.4). Subsequently, in Sec. 3, we apply the model to tunneling dynamics, analysing the time profiles of autocorrelation functions (Sec. 3.1) and Wigner quasiprobability distributions (Sec. 3.2). In particular, we determine the main frequencies with which the above-mentioned quantum bridges propagate and their dependence on the initial wave packet. Finally, in Sec. 4, we conclude the paper and discuss possible future directions.

Methods
We will consider the evolution of a time-dependent wave packet Ψ(x, t) using the basis of eigenstates ψ n (x) that solve the time-independent Schrödinger equation (TISE) with the Hamiltonian defined byĤ The binding potential is a member of the wider family of symmetric hyperbolic potentials of the form where V (m) 0 specifies the depth of the potential and d its peak location. For m = 1, 2 they produce a double-well (bistable) potential and for m = 0 they reduce to the (single-well) Pöschl-Teller potential, which is exactly solvable [25].
In the eigenbasis of the TISE, where Λ n = Ψ(x, 0)ψ * n (x)dx (6) are overlap integrals between the initial wavepacket Ψ(x, 0) and eigenfunctions ψ n (x) of the hyperbolic double well. The goal of the present investigation is to determine Λ n by analytical means. The 1D time-independent Schrödinger equation (TISE) for the potential given in Eq. (3) reads as with dimensionless parameters z = x/d, U 0 = 2mV 0 / 2 , = 2mE/ 2 . Note that the potential V (x) has even parity, which implies the existence of even/odd parity wavefunctions.

Reduction of Schrödinger's to Heun's equation
It was shown [13] that for even parity wavefunctions the above equation may be reduced to the Heun confluent differential equation by introducing the new variable ξ = 1/ cosh 2 z (with 0 < ξ ≤ 1 as −∞ < z < +∞), that is, where the even-parity solutions to TISE are of the following form: At this point we should make few technical comments. First, note that the x → ξ mapping is not injective and uniquely represents only the half of the −∞ < x < +∞ range. However, the other half of the range is just the symmetric copy of the first one -this mapping intrinsically constrains the wavefunctions to be even in x-variable space. Second, as Heun's confluent equation arises from Schrodinger's equation, some of the parameters out of α, β, γ, δ, µ, v are dependent on each other: in fact only α (property of the depth and width of the potential) and β (free parameter which determines the allowed energies) are independent.
Hence, we may write H(α, β, γ, δ, η, ξ) as the following infinite power series involving only two parameters 2 with the radius of convergence |ξ| < 1 given by Poincaré-Perron theorem [26]. The above power series is supplemented with the three-term recurrence relation with initial conditions v 0 = 1, v −1 = 0 and parameters The solution to Heun's differential equation is given in terms of the infinite power series with the coefficients determined by the above three-term recurrence relation. Unfortunately, it is not possible to find the explicit formula for v n solving the recurrence relation in Eq. (11). Instead the solution may be provided in terms of the holonomic sequence 3 . For a fixed α parameter, v n (α, β) are polynomials in β with their degree increasing with n.
Examples of these polynomials are displayed in Fig. 1 for the even wavefunctions. This figure shows a rather surprising fact that the different degree polynomials in β have their roots for almost exactly the same β values. This is crucial as it means that if β is chosen such that the corresponding N th order polynomial attains 0, all higher order polynomials n > N will be "very close" to 0 too, "very close" being further quantified as O(1/N ).
The infinite power series given in Eq. (10) may be terminated to a polynomial of degree N if both C N +2 = 0 and v N +1 = 0 conditions are simultaneously satisfied [13,26]. In such case the system becomes quasi-exactly-solvable and a subset of its eigenvalues may be found explicitly. However, imposing two equations on one free parameter β in the model implies that the other equation must constrain the value of α which represents well location and depth of the potential. This means that, for a fixed well location, terminating to a polynomial approach will be permissible only for selected values of its depth and viceversa. Unfortunately the above-mentioned constraints did not result in a choice of physically relevant parameters and only one eigenvalue per choice of α may be found [13]. This effectively precludes the calculation of all overlap integrals between the bound-states and the arbitrary initial wavepacket placed in the hyperbolic double-well.

Quantisation condition
However, it is evident that Schrödinger's equation should provide us with the solution for a full-range of depths and well locations of the potential. Instead of working in a quasi-exactlysolvable framework and terminating the Heun power series to a polynomial, we propose an alternative approach. We suggest that the entire eigenspectrum may be found by ensuring that the infinite power series converges to 0 sufficiently quickly such that the wavefunctions are still square integrable. This possibility stems from the asymptotic (discarding 1/n 2 terms) behaviour of the holonomic sequence v n : it may be found empirically that the values of its terms for large n significantly depend on the energy quantisation parameter β.
We propose that, for a given parameter α of the potential, the admissible values of an energy quantisation parameter β = β crit (with β > 0) correspond to the roots of the polynomial v n (β crit ) for large value of n (or more strictly as n → ∞). Thus, in practice, the problem of finding the allowed energies in the hyperbolic well problem boils down to finding the roots of a certain (usually) high degree polynomial in β. Based on the above claim we find the energy eigenvalues via numerical root-finding methods. This quantisation criterion forms a more efficient alternative to a typically used condition based on Wronskians [19,20,28] and is similar to the one achieved in [29] on different grounds. Furthermore, whenever α is chosen such that the infinite power series terminates, it is clear from Table 1 that the eigenvalues found using the above-mentioned claim are arbitrarily close to the ones given by the explicit analytical formula.
In other words, we propose that the sequence of the power series coefficients {v n (α, β crit )} in Eq. (10) decreases with n quickly enough for allowed values of β such that the wavefunc-  Table 1: Comparison of the eigenenergies obtained using the quantisation condition in section 2.2 and analytic formula obtained by [13]. The precision displayed corresponds to 20 significant figures. tion is square-integrable in the range ξ = 0 to ξ = 1. In a typical setup this is equivalent to demanding that ψ even (x → ±∞) → 0 which corresponds to ψ even (ξ → 0) → 0. Interestingly, in the present case, ψ even (ξ → 0) requirement is trivial. On the other hand the ξ → 1 limit lies just "on the edge" of the radius of convergence. The ξ → 1 (corresponding to x = 0) limit may be investigated using Abel's theorem [30]. The value of the power series as ξ → 1 should approach ∞ n v n (α, β) provided that ∞ n v n (α, β) converges, which requires choice of β (for a fixed α) such that v n (β) < 1/n for large n (by direct comparison test). However, it should be noted that convergence of ψ even (ξ) at ξ = 1 is by no means a sufficient criterion for wavefunction square-integrability. Unfortunately, as v n coefficients are not given by an explicit formula, it seems to be burdensome to evaluate the square-integrability constraint directly to find the admissible values of β -especially that the asymptotic (large n) behaviour of terms v n as functions of its parameters, to our best knowledge, is not well-understood [26]. Therefore, we rely on the indirect arguments presented below.

Argument based on the recurrence relation
To motivate the above-stated claim we propose the following argument based on the recurrence relation (11). Consider the particular term n = N of Eq. (11). By β crit n we will denote such value of β, for which the v n (α, β crit n ) = 0. Suppose that v N −1 (α, β crit N −1 ) = 0. Then, using Eq. (11), we obtain Evaluating the right-hand-side for large N (i.e when α/N α α + 2β crit Note that −1 < α/ N + β crit N −1 < 0 provided that α < 0 and β < |α| which is indeed fulfilled as for the normalisable solution to Schrodinger equation we require E > min x V (x). As it was stated in the previous section, the necessary criterion for square-integrability is that the sequence {v n } must converge to 0 for large n quicker than 1/n for large n. Therefore, we have that: with α/N 1 hence mimicking the termination of the power series to a polynomial through zeroing two subsequent-terms in the three-term recurrence relation [31]. Note that as in principle we can make N to be arbitrarily large, the error associated with the non-exact truncation of the power series can be made arbitrarily small. Furthermore, it may be readily seen that for β chosen such that the polynomial v N −1 (β) = 0, all polynomials for n ≥ N − 1 can be factorised into where q n (β) is another polynomial. Such factorisation property resembles the result from the theory of quasi-exactly-solvable models due to [32]. However, in contrast to [32], the zeroing of the "critical polynomial" v N −1 (β) does not imply that all subsequent terms will vanish but instead they will pick up a very small factor of α/N 1.

Argument based on smoothness of the wavefunction
A more strict argument is given by the smoothness of the wavefunction. Although the x → ξ mapping intrinsically constrains the wavefunctions to be even, it is not guaranteed that the wavefunction produced by joining of the two half-space wavefunctions will be 'smooth'. However, it is reasonable to demand that their derivatives should be continuous everywhere [33], i.e., that ψ (x = 0) = 0 (corresponding to ξ = 1). This can be expanded to produce The first term of the expression will always be 0 but the second term needs to be finite for the product to be 0. For large n we can omit the constants in the second term which becomes: ∞ n=0 v n (α, β)n. Now this expression will converge if v n goes to 0 quicker than 1/n 2 (by direct comparison test). In other words, the wavefunction ψ even (x(ξ)) will be acceptable only for such β and certain constant c that |v n | < c/n 2 for sufficiently 4 large n . This is a stronger condition for large n than previously stated.
Empirical results in Fig. 2 indeed confirm that the v n (β = β crit ) < 1/n 2 with β crit given by the quantisation condition from section 2.2. Furthermore, values of β slightly away from β crit do not fulfill this criterion which suggests that the v n < 1/n 2 bound is already tight.
Furthermore, the quantisation condition may be rewritten by using a well-established link between continued fractions and three-term recurrence relations. Interestingly, such formulation of the quantisation condition is closely related to a quantisation condition proposed by Manning (1935) ( [18], p. 137, Eq. (7)) for the potential bearing his name -also a member of a hyperbolic double-well family. The infinite continued-fraction formulation is very convenient for numerical implementations. For details see Appendix A.

Wavefunctions and accuracy of a non-exact series truncation
Here we illustrate that the error attained with the non-exact finite-order truncation of the infinite power series Eq. (10) is negligible. The first clue comes from the factorisation property discussed in section 2.2.1: all terms beyond N − 1 pick up an additional α/N 1 factor if N is sufficiently large -and hence should contribute very little to the shape of the Heun function.
On a more practical side, for the parameters chosen (V 0 = 74.785, d = 1 ⇒ α = −12.229), based on Fig. 3, it is clear that for terms n 12 for the even wavefunctions and n 28 for the odd wavefunctions the truncation error may be safely neglected -thus making the non-exact truncation to be extremely accurate and computationally feasible. The set of all bound-state eigenfunctions for parameters V 0 = 74.785, d = 1 are shown in Fig. 4. The total number of bound states is discussed in the following section.

Number of bound states
It is not straightforward to predict how many bound states we should expect as functions of V 0 and d of the potential, without explicitly invoking the proposed quantisation condition.
Here, we instead employ the theoretical lower and upper bounds on the number of bound states B. Ref. [34] has given the following upper-bound on the number of bound states in 1D, which was later used in context of hyperbolic-well potentials by [35] B . Furthermore, we recall the well-known theorem stating that the arbitrarily weak potential in one and two dimensions fulfilling V (x) ≤ 0 for all x and +∞ −∞ V (x)d n x < 0 (n = 1, 2) will have a bound state (see e.g. [34] p. 2) for reference). Therefore, we conclude that the number of bound states in a hyperbolic double-well potential is bounded by It may be clearly seen that for the hyperbolic-double-well potential we expect the number of bound states to be finite, with an upper bound growing like (V 0 d) 1/2 . For parameters of the potential V 0 = 74.785 and d = 1 this results in the number of bound states 1 ≤ B ≤ 10 which provides rather tight bounds on the actual number of bound states found from a quantisation condition: B = 6.

Initial wavepackets, overlap integrals and temporal evolution
Having found the eigenfunctions and eigenvalues of the TISE, we now formulate the initial wavepackets to be placed in a hyperbolic double-well. The purely even/odd (when mapped c=7, Ω=1/10 c=20, Ω=1/10  Figure 5: Graphs of the ψ DLE (x, W, τ ) (left) and ψ DLO (x, W, τ ) (right). Note that in both cases width and peak location of the wavepackets may be modified independently.
to the x-space) wavepackets may be divided in ξ-and ζ-spaces respectively. We would like to benefit from the relatively simple forms of the eigenfunctions in ξ-and ζ-spaces and to devise the simple purely even/odd initial wavepackets in ξ-and ζ-spaces. In this way the intricate 5 overlap integrals from x-spaces can be evaluated in ξ-and ζ-spaces by introducing the appropriate measures (weight functions) q(ξ) and Q(ζ) to the integrals. In ξ-space, noting that dξ/dz = −2 sinh (z)/ cosh 3 (z) we obtain whereas in ζ-space we have: dζ/dz = 1/ cosh(z) 2 = 1 − tanh(z) 2 = 1 − ζ 2 and therefore
• ψ DLE (x, c, Ω) contains two parameters which independently specify width and location of the peak. This is a result of a particularly simple relation between the location of a peak of a β-distribution as a function of parameters p and c: Therefore we set p = Ωc with Ω solely specifying the location of the peak and c solely specifying the wavepacket width.
The odd parity wavepacket is given by with W < 1 2τ 2 and 0 < τ < 1 and Γ(u) denoting a complete gamma function and 1F1 (a, b, c) regularised confluent hypergeometric function. Its choice has been motivated by the properties stated below: • The functional form of the odd parity wavepacket is proportional to ζ 1 − ζ 2 P e −W ξ 2 .
Such form may be motivated by noting that the odd wavepacket should have zeros at ζ = ±1 (corresponding to x = ±∞) and at ζ = 0 (corresponding to x = 0). Constraint W < 1 2τ 2 stems from demanding finite value of ψ DLO at ζ = ±1 (and hence P > 0) along with noting that P is a decreasing function of W for −1 < τ < 1 (see bullet point below).

Arbitrary wavepackets and a temporal evolution
Wavefunctions φ(ξ) from ξ-space and φ(ζ) from ζ-space are intrinsically mapped to even and odd x-space wavefunctions respectively. To produce arbitrary localised states we note that because functions in ξ-and ζ-spaces are orthogonal to each other when evaluated in xspace, we can just form the linear combinations of the initial wavepackets from both spaces to get neither purely even nor purely odd initial wavepacket in x-space. In such case, for the wavepacket of the form Ψ g (ξ, ζ, 0) = cos (∆)ψ DLE (ξ) + sin (∆)ψ DLO (ζ), the overlap integrals become , 0)ψ n (x)dx = Λ n cos(∆) for n = 0, 2, 4 Λ n sin(∆) for n = 1, 3, 5 where we have used the fact that ψ DLE (ξ(x)) and ψ DLO (ζ(x)) are respectively purely even and odd in an x-space. Furthermore, it should be noted, that as it is possible to calculate the overlap integrals for arbitrary values of (c, Ω) and (W, τ ) (subject only to constraints imposed in the previous sections), almost an arbitrary wavepacket in the x-space may be formed by making the linear combinations of the purely even (or odd) wavepackets with the fixed width c (or W ) and varying peak location Ω (or τ ).

Applications to tunneling dynamics
Next we will apply the analytical model developed here to the tunneling dynamics of an electronic wave packet, with focus on non-adiabatic temporal evolution. Our motivation is related to the presence of momentum gates, which have been first identified in the context of strong-field enhanced ionisation in position-momentum phase space using Wigner quasiprobability distributions [21]. Momentum gates are lines of approximately constant momentum through which there is a direct intra-molecular quasiprobability flow from one molecular centre to the other. In [21], they have been attributed to the non-adiabatic effect of a transient electron localisation at one of the wells due to the presence of a strong laser field. Such behaviour was further expounded by [22] who have shown that the timedependent field is not a necessary prerequisite for the momentum gates to occur and that the strong quasi-probability transfers may occur through "quantum bridges". These are highly non-classical, cyclic structures that form due to quantum interference. The aim of this section is to quantify this evolution for different initial wave packets, both in time and phase space.

Temporal evolution of the wavepacket
Now, we exploit the ability to calculate the overlap integrals in the hyperbolic-double well. This may be done in few steps: 1. We fix the parameters of the hyperbolic-double-well potential and find the eigenvalues E n (and hence eigenfrequencies ω n ) based on the quantisation condition provided in section 2.2.
2. We find the eigenfunctions given by Eqs. (9) and (12) for the values of β (allowed energies) found previously.
3. We devise the initial even/odd wavepackets according to Eqs. (18) and (19). 4. We calculate the overlap integrals in terms of incomplete gamma functions Γ(a, u) as detailed in section 2.4.2.
The results of applying this procedure are presented here for two different sets of parameters, each corresponding to a different limit behaviour. The time evolution of a wavepacket can be inferred using an autocorrelation function where Λ n are the overlap integrals defined in Eq. (6). Therefore Thus, any time dependence of |a(t)| 2 will stem from the differences in eigenenergies. Note that if |Λ n | = 0 and |Λ m | = 0 only for one pair of n and m with n = m then |a(t)| 2 will oscillate with a single frequency. Otherwise, the time evolution will be more involved.
First we devise an even-parity, delocalised wavepacket. The absolute value squared of the autocorrelation function |a(t)| 2 is displayed in Fig. 6 for wavepackets of different widths. They have been computed analytically using the method developed above, and numerically using the method in [22]. The agreement is excellent, with the analytical and numerical curves being practically indistinguishable and the temporal behaviours depending critically on the width.
In Fig. 6(a) this behaviour is quite intricate with two main frequencies: ω 20 = 4.73 a.u. and ω 40 = 7.46 a.u.. This is due to the coupling of n = 0 with n = 2, and n = 0 with n = 4 eigenstates (as Λ 2 and Λ 4 are small the n = 2 with n = 4 coupling may be safely neglected). In contrast, in Fig. 6(b), one can identify a single frequency for |a(t)| 2 , namely ω 20 = 4.73 a.u., which corresponds to only one pair of states with non-vanishing overlap integrals: Λ 0 and Λ 2 . Finally, the straight horizontal line in panel (b) corresponds to an initial wavepacket being very close to an eigenstate. As expected, this leads to a constant |a(t)| 2 within the precision used here. Minor discrepancies between the analytical and numerical results are related to the former not including overlaps with scattering states.
This critical behaviour is also observed for initially localised wavepackets such as those presented in Fig. 7(a). The corresponding values of |a(t)| 2 are displayed in Fig. 7(b) and are quite distinct. The blue curves in both panels illustrate the scenario in which only the overlap integrals Λ 0 and Λ 1 are non-vanishing. In contrast, the red curves show a slightly different wavepacket which gives rise to several different frequencies in the modulus squared of the autocorrelation function. One should note that for localised wave packets there are contributions from both even and odd eigenstates, which may lead to pulsated high-frequency oscillations enveloped by a slow oscillation.
We produce the localised wavepackets by making different linear combinations of initial wavepackets: ψ DLE (ξ) and ψ DLO (ζ). Various modes of behaviour are displayed in Fig. 7. Note that as the initial wavepackets (for the choices of parameters made) strongly overlap with the ground/first-excited states any short-scale oscillations in the |a(t)| 2 are enveloped with the long-scale oscillation of a period T = 2π/(E 1 − E 0 ) ≈ 520(a.u.).

Phase-space dynamics
Next we will investigate the wave packet's phase space evolution, with emphasis on the quantum bridges and their periodic motion. For that purpose, we will employ Wigner quasiprobability distributions. They are given by where the position and momentum coordinates are represented by x and p, respectively. Eq. space if integrated over the momentum or position coordinates, respectively. One should note, however, that Eq. (21) can be negative making it is a quasiprobability distribution. For more details on quantum systems in phase space see, e.g., [36]. In the analytical model W (x, p, t) may be calculated by numerical integration of Eq. (21), with the temporal evolution of the wavepacket Ψ(x, t) given by Eq. (5). The wavepacket, the eigenenergies and the eigenfunctions are calculated analytically as discussed in the previous sections. Throughout, we will focus on the scenario for which the quantum bridges are strong, namely initially delocalised wave packets and intermediate internuclear separations. The results comparing the present analytical model and the numerical results in [22] are displayed in Figs. 8 and 9. Fig. 8 corresponds to an initial wave packet leading to a single oscillation frequency in the autocorrelation function, while in Fig. 9 a more involved scenario with superimposed oscillations is explored. Overall, the agreement between the numerical and analytical results is excellent, which once more shows that the present model is reliable and, in contrast to the numerical approach in [22], can be used to determine the temporal evolution of the quantum bridges exactly.
In Fig. 8, we display the Wigner quasiprobability distribution computed using the initial wavepacket in Fig. 6(b). The figure shows a quasiprobability flow from one centre to the other, with a strong "quantum bridge" near p = 0. As the time flows, there is a motion of frequency ω 20 = 4.73 a.u., which corroborates the statement that only the overlap integral between the ground and second excited state is relevant to the problem at hand. The plot corresponds to almost a whole period of the autocorrelation function, and illustrate an oscillating behavior in the Wigner quasiprobability distribution. The bridges become slanted, change slope and then return to their original configuration at T ≈ 1.    6(a). The quasiprobability flow behaves in a much more convoluted way, with additional maxima near the quantum bridge and in both wells. For longer times, there will also be tails in the Wigner functions moving away from the potential wells, which indicate an overlap with a delocalised eigenstate, or in some cases ionisation. These tails are visible in the bottom panels of Fig. 9. For a detailed discussion of tails of Wigner functions in the context of strong-field ionisation see our previous publications [37,22]. From the autocorrelation function, we expect that the frequencies ω 20 and ω 40 will play a role. This convoluted behavior will be discussed in the rightmost column of Fig. 9, in which, instead of constructing Wigner quasiprobability distributions using the full analytical wavefunction, we consider only a coherent superposition between the ground and second excited state, where the overlap integrals Λ n (n = 0, 2) are given by Eq. (6). The partial Wigner quasiprobability flow mirrors the overall behavior reported in the central column of Fig. 9 except for the substructure and the tails. It determines the existence of the quantum bridges and their slopes, whose time evolution has the frequency ω 20 . This shows the dominance of this specific coupling and is expected, as tunneling should be dominated by the lower frequency. However, a modulation is introduced due to the non-vanishing overlap between the ground and the fourth excited state and its higher frequency ω 40 . Furthermore, the tails are absent in the partial results. This is due to the missing overlap integral with the fourth excited eigenstate, which is significantly broader (see Fig. 4).

Conclusions
In this work, we present an analytical method for solving Schrödinger's equation in a hyperbolic double well potential of any height and width. Our approach allows us to determine the entire eigenspectrum and corresponding eigenfunctions for the system up to an arbitrary precision, in contrast to finding only few individual eigenstates and eigenenergies in a related approach of quasi-exactly-solvable models [13,23,38]. By means of the exchange of variables x → ξ and x → ζ along with exploiting the parity of the potential, Schrödinger's equation for the hyperbolic-double well is reduced to Heun's equation [13] with the resulting wavefunction involving Heun's infinite power series. Instead of truncating this series to a polynomial [26,13] we avoid the problem of constraining the height/width of the potential by focusing on the series' convergence. This leads to a quantisation condition which reduces a problem of finding allowed energies to finding roots of a high-degree polynomial with coefficients generated from a three-term recurrence relation. The proposed quantisation condition displays some similarities with the results of the theory of quasi-exactly-solvable models, in particular sharing the same polynomial factorisation property [32]. However, it is more general as it gives a whole spectrum instead of a small subset of eigenvalues.
Using the initial wavepackets with independently tunable width and peak location, we calculate the overlap integrals with the system's eigenstates in terms of incomplete gamma functions. This allows us to analytically evaluate temporal evolution of the wavepacket as a function of its initial parameters. This method is then employed to study tunelling through a central barrier for different initial wave packets for different coherent superpositions involving two or more eigenstates. Apart from an excellent agreement with the numerical model in [22], which was used as a benchmark, this analytical model provides far more insight about the system's dynamics. Specifically, the autocorrelation functions and Wigner quasiprobability distributions exhibit a periodic motion that can be precisely determined using the system's eigenfrequencies. These dynamics are strongly dependent on the width of the initial wavepacket, and the time-independent overlap integrals obtained for an eigenfunction basis has greatly facilitated our studies. In addition to that, the present phase-space studies support the conclusions in [22] that the intra-center quasiprobability flows caused by quantum interference, dubbed 'quantum bridges' in our previous publication, have their time evolution determined by frequencies intrinsic to the system, instead of a non-adiabatic response to an external driving field as proposed in [21]. Moreover, for the specific, field-free case studied in this article, we have determined such frequencies exactly for a hyperbolic double-well, thus going beyond the rough estimates in [22].
The analytical model developed here may form a basis for investigating a wide-range of static or time-dependent perturbative effects and be helpful in testing predictions of more realistic but non-analytically-solvable models of a double-well. In particular, the model could analytically address the issue of finding the optimal parameters for enhanced ionisation [22] in a time-dependent field. For that purpose, it will be necessary to overcome a series of obstacles. First, the model developed in this article is strongly reliant on parity and inversion symmetry. Adding even a static field would break this symmetry and require changes in the way the eigenstates are calculated. Second, ionisation would require the computation of continuum states, which are not yet available in the present model. Third, a time dependent field would imply that the time-dependent Schrödinger equation may no longer be reduced to an eingenvalue equation. Hopefully, a low enough driving-field frequency may allow for a quasi-static picture with an effective potential and time-dependent dressed states.
The quantisation condition proposed here may be successfully applied to a wider class of potentials than the one given by Eq. (3). Although the arguments developed in Section 2.2.1 exploit the specific form of the recurrence relation, it should be possible to extend it to the other potentials for which the Schrödinger equation may be solved in terms of Heun's infinite power series generated from a three-term recurrence relation. In particular we verified that for the distinct symmetric hyperbolic potential proposed by [19] the quantisation condition predicts the eigenvalues of E = −1, −0.19113 and E = −1.0048, −0.25 for parameters {V 1 = 1, V 2 = −6, V 3 = 6} and {V 1 = 1, V 2 = −7, V 3 = 27/4} respectively, in consonance with the ones earlier reported ( [19], p. [4][5]. Furthermore, it appears that the proposed condition may be also applied to the asymmetric hyperbolic double-well potential [20], predicting the eigenvalues of E = 0.311, 2.434, 3.875 for parameters {w 1 = 15, w 2 = 12, w 3 = 1} in agreement with the Wronskian's method used by [20]. However, for asymmetric wells we do not expect the wavefunctions to have even/odd parities, Hence, a modified procedure would have to be applied to find the suitable initial wavepackets used for temporal evolution. Such asymmetric double-well potential could, for example, model the dynamical behaviour of the wavepacket in heteronuclear molecule setups.
for n ≥ 1 and b 0 = 0. Given value of α we can numerically search for such β which fulfills the above condition by a process of successive approximations. The proof of the above statement is presented below.
The value of the infinite continued fraction of the form (following notation used by ( [39],p. 28)) may be written as In such case, M n and L n fulfill the following three-term recurrence relations ( [39], p. 28): M n = b n M n−1 + a n M n−2 L n = b n L n−1 + a n L n−2 differing only by initial conditions: M −1 = 1, M 0 = 0, L −1 = 0, L 0 = 1. At this point we recognise that recurrence relation for L n is equivalent to Heun's recurrence relation for v n (Eq. 2.6) if we choose a n (α, β) = C n (α, β) A n (α, β) and b n (α, β) = B n (α, β) A n (α, β) for n ≥ 1. Hence, we conclude that searching for such β that v n (α, β) = 0 for large n based on (Eq. 2.6) is equivalent to searching for a root of the 1/x infinite continued fraction (corresponding to a n → ∞ limit). From there it may be easily observed that: which may be solved for β by root-finding methods with level of precision set by number of terms used to approximate the infinite continued fraction.