Preprint
Article

On complex dynamics in a Suris's integrable map

Altmetrics

Downloads

80

Views

14

Comments

0

A peer-reviewed article of this preprint also exists.

This version is not peer-reviewed

Submitted:

29 March 2024

Posted:

01 April 2024

You are already at the latest version

Alerts
Abstract
Quantum tunneling in a two-dimensional integrable map is studied. The orbits of the map are all confined to the curves specified by the one-dimensional Hamiltonian. It is found that the behavior of tunneling splitting for the integrable map and the associated Hamiltonian system is qualitatively the same, with only a slight difference in magnitude. However, the tunneling tails of the wave functions, obtained by superposing the eigenfunctions that form the doublet, exhibit significant difference. To explore the origin of the difference, we observe the classical dynamics in the complex plane and find that the existence of branch points appearing in the potential function of the integrable map could play the role for yielding non-trivial behavior in the tunneling tail. The result highlights the subtlety of quantum tunneling, which cannot be captured in nature only by the dynamics in the real plane.
Keywords: 
Subject: Physical Sciences  -   Mathematical Physics

1. Introduction

Quantum tunneling is a subtle phenomenon because it cannot be captured by any power series, such as the WKB (Wentzel-Kramers-Brillouin) expansion, in terms of the Planck constant. This unique aspect can also be understood by recognizing that the power series expansion diverges, meaning that the tunneling effect is fundamentally non-perturbative [1].
Recent advances in resurgence theory have revealed that divergent series contain a wealth of information, even though they are divergent, and have provided a mathematical framework for dealing with exponentially small quantities [1,2,3,4,5]. Historically, the instanton method, proposed in the field theory and the theory of chemical reactions, is established when describing tunneling effects of the system with a single degree of freedom [6,7,8,9,10]. The instanton is a path running along the imaginary time, so it is not the solution in the real plane. Along the same lines, it has long been recognized that it is necessary to consider complex paths in order to treat exponentially small effects. In resurgence theory, complex paths, including the instanton path, can be formulated as singularities in the Borel plane, which is obtained by applying the Borel transform to given divergent series [1,2,3,4,5].
This paper discusses the subtlety of tunneling effects in a class of integrable maps for which one would expect nothing special since they are completely integrable. However, as shown below, this is not necessarily the case. The work presented here can be regarded as one of a series of studies on tunneling effects in nonintegrable systems based on the complex semiclassical analysis [11,12,13]. Since only energetic barrier tunneling but also dynamical tunneling [14,15,16] are exponentially small effects, the importance of complex paths is obvious and the use of complex orbits is unavoidable if one intends to extract components originating from tunneling effects. The integrable map studied here shares a certain characteristic with nonintegrable maps, although it is integrable.
To study tunneling effects in nonintegrable systems, discrete-time dynamical systems, commonly referred to as maps, are often used as model systems. This is supported by several reasons: among them, is that the behavior of complex paths has been well studied for discrete-time dynamical systems compared to continuous-time (Hamiltonian) dynamical systems [17,18,19]. However, we should keep in mind that it is not immediately clear whether the results derived for discrete-time dynamical systems, especially those based on complex dynamics, are applicable to continuous-time flow systems. The difficulty arises from the lack of a formal connection between the complex dynamics for the maps and continuous-time flow systems in which time is complexified.
On the other hand, one can find nonlinear integrable Hamiltonian flow systems, while nonlinear integrable maps are rather rare. The instanton is a path describing the tunneling transition in continuous Hamiltonian flow systems, and it is sometimes used as a reference when comparing the behavior of tunneling in integrable and nonintegrable systems [20]. In contrast, to the authors’ knowledge, there has not yet been an analysis of the tunneling effect for the map that is obtained by perturbing an integrable map. In most analyses, a continuous-time Hamiltonian system is taken as the integrable limit [21,22]. Such a treatment is somewhat self-inconsistent, and it is advisable to seek alternative solutions, if possible.
To fill this gap, here we study tunneling for a class of integrable maps found by Suris [23]. The integrability of the map is defined similarly to the case of continuous-time Hamiltonian systems. In particular, to verify the integrability for the two-dimensional symplectic map, it is sufficient to find a constant of motion under the discrete-time evolution. In Ref. [23], it was also shown that the maps with different types of potential function provide `Hamiltonians’, and it is easy to check that they are conserved under the time evolution.
Here we compare quantum tunneling for the integrable map and the associated `Hamiltonian’. We quantize the map as usual by introducing the one-step unitary operator, while Weyl quantization is applied to the continuous Hamiltonian. For a given classical Hamiltonian, recall that there are infinitely many quantizations resulting from the choice of the operator ordering. The differences are known to be of the order of 2 , which is larger than exponentially small effects. Thus, even within quantizations of a continuous Hamiltonian, each of which is determined by the chosen operator ordering, tunneling effects are in principle uncontrollable. It is therefore not surprising that the difference in quantization between integrable maps and the associated integrable Hamiltonian systems can lead to different tunneling properties.
Our strategy in this paper is to perform high-precision numerical calculations of quantum systems and to observe the behavior of the complex dynamics generated by both the integrable map and the corresponding Hamiltonian flow. We do not perform any semiclassical analysis here, but the comparison of the complex classical dynamics with the phase space profile of the time-evolved wave packet reveals non-trivial signatures in the tunneling components, even though the map is integrable.
Recently, the ultra-near integrable system has been studied in order to explore the difference of tunneling between integrable and nonintegrable systems [13,24]. Ultra-near integrable systems are defined as those systems in which the classical phase spaces do not exhibit any invariant structures inherent in nonintegrability in the size of the Planck cell. Even without islands of stability or chaotic seas, the tail of the wave function generates non-monotonic step structures [24], and the ergodicity of complex dynamics has been shown to play a key role in reproducing non-trivial tunneling behavior in nonintegrable maps [13]. Since there are no symptoms in the real phase space, the origin of the non-trivial tunneling tails comes down to the question of how chaos is born in the complex plane. The study of the integrable map is expected to be helpful in understanding this question.
The organization of this report is as follows. In Section 2, we introduce the symplectic map which was found to be integrable by Suris [23], together with the Hamiltonian to which the orbits generated by the map are confined. We then define the unitary operator that provides the quantum dynamics for the map. Section 3 gives numerical results showing the difference in tunneling signatures between the integrable and the corresponding Hamiltonian system. Although the tunneling splitting exhibits only a tiny difference, the tunneling tails for the localized states constructed from the eigenfunctions forming the doublet differ significantly. In Section 4 we observe the behavior of the classical dynamics in the complex plane. This is only a preliminary study towards the semiclassical understanding of the non-trivial tunneling revealed in the integrable map, but the result provides important implications for understanding the observed phenomenon. What distinguishes the two systems is the fact that the solution of the Hamiltonian flow has no singularities in the complex plane, while the derivative of the potential function has branch points, which lead to multivalued dynamics in the integrable map. We also discuss issues that are expected to arise in the development of semiclassical analysis in the complex domain. sec:conclusion summarizes the results obtained in the present report and gives some arguments, especially in relation to quantum tunneling in nonintegrable systems.

2. Model

To begin with, we introduce a two-dimensional area-preserving map f : T 2 T 2 ,
f : q k + 1 p k + 1 = q k + p k V ( q k ) p k V ( q k ) .
Suris [23] has shown that there exists a conserved quantity for the map:
H ( q k + 1 , p k + 1 ) = H ( q k , p k )
for some specific choices of the potential function V ( q ) . Among several possibilities, we here take a potential function defined by [25]
V λ ( q ) : = 1 π 0 q arctan λ sin 2 π q 1 + λ cos 2 π q d q .
Note that the potential function V λ ( q ) is periodic with period 1. For this potential, the Hamiltonian
H λ ( q , p ) : = cos 2 π p λ cos 2 π q + cos 2 π ( q p ) ,
satisfies the condition (2) for any λ . Below, we denote the map with the potential (3) by f λ to specify the map depends on the parameter λ .
Since the condition (2) holds for k Z , H λ can be regarded as a conserved quantity (Hamiltonian) for the map f λ . In addition, the Hamiltonian flow F λ ( q ( t ) , p ( t ) ) generated by the Hamiltonian equations of motion,
q ˙ ( t ) = H λ ( q , p ) p , p ˙ ( t ) = H λ ( q , p ) q ,
is obviously integrable. Therefore, the map f λ does not exhibit chaotic motion, and can be said to be integrable. In this report, we refer f λ as the Suris’s integrable map.
Figure 1(a) shows the functions shows the functions V λ ( q ) and V λ ( q ) for different values of λ 1. The profile of V λ ( q ) for | λ | < 1 looks a "skewed" sine function. As λ ± 1 , V λ ( q ) tends to a piecewise linear function, which is discontinuous at q = 1 / 2 for λ = 1 and q = 0 for λ = 1 . A fixed point ( q , p ) = ( 0 , 0 ) is elliptic-type if λ > 0 and hyperbolic if λ < 0 . For | λ | > 1 there are several discontinuous points in q [ 0 , 1 ) . In this report we take λ ( 1 , 0 ) .
Figure 1(b) and Figure 1(c) give phase space portraits for λ = 1 / 3 and 2 / 3 . Each dot and solid curve of the same color represents a single orbit k ( q 0 , p 0 ) k N starting from different initial conditions ( q 0 , p 0 ) and the corresponding contour curves of H λ ( q , p ) = H λ ( q 0 , p 0 ) , respectively. There are stable fixed points at ( q , p ) = ( ± 1 / 2 , 0 ) and an unstable fixed point at ( q , p ) = ( 0 , 0 ) . In addition, the map has stable periodic points with period 2 at ( q , p ) = ( arccos ( λ 2 ) / 2 π , ± arccos ( λ 2 2 1 ) / 2 π ) and unstable periodic points with period 2 at ( q , p ) = ( 0 , ± 1 / 2 ) (when we take mod 1 for q). Separatrices starting from ( q , p ; E ) = ( 0 , 0 ; 1 2 λ ) and ( q , p ; E ) = ( 0 , 1 / 2 ; 1 ) are shown in black dashed and black dash-dotted curve, respectively. As expected, the orbits f λ k ( q 0 , p 0 ) k N for any ( q 0 , p 0 ) R 2 lie on the contour curves of H λ ( q , p ) = H λ ( q 0 , p 0 ) .
We next introduce quantum dynamics associated with the classical map f λ . We here assume the periodic boundary condition for both q- and p-directions, which leads to the finite-dimensional Hilbert space whose dimension is given by, say, N, and apply the canonical quantization to the one-step time evolution, which gives the quantum map U ^ : | ψ K | ψ k + 1 [27,28],
U ^ λ : = e i p ^ 2 e i V λ ( q ^ ) .
Floquet theorem leads to the eigenvalue equation,
U ^ λ | ψ n = e i E n | ψ n .
Correspondingly, we introduce the quantized Hamiltonian H ^ λ ( q ^ p ^ ) by adopting the canonical quantization for the classical Hamiltonian (5). Note here that we take theWeyl’s ordering for the product of q and p. The eigenvalue equation for H ^ λ then reads
H ^ λ | Ψ n = E n | Ψ n .
In the following, the quantum number m of the quasi-eigenstate | ψ m is determined by the quantum number of the eigenstates | Ψ n attaining the maximal overlap | ψ m | Ψ n 2 .
In this paper, numerical computations for the quantum map are performed by using arbitrary precision arithmetics implemented in Mathematica, MATLAB with ADVANPIX toolbox [29].

3. Difference between | Ψ n and | ψ n

The discrete map f λ and the continuous Hamiltonian flow F λ share the same constant of motion H λ ( q , p ) , but the quantization procedures are different: the former is quantized by introducing the single-step unitary operator U ^ λ while the Weyl’s quantization is applied to the Hamiltonian H λ . Thus, it is by no means obvious whether the eigenvalues and eigenfunctions are identical and have the same properties.
In particular, we are interested here in tunneling effects, so the difference in quantization can be crucial. This is because, as mentioned above, exponentially small effects slip through ordinary semiclassical arguments, which are made by expressing related quantities in terms of the expansion of . Even within quantizations for the Hamiltonian H λ there is ambiguity about the operator ordering, leading to an uncertainty of the order of O ( 2 ) , so it is all the more unclear what will happen if the quantization method differs from each other.
Here we examine tunneling splittings generated in the nearly degenerate energy levels of the quantum map (6). To obtain a double-well-like setting, we impose the periodic boundary condition on the phase space region ( q , p ) = [ 1 , 1 ) × [ 1 / 2 , 1 / 2 ) . The periodicity of V ( q ) leads to the point-symmetric phase space around ( 0 , 0 ) (see Figure 1), and the `libration’ motion appears for 1 + 2 λ < q , p ) < 1 2 λ . Due to the symmetry with respect to the origin, the energy levels of the quantum map produce tunneling splittings, similar to the one-dimensional double-well potential system.
Figure 2(a) gives some eigenfunctions ψ n ( q ) and Ψ n ( q ) for = 1 / 200 π in semi-logarithmic scale, and Figure 2(b) shows the corresponding Husimi representation for | Ψ n in normal scale. There is no noticeable difference between ψ n ( q ) and Ψ n ( q ) , and the Husimi representation has support on the contour curves q , p ) = E n . Figure 3(a) and Figure 3(b) plot the tunneling splitting defined by
Δ E 0 = | E 1 E 0 | , Δ E 0 = E 1 E 0 ,
as a function of 1 / and λ , respectively. In contrast to the almost perfect agreement of the eigenfunctions shown in Figure 2(a), we can find a few orders of magnitude difference between the tunneling splittings Δ E 0 and Δ E 0 , while the parameter dependence of Δ E 0 and Δ E 0 on 1 / and λ is qualitatively the same; both are expected from the semiclassical analysis for integrable systems.
A tiny difference in the magnitude implies that some discrepancy should be hidden in the eigenfunctions. Figure 3(c) plots the magnitude 1 | ψ n | Ψ n | 2 as a function of λ , which shows that the difference between | ψ n and | Ψ n does indeed exist, although it is exponentially small, less than 10 8 for λ ( 1 , 0 ) . Note also that the difference 1 | ψ n | Ψ n | 2 increases as λ tends to 1.
The tunneling splitting appears as a result of the overlap of the wave function localized at the left ( q < 0 ) and the wave function localized at the right ( q > 0 ) :
| ψ L ( R ) : = 1 2 ( | ψ 1 ± | ψ 0 ) , | Ψ L ( R ) : = 1 2 ( | Ψ 1 ± | Ψ 0 ) .
The Figure 4(a) and Figure 4(b) present the wave functions Ψ L ( q ) and ψ L ( q ) and the corresponding Husimi representation, respectively. Let q × and q + × be the turning points of the manifold q , p ) = E 0 projected onto the q-axis. As seen in Figure 4(a), | Ψ L ( q ) | 2 has an exponentially decaying tunneling tail in the region q ( q × , q + × ) , as predicted by the WKB argument.
The Husimi representation provides more detailed information. As shown in Figure 4(b i), the amplitude of the Husimi representation for | Ψ L decays exponentially along the instanton curve The instanton curve is obtained here by the purely imaginary time evolution t i R of the Hamiltonian flow F λ starting from the turning point q × .
| ψ L ( q ) | 2 also shows an exponentially decaying tunneling tail from q = q × down to q = 0 . However, | ψ L ( q ) | 2 stops decaying just beyond q = 0 and then forms a plateau with an oscillatory pattern in the region q > 0 . As shown in Figure 4(b ii), the Husimi representation of | ψ L tells us that the localization along the separatrix starting at ( q , p ) = ( 0 , 0 ) is responsible for the oscillatory plateau observed in the plot of | ψ L ( q ) | 2 .
The plateau found in | ψ L ( q ) | 2 and the localized component along the separatrix are reminiscent of “tunneling across the separatrix”. Note that its importance for our understanding of the anomalous enhancement of the tunneling probability has been pointed out in Ref. [22]. However, they are not necessarily the same phenomenon. The main difference is that similar oscillation patterns in nonintegrable systems penetrate into the region q < 0 , while the plateau found in the Suris’s integrable map does not. To confirm this, we examine the interval [ q , q + ] in which | ψ L ( q ) | 2 is 10 times larger than | Ψ L ( q ) | 2 . Here, the edges q and q + are defined by
q : = min q [ 1 , 1 ) : log 10 ψ L ( q ) Ψ L ( q ) 2 > 1 ,
q + : = max q [ 1 , 1 ) : log 10 ψ L ( q ) Ψ L ( q ) 2 > 1 .
Figure 5 plots q as a function 1 / . Dots and a black dashed curve in fig:deviation indicate numerical data from q and a linear fit to the obtained data, respectively. The plot shows that the value of 1 / q increases linearly with 1 / , meaning that q 0 for 1 / . Note in particular that the decay rate of 1 / q with 1 / is extremely slow compared to nonintegrable cases [22].

4. Toward Semiclassical Understanding for Plateau with Oscillatory Pattern

While the WKB (smiclassical) analysis for eigenstates of integrable Hamiltonian systems is well established [30], there is no semiclassical prescription to analyze eigenstates of nonintegrable systems. On the other hand, the semiclassical analysis in the time domain is available even for nonintegrable systems, and the methodology incorporating the complex domain is now in a usable state [12,31,32,33,34]. Here, we discuss the origin of the plateau accompanied by oscillations observed in the quantum Suris’s integrable map by taking the time-domain semiclassical approach.
Among possible options for initial conditions, we choose the state | Ψ L , an eigenstate of the Hamiltonian H ^ λ , since we are particularly interested in the discrepancy between | ψ L and | Ψ L . The Figure 6(a) and Figure 6(b) illustrate the propagation of the wave function U ^ λ k | Ψ L in the q-representation and the Husimi representation, respectively. We can see that the amplitude of the wave function U ^ λ k | Ψ L for q [ q , q + ] increases rapidly with time, and after four steps it already reaches a profile similar to ψ L ( q ) . Note that the saturated amplitude of U ^ λ k | Ψ L in the interval q [ q , q + ] is slightly higher than the amplitude of ψ L ( q ) .
The Husimi representation shows us an explicit wave packet propagation in phase space. As time evolves, the dominant part of the wave function U ^ λ k | Ψ L passes through the fixed point ( q , p ) = ( 0 , 0 ) and propagates along the separatrix. After five steps, the amplitude spreads almost uniformly along the separatrix. This tells us that the plateau with oscillation observed in ψ L ( q ) appears as the result of the tunneling transport along the separatrix associated with the fixed point ( q , p ) = ( 0 , 0 ) .

4.1. Complex Dynamics of F λ and f λ

In this section we explore the classical dynamics F λ and f λ in the complex domain. To do this, we introduce notations for the real and imaginary parts as q= x + iy, and p = u +iv ,respectively.
First we observe the solution curves in the complex plane under the flow F λ . They lie on the equi-energy surface specified by the condition H λ = E n R . To obtain these curves, we first find the turning points q ± × , marked by the black "×" in Figure 7, and then solve the equations of motion for the Hamiltonian H λ along the imaginary time axis i R to connect the turning points q ± × R . Recall that this solution curve, shown by the thick black curve, is nothing more than the instanton. Along the instanton curve we place equally spaced initial conditions, from each of which we solve the equations of motion along the real time axis R , yielding a set of closed curves shown in green.
As for the dynamics of the map H λ in the complex plane, the multivaluedness of the potential function, V λ generates non-trivial behavior. To see this, we focus on the dynamics at the points ( x , u ) = ( 0 , 0 ) , ( 0 , ± 1 / 2 ) , and (±1/2, ±1/2), respectively. The equations of motion of the Hamiltonian flow F λ are written explicitly as
y ˙ ( t ) = 2 π α 1 sinh 2 π v + α 2 λ sinh 2 π ( y v ) , v ˙ ( t ) = 2 π λ α 3 sinh 2 π y + α 4 sinh 2 π ( y v ) , x ˙ ( t ) = u ˙ ( t ) = 0 ,
where (α1234)=(1,-1,1,1) for (x,u)=(0,0),(α1234)=(-1,-1,-1,1) for (x, u)=(±1/2,±1/2), and (α1234)=(-1,1,1,-1) for (x, u)=(0,±1/2). Since x ( t ) , u ( t ) = const . , the solution curves shown in Equation (13) are restricted to the purely imaginary (y,v)-plane. Note that the fixed points (y, v) = (0, 0) of Equation (13) are all hyperbolic, and the right-hand sides of Equation (13) consist of monotonic and entire functions, so that the solution curves around (y, v)=(0,0) diverge monotonically to y(t), v(t) → ±∞ under forward time evolution (see Figure 8).
The greens dots in Figure 7 show the classical orbits k ( q 0 , p 0 ) k Z for which the initial points ( q 0 , p 0 ) are set along the instanton. We observe that all the orbits k ( q 0 , p 0 ) k Z preserve the energy E as expected, but there are exceptions: the orbits starting at ( x 0 , u 0 ) = ( 0 , 0 ) are not bounded on the Hamiltonian flow F λ . As explained below, the multivaluedness of the function V λ ( q ) for q C leads to such anomalous behavior. In Figure 7(b), the cyan dots show the projection of the itinerary of the orbit whose initial condition is placed on the instanton orbit onto the ( x , u ) -plane. We can see that the real part ( x k , u k ) of the orbit jumps among ( x k , v k ) = ( 0 , 0 ) , ( ± 1 / 2 , ± 1 / 2 ) , and ( 0 , ± 1 / 2 ) , while the real part ( x ( t ) , u ( t ) ) does not change under the Hamiltonian flow F λ [see eq:Heq].
This itinerary of the orbit can be understood by considering the derivative V λ ( q ) of the potential function in the complex q-plane. Introducing the new variable,
z = λ sin 2 π q 1 + λ cos 2 π q ,
we find that the derivative V λ ( z ) can be expressed as
V λ ( z ) = 1 π arctan z ,
where the arctangent function has an alternative expression,
arctan z = 1 2 i log 1 + i z 1 i z = d z 1 + z 2 ,
which obviously has the branch points on the z-plane at z = ± i , originating from the logarithmic function. Correspondingly, on the q-plane V λ ( q ) has the branch points at q = ± i 1 2 π log ( λ ) = : ± i y .
Now, we consider the itinerary of the orbit starting from ( x , y , u , v ) = ( 0 , y 0 , 0 , v 0 ) . In other words, we set the real part of the initial points to ( x 0 , u 0 ) = ( 0 , 0 ) , but take various values for the imaginary parts ( y , v ) = ( y 0 , v 0 ) , assuming that | y 0 | < y . Since the V λ ( q ) is a multivalued function in the logarithmic type, the one-step iteration leads to an infinite number of images. Several studies have explored the complex dynamics associated with multivalued functions [35], but here we consider only the one-step iteration to get rid of the complication caused by the multi-step iteration.
Since the stability at the origin ( y , v ) = ( 0 , 0 ) is hyperbolic [see fig:Hflow(a)], for any given v 0 , there exists an initial condition | y 0 | < y such that the next step satisfies the condition | y 1 | > y . Recall that the Hamiltonian flow follows eq:Heq, so in the time evolution the x ( t ) remains constant ( x ( t ) = 0 ) and only the value of y ( t ) changes. Therefore, in the ( x , y ) -plane, the point ( 0 , y 0 ) passes over the branch points ( 0 , ± y ) before and after being mapped to the point ( 0 , y 1 ) . After passing one of the branch points, the mapped point gains either 1 / 2 + m and 1 / 2 m ( m , m Z ) , reflecting the multivalued nature of the potential function (11). As a result, the initial point ( x 0 , u 0 ) = ( 0 , 0 ) is mapped to the lattice points ( 1 / 2 + m , 1 / 2 + m ) ( m , m Z ) . However, these lattice points can be identified as ( 1 / 2 , ± 1 / 2 ) and ( 1 / 2 , ± 1 / 2 ) by applying the periodic boundary condition. So the real parts ( x 0 , u 0 ) = ( 0 , 0 ) are mapped to ( x 1 , u 1 ) = ( 1 / 2 , ± 1 / 2 ) and ( 1 / 2 , ± 1 / 2 ) .
For the initial points ( 0 , y 0 ) in the ( x , y ) -plane that do not pass the branch points within the single step iteration, they simply move along the solution curve governed by eq:Heq. The points generated by the multivalued nature of the map are again identified by applying the periodic boundary condition. Depending on the initial condition ( x 0 , u 0 ) , the single step iteration in the ( y , v ) -plane gives the map from ( y 0 , v 0 ) to ( y 1 , v 1 ) on the same flow curve. Note that the transition between different ( x , u ) -planes does not occur because we here assume that the initial points ( 0 , y 0 ) in the ( x , y ) -plane do not pass the branch points.

4.2. One-step propagation

As illustrated in Figure 6, the plateau with oscillation observed in the eigenfunction | ψ n is reproduced by a wave function with relatively short time steps. This should make it possible to develop the semiclassical analysis in the time domain.
As explained in the previous section, the multivaluedness of the map leads to the transition that never occurs in the integrable Hamiltonian flow. The orbits of the Hamiltonian flow F λ starting from the unstable fixed point ( x ( 0 ) , u ( 0 ) ) = ( 0 , 0 ) do not change the real part, i.e., ( x ( t ) , u ( t ) ) = ( 0 , 0 ) ( t > 0 ) independent of the initial conditions ( y 0 , v 0 ) , while the orbits generated by the map f λ jump from (x0,u0) = (0,0) to either (x1,u1) = (1/2, ±1/2) or (−1/2, ±1/2) when they pass over the branch points, otherwise they stay at (x1, u1)=(0,0) Such a transition should explain the propagation of the wave function along the separatrix, but some issues need to be clarified before proceeding.
The first one is that the orbit in the complex plane is bifurcated even in one-step propagation, reflecting the multivaluedness of the map. The transition from ( 0 , 0 ) to ( 1 / 2 , 1 / 2 ) is observed in the one-step propagation of the quantum wave packet [see Figure 10], while the orbit can jump from ( 0 , 0 ) to ( ± 1 / 2 , 1 / 2 ) as explained in the previous section. We expect that this can be solved by properly treating the Stokes phenomenon when performing the saddle point approximation. If the path from ( 0 , 0 ) to ( ± 1 / 2 , 1 / 2 ) gives rise to a divergent solution, then it should be removed from the final solution due to the Stokes phenomenon. For multi-step analysis, however, the multivaluedness of the map must be treated in each step, which makes semiclassical analysis difficult. As can be seen from Figure 6, the plateau along the separatrix becomes more visible as time evolves. It may be a relatively short time, but we have to perform a semiclassical calculation with multiple time steps. This requires a more sophisticated approach to the Stokes phenomenon.
On the other hand, Figure 10 shows that the plateau can also be made visible by increasing the parameter λ . This may give us the impression that the one-step semiclassical analysis is sufficient to explain the observed phenomenon, and one can get rid of multi-step calculations. However, the problem is not simple, because as the parameter λ is increased, the potential V λ ( q ) and its derivative V λ ( q ) become sharper, as shown in Figure 1(a). This reminds us of the diffraction effect [36], which is expected to occur when the potential function has discontinuities. We should include higher-order terms in the saddle approximations, or more legitimately develop uniform approximation type arguments to cope with such situations.

5. Conclusion and Discussion

In this report, we have studied quantum tunneling for a nonlinear integrable map found by Suris, with a particular focus on the comparison with the corresponding continuous Hamiltonian system. Although the Suris’s map is completely integrable in the sense of Liouville-Aronld in Hamiltonian systems, and the orbits in the real plane generated by the Suris’s map lie on the equi-energy curve for the continuous Hamiltonian, our numerical computations performed with arbitrary precision arithmetic have revealed that the characteristics of the tunneling tails in the eigenfunctions differ from each other. This means that the behavior of the tunneling tail cannot be deduced from the classical dynamics in real phase space.
The discrepancy found here between the map and the corresponding Hamiltonian system is not surprising. As noted in the introduction, even for the same Hamiltonian, there is an ambiguity in the quantization arising from the choice of the operator ordering, while the magnitude of the tunneling effects is in general exponentially small with respect to the Planck constant , which cannot be captured by any power series expansion.
A lesson we learn from this example is that it is necessary to investigate the dynamics in the complex plane to understand the tunneling effect properly. Here, we have not performed a fully semiclassical calculation with complex orbits, but it would be reasonable to expect that the difference in the behavior of the complex dynamics gives rise to the difference in the tunneling tail. The branch points, appearing in the potential function of the map, produce multivaluedness of the complex map whose treatment is one of the topics in the theory of complex dynamical systems [35]. We have shown that the real part of the orbits at the fixed point of the Hamiltonian flow remains the same under the dynamics, while the orbit changes its real part discontinuously before and after crossing the branch point. This must be the origin of the localization of the wave function along the separatrix for the fixed point, and also reminiscent of the tunneling transition across the separatrix found in Refs. [21,22].
The map f λ is slightly complicated compared to the flow F λ due to the presence of singularities producing multivaluedness of the dynamics. This fact reminds us of the so-called Poinlev`e conjecture: in integrable systems, singularities appearing in the time plane are at most limited to the poles [37]. This in turn suggests that the system can be nonintegrable if branch points appear in the complex time plane of the solutions. One can say that the Suris’s integrable map is completely integrable as far as the dynamics in the real plane is concerned, whereas it is one step outside of integrable systems. Our present result shows that this difference is reflected in the tunneling effect.
It is important to note that such subtlety is not limited to the system studied in this report. If the system is close enough to the integrable limit, the phase space is almost covered by Kolmogorov-Arnold-Moser curves. For the classical-to-quantum correspondence, the size of the Planck cell is an important spatial scale; even if invariant structures such as nonlinear resonances or stochastic layers appear in phase space, it does not matter if their sizes are smaller than the size of the Planck cell under consideration. This is a fundamental working hypothesis when considering the manifestation of classical chaos in the corresponding systems [38]. According to this principle, the nearly integrable system, in which the invariant structures inherent in the nonintegrability are all on the sub-Planck cell scale, should be identified with the completely integrable system.
This must hold true in the real plane, but not necessarily in the complex plane. This is because the phase space profile for the nonintegrable system, no matter how close it is to an integrable system, is dramatically different from that of the completely integrable system. The KAM curves are analytic in the real plane, but they have natural boundaries in the complex plane [39,40,41,42]. Furthermore, it is highly probable that homoclinic intersections between stable and unstable manifolds occur beyond the natural boundaries [12]. Chaos in the complex plane developed there is difficult to `see’ from the phase space profile in the real plane, but tunneling transport is driven by chaotic orbits in such a deep complex plane [13].
Quantum tunneling is a phenomenon that does not exist in classical mechanics. There are approaches that try to understand quantum tunneling in terms of invariant objects in the real phase space [43,44]. These implicitly assume that tunneling tails of wave functions in nonintegrable systems can be expressed in terms of tunneling tails of wave functions supported by invariant manifolds in the real phase space. In other words, the broadening of the wave functions [45] is assumed to be sufficient to capture the signature of nonintegrable tunneling. However, there is a situation in which such an approach does not work [13], and there one needs to make full use of chaos in the complex plane. Chaos in the complex plane appears beyond natural boundaries and thus is not reachable from the real plane. We believe that this fact is not necessarily limited to the situation where no visible invariant structures appear compared to the Planck cell as reported in Ref. [13], but the mechanism found there works in tunneling processes in nonintegrable systems in general. The phase spaces used so far to study the tunneling in nonintegrable systems are too complex to specify purely quantum phenomena, which could slip through any kind of analysis based on classical dynamics in the real plane. The question of what are the unique characteristics of tunneling effects in nonintegrable systems should be revisited.

Author Contributions

Conceptualization, Y.H. and A.S.; Investigation, Y.H.; writing—original draft, Y.H. and A.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by JSPS KAKENHI Grants No.17K05583 and 22H01146.

Data Availability Statement

Not applicable.

Acknowledgments

The authors are grateful to Yutaka Ishii for helpful comments on the complex dynamical systems.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Voros, A. The return of the quartic oscillator. The complex WKB method. Annales de l’I.H.P. Physique théorique 1983, 39, 211–338. [Google Scholar]
  2. Dingle, R. Asymptotic expansions: their derivation and interpretation; Vol. 521, Academic Press London, 1973.
  3. Écalle, J. Les fonctions résurgentes:(en trois parties); Vol. 1, Université de Paris-Sud, Département de Mathématique, Bât. 425, 1981.
  4. Delabaere, E.; Dillinger, H.; Pham, F. Exact semiclassical expansions for one-dimensional quantum oscillators. Journal of Mathematical Physics 1997, 38, 6126–6184. [Google Scholar] [CrossRef]
  5. Kawai, T.; Takei, Y. Algebraic analysis of singular perturbation theory; Vol. 227, American Mathematical Society, 2005.
  6. Callan Jr, C.G.; Coleman, S. Fate of the false vacuum. II. First quantum corrections. Physical Review D 1977, 16, 1762. [Google Scholar] [CrossRef]
  7. Coleman, S. Fate of the false vacuum: Semiclassical theory. Physical Review D 1977, 15, 2929. [Google Scholar] [CrossRef]
  8. Coleman, S. Aspects of Symmetry: Selected Erice Lectures; Cambridge University Press, 1988.
  9. George, T.F.; Miller, W.H. Complex-Valued Classical Trajectories for Linear Reactive Collisions of H + H2 below the Classical Threshold. The Journal of Chemical Physics 1972, 56, 5722–5723. [Google Scholar] [CrossRef]
  10. Miller, W.H. Quantum mechanical transition state theory and a new semiclassical model for reaction rate constants. The Journal of Chemical Physics 1974, 61, 1823–1834. [Google Scholar] [CrossRef]
  11. Shudo, A.; Ikeda, K.S. Complex classical trajectories and chaotic tunneling. Physical Review Letters 1995, 74, 682. [Google Scholar] [CrossRef] [PubMed]
  12. Shudo, A.; Ikeda, K.S. Complex Semiclassical Approach to Chaotic Tunneling. Dynamical Tunneling: Theory and Experiment.
  13. Koda, R.; Hanada, Y.; Shudo, A. Ergodicity of complex dynamics and quantum tunneling in nonintegrable systems. Physical Review E 2023, 108, 054219. [Google Scholar] [CrossRef] [PubMed]
  14. Davis, M.J.; Heller, E.J. Quantum dynamical tunneling in bound states. The Journal of Chemical Physics 1981, 75, 246–254. [Google Scholar] [CrossRef]
  15. Creagh, S.C. Tunneling in two dimensions. In Tunneling in Complex Systems; Tomsovic, S., Ed.; World Scientific, 1998; p. 35.
  16. Keshavamurthy, S.; Schlagheck, P. Dynamical tunneling: theory and experiment; CRC Press, 2011.
  17. Milnor, J. Dynamics in One Complex Variable; Vol. 160, Princeton University Press, 2011.
  18. Beardon, A.F. Iteration of rational functions: Complex analytic dynamical systems; Vol. 132, Springer Science & Business Media, 2000.
  19. Morosawa, S.; Nishimura, Y.; Taniguchi, M.; Ueda, T. Holomorphic dynamics; Vol. 66, Cambridge University Press, 2000.
  20. Takahashi, K. Semiclassical Analysis of Multidimensional Barrier Tunneling. Dynamical Tunneling: Theory and Experiment.
  21. Hanada, Y.; Shudo, A.; Ikeda, K.S. Origin of the enhancement of tunneling probability in the nearly integrable system. Physical Review E 2015, 91, 042913. [Google Scholar] [CrossRef] [PubMed]
  22. Hanada, Y.; Ikeda, K.S.; Shudo, A. Dynamical tunneling across the separatrix. Physical Review E 2023, 108, 064210. [Google Scholar] [CrossRef] [PubMed]
  23. Suris, Y.B. Integrable mappings of the standard type. Functional Analysis and Its Applications 1989, 23, 74–76. [Google Scholar] [CrossRef]
  24. Iijima, R.; Koda, R.; Hanada, Y.; Shudo, A. Quantum tunneling in ultra-near-integrable systems. Physical Review E 2022, 106, 064205. [Google Scholar] [CrossRef]
  25. Meiss, J.D. Symplectic maps, variational principles, and transport. Reviews of Modern Physics 1992, 64, 795–848. [Google Scholar] [CrossRef]
  26. Takahasi, H.; Mori, M. Double Exponential Formulas for Numerical Integration. Publications of the Research Institute for Mathematical Sciences 1973, 9, 721–741. [Google Scholar] [CrossRef]
  27. Berry, M.V.; Balazs, N.L.; Tabor, M.; Voros, A. Quantum maps. Annals of Physics 1979, 122, 26–63. [Google Scholar] [CrossRef]
  28. Casati, G.; Chirikov, B.V.; Izraelev, F.M.; Ford, J. Stochastic behavior of a quantum pendulum under a periodic perturbation. In Proceedings of the Stochastic Behavior in Classical and Quantum Hamiltonian Systems; Casati, G.; Ford, J., Eds., Berlin, Heidelberg, 1979; Lecture Notes in Physics; pp. 334–352. [Google Scholar] [CrossRef]
  29. Advanpix LLC.. Multiprecision Computing Toolbox for MATLAB. http://www.advanpix.com/.
  30. Berry, M.V.; Mount, K.E. Semiclassical approximations in wave mechanics. Reports on Progress in Physics 1972, 35, 315. [Google Scholar] [CrossRef]
  31. Shudo, A.; Ishii, Y.; Ikeda, K.S. Julia sets and chaotic tunneling: I. Journal of Physics A 2009, 42, 265101. [Google Scholar] [CrossRef]
  32. Shudo, A.; Ishii, Y.; Ikeda, K.S. Julia sets and chaotic tunneling: II. Journal of Physics A 2009, 42, 265102. [Google Scholar] [CrossRef]
  33. Shudo, A.; Ikeda, K.S. Stokes geometry for the quantum Hénon map. Nonlinearity 2008, 21, 1831. [Google Scholar] [CrossRef]
  34. Shudo, A.; Ikeda, K.S. Toward pruning theory of the Stokes geometry for the quantum Hénon map. Nonlinearity 2016, 29, 375. [Google Scholar] [CrossRef]
  35. Bullett, S.; Lomonaco, L.; Siqueira, C. Correspondences in complex dynamics. New Trends in One-Dimensional Dynamics: In Honour of Welington de Melo on the Occasion of His 70th Birthday IMPA 2016, Rio de Janeiro, Brazil, November 14–17, 14 November. [CrossRef]
  36. Ishikawa, A.; Tanaka, A.; Ikeda, K.S.; Shudo, A. Diffraction and tunneling in systems with mixed phase space. Physical Review E 2012, 86, 036208. [Google Scholar] [CrossRef] [PubMed]
  37. Ablowitz, M.J.; Ramani, A.; Segur, H. A connection between nonlinear evolution equations and ordinary differential equations of P-type. II. Journal of Mathematical Physics 1980, 21, 1006–1015. [Google Scholar] [CrossRef]
  38. Haake, F. Quantum signatures of chaos; Springer, 1991.
  39. Greene, J.M.; Percival, I.C. Hamiltonian maps in the complex plane. Physica D 1981, 3, 530–548. [Google Scholar] [CrossRef]
  40. Percival, I.C. Chaotic boundary of a Hamiltonial map. Physica D 1982, 6, 67–77. [Google Scholar] [CrossRef]
  41. Berretti, A.; Chierchia, L. On the complex analytic structure of the golden invariant curve for the standard map. Nonlinearity 1990, 3, 39. [Google Scholar] [CrossRef]
  42. Berretti, A.; Marmi, S. Standard map at complex rotation numbers: creation of natural boundaries. Physical Review Letters 1992, 68, 1443. [Google Scholar] [CrossRef] [PubMed]
  43. Tomsovic, S.; Ullmo, D. Chaos-assisted tunneling. Physical Review E 1994, 50, 145. [Google Scholar] [CrossRef] [PubMed]
  44. Brodier, O.; Schlagheck, P.; Ullmo, D. Resonance-assisted tunneling. Annals of Physics 2002, 300, 88–136. [Google Scholar] [CrossRef]
  45. Balazs, N.; Voros, A. Wigner’s function and tunneling. Annals of Physics 1990, 199, 123–140. [Google Scholar] [CrossRef]
1
The integration of V λ ( q ) is numerically evaluated by "double exponential quadrature" [26] implemented in Mathematica until numerical convergence is achieved to at least with 200 digits of mantissa.
Figure 1. (a) Plot of the functions V λ ( q ) (upper) and V λ ( q ) (lower) for λ = 1 / 3 , 2 / 3 , and 1 . Phase space portrait for (b) λ = 1 / 3 and (c) λ = 2 / 3 . In (b) and (c), dots and solid curves represent the orbits for f λ and the associated contour curve H λ ( q , p ) , respectively. Black dashed and black dash-dotted curves represent the separatrix starting from (q,p)=(0,0) and (0, ±1/ 2), respectively.
Figure 1. (a) Plot of the functions V λ ( q ) (upper) and V λ ( q ) (lower) for λ = 1 / 3 , 2 / 3 , and 1 . Phase space portrait for (b) λ = 1 / 3 and (c) λ = 2 / 3 . In (b) and (c), dots and solid curves represent the orbits for f λ and the associated contour curve H λ ( q , p ) , respectively. Black dashed and black dash-dotted curves represent the separatrix starting from (q,p)=(0,0) and (0, ±1/ 2), respectively.
Preprints 102573 g001
Figure 2. (a) Eigenstates ψ n ( q ) (cyan) and Ψ n ( q ) (black) with = 1 / 200 π for (i) n = 0 , (ii) n = 50 , (iii) n = 150 , and (iv) n = 199 . (b) Husimi representation of ψ n ( q ) in normal scale for (i) n = 0 , (ii) n = 50 , (iii) n = 150 , and (iv) n = 199 . The black dotted curve shows the contour curve with q , p ) = E n . The black broken curve shows the separatrix starting from ( q , p ) = ( 0 , 0 ) .
Figure 2. (a) Eigenstates ψ n ( q ) (cyan) and Ψ n ( q ) (black) with = 1 / 200 π for (i) n = 0 , (ii) n = 50 , (iii) n = 150 , and (iv) n = 199 . (b) Husimi representation of ψ n ( q ) in normal scale for (i) n = 0 , (ii) n = 50 , (iii) n = 150 , and (iv) n = 199 . The black dotted curve shows the contour curve with q , p ) = E n . The black broken curve shows the separatrix starting from ( q , p ) = ( 0 , 0 ) .
Preprints 102573 g002
Figure 3. (a) Plot of the tunneling splittings Δ E 0 (solid) and Δ E 0 (dashed) vs of 1 / . (b) Plot of the tunneling splittings Δ E 0 (solid) and Δ E 0 (dashed) vs λ . (c) Plot of 1 | ψ 0 | Ψ 0 | vs λ .
Figure 3. (a) Plot of the tunneling splittings Δ E 0 (solid) and Δ E 0 (dashed) vs of 1 / . (b) Plot of the tunneling splittings Δ E 0 (solid) and Δ E 0 (dashed) vs λ . (c) Plot of 1 | ψ 0 | Ψ 0 | vs λ .
Preprints 102573 g003
Figure 4. (a) Wave functions ψ L ( q ) (solid curves) and Ψ L ( q ) (dashed curves) for different values of 1 / . The "" mark is put to indicate the deviation point q . (b) Husimi representation of (i) Ψ L ( q ) and (ii) ψ L ( q ) for 1 / = 200 π in the log 10 scale. The dashed gray curve represents the separatrix. Closed cyan curves show the contour curve associated with the ground state energy level E 0 = 1.639 057 02 . Black and white "×" symbols indicate the position of the turning points for q ( t ) ; ( q , p ) = ( ± 0.450 847 , ± 0.012 215 ) , respectively. (b i) The cyan curve connected to the two turning points shows the instanton curve projected onto the ( x , u ) -plane.
Figure 4. (a) Wave functions ψ L ( q ) (solid curves) and Ψ L ( q ) (dashed curves) for different values of 1 / . The "" mark is put to indicate the deviation point q . (b) Husimi representation of (i) Ψ L ( q ) and (ii) ψ L ( q ) for 1 / = 200 π in the log 10 scale. The dashed gray curve represents the separatrix. Closed cyan curves show the contour curve associated with the ground state energy level E 0 = 1.639 057 02 . Black and white "×" symbols indicate the position of the turning points for q ( t ) ; ( q , p ) = ( ± 0.450 847 , ± 0.012 215 ) , respectively. (b i) The cyan curve connected to the two turning points shows the instanton curve projected onto the ( x , u ) -plane.
Preprints 102573 g004
Figure 5. The dots represent the deviation point q as a function of 1 / in case of λ = 1 / 3 . The dashed line is obtained by applying linear regression to the plot of 1 / q vs 1 / .
Figure 5. The dots represent the deviation point q as a function of 1 / in case of λ = 1 / 3 . The dashed line is obtained by applying linear regression to the plot of 1 / q vs 1 / .
Preprints 102573 g005
Figure 6. Plot of U ^ λ k | Ψ L (a) in the q-representation and (b) in the Husimi representation.
Figure 6. Plot of U ^ λ k | Ψ L (a) in the q-representation and (b) in the Husimi representation.
Preprints 102573 g006
Figure 7. The solution curves obtained under the continuous Hamiltonian flow F λ . They are confined in the equi-energy surface given by the condition H λ = E 0 = 1.639 057 02 . The curves are projected onto (a) the ( x , u , v ) -space and (b) the ( x , u ) -plane. In (a) and (b) the black closed curve and the black semicircular curve represent the equi-energy curves in the real plane and the instanton curve, respectively. The symbol "×" stands for the turning point of x ( t ) . (a) The cyan line tending to infinity shows a solution curve of H λ starting from the point ( x , u ) = ( 0 , 0 ) , marked by the cyan point on the instanton. (b) The cyan dots show the orbits of f λ starting from the point ( x , u ) = ( 0 , 0 ) on the instanton.
Figure 7. The solution curves obtained under the continuous Hamiltonian flow F λ . They are confined in the equi-energy surface given by the condition H λ = E 0 = 1.639 057 02 . The curves are projected onto (a) the ( x , u , v ) -space and (b) the ( x , u ) -plane. In (a) and (b) the black closed curve and the black semicircular curve represent the equi-energy curves in the real plane and the instanton curve, respectively. The symbol "×" stands for the turning point of x ( t ) . (a) The cyan line tending to infinity shows a solution curve of H λ starting from the point ( x , u ) = ( 0 , 0 ) , marked by the cyan point on the instanton. (b) The cyan dots show the orbits of f λ starting from the point ( x , u ) = ( 0 , 0 ) on the instanton.
Preprints 102573 g007
Figure 8. Flow of Equation (13) on the ( y , v ) -plane for λ = 1 / 3 . Vertical dashed lines represent y = y for λ = 1 / 3 .
Figure 8. Flow of Equation (13) on the ( y , v ) -plane for λ = 1 / 3 . Vertical dashed lines represent y = y for λ = 1 / 3 .
Preprints 102573 g008
Figure 9. (a) The function V λ ( q ) in the ( x , y ) -plane with λ = 1 / 3 . The argument of V λ ( q ) and the absolute value | V λ ( q ) | are distinguished by color and brightness, respectively. The blue diamond "⋄" symbol represents the branch point of V ( q ) . The white line starting from the branch point represents the branch cut. The black curve shows a contour curve Im ( V λ ( q ) ) = 0 .
Figure 9. (a) The function V λ ( q ) in the ( x , y ) -plane with λ = 1 / 3 . The argument of V λ ( q ) and the absolute value | V λ ( q ) | are distinguished by color and brightness, respectively. The blue diamond "⋄" symbol represents the branch point of V ( q ) . The white line starting from the branch point represents the branch cut. The black curve shows a contour curve Im ( V λ ( q ) ) = 0 .
Preprints 102573 g009
Figure 10. (a) Plot of the wave functions U ^ λ | Ψ L (red curve) and | Ψ L (black dashed curve) for (i) λ = 0.3 , (ii) λ = 0.5 , (iii) λ = 0.7 , and (iv) λ = 0.9 . (b) Husimi representation of U ^ λ | Ψ L . The cyan dotted line represents a line ( q , p ) = ( 0 , 0 ) to ( q , p ) = ( 1 / 2 , 1 / 2 ) for a guide. The gray solid curve and the gray dashed curve represent a contour curve of q , p ) = E 0 and the separatrix starting from ( q , p ) = ( 0 , 0 ) , respectively.
Figure 10. (a) Plot of the wave functions U ^ λ | Ψ L (red curve) and | Ψ L (black dashed curve) for (i) λ = 0.3 , (ii) λ = 0.5 , (iii) λ = 0.7 , and (iv) λ = 0.9 . (b) Husimi representation of U ^ λ | Ψ L . The cyan dotted line represents a line ( q , p ) = ( 0 , 0 ) to ( q , p ) = ( 1 / 2 , 1 / 2 ) for a guide. The gray solid curve and the gray dashed curve represent a contour curve of q , p ) = E 0 and the separatrix starting from ( q , p ) = ( 0 , 0 ) , respectively.
Preprints 102573 g010
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

Subscribe

© 2024 MDPI (Basel, Switzerland) unless otherwise stated