Quantum tunnelling, real-time dynamics and Picard-Lefschetz thimbles

We follow up the work, where in light of the Picard-Lefschetz thimble approach, we split up the real-time path integral into two parts: the initial density matrix part which can be represented via an ensemble of initial conditions, and the dynamic part of the path integral which corresponds to the integration over field variables at all later times. This turns the path integral into a two-stage problem where, for each initial condition, there exits one and only one critical point and hence a single thimble in the complex space, whose existence and uniqueness are guaranteed by the characteristics of the initial value problem. In this paper, we test the method for a fully quantum mechanical phenomenon, quantum tunnelling in quantum mechanics. We compare the method to solving the Schr\"odinger equation numerically, and to the classical-statistical approximation, which emerges naturally in a well-defined limit. We find that the Picard-Lefschetz result matches the expectation from quantum mechanics and that, for this application, the classical-statistical approximation does not.


Introduction
The computation of real-time correlators in quantum mechanics may be formally phrased through a path integral formulation of the time evolution of the density matrix. Given that a correlator is expressed through (1.1) one may rewrite this as a path integral (1. 2) The operatorÔ is assumed to be expressed in terms of some fields φ, that now live on the Schwinger-Keldysh contour C [1, 2] in the complex time-plane. This means that their time index starts at the initial time t = 0, runs to some late time, and returns to t = 0 (see below). ρ 0 is a representation of the initial density matrix, L is the Lagrangian expressed in terms of the fields φ and Z is the partition function. This integral can in certain simple cases be computed explicitly, and advanced perturbative techniques allows the computation of real-time correlators, say in equilibrium. However, in the cases where non-perturbative or out-of-equilibrium information is of interest, the situation becomes very challenging. Numerical simulations must then often be employed, and since the integrand is a highly oscillatory complex phase (the Lagrangian L is usually real), convergence of numerical integration is in general poor due to this "sign problem".
Hence, to the extent that the initial state permits it 1 , all the thimbles will be included in the total path integral by sampling the ensemble of initial conditions one by one. In [35] we implemented this in quantum mechanics (0-dimensional field theory) for a simple scalar model, and computed the non-equal time correlator.
Still, it would be interesting to stress-test our method for a truly quantum mechanical phenomenon, in a non-trivial system. And so in the present paper we study quantum mechanical tunnelling (for applications of Lefschetz thimble to tunnelling problem from a more theoretical viewpoint, see [36][37][38]), precisely because we expect this process to challenge the numerical path integral. Furthermore, we expect the difference between the thimble approach and the classical-statistical approximation to be much more pronounced. In our real-time thimble approach, we sample the initial conditions and solve for the classical trajectories. The further identification and sampling along the thimbles provides the quantum contribution, but omitting this step allows us to recover the standard classical-statistical averaging. This allows a straightforward comparison between the two methods.
The structure of the paper is as follows: In section 2, we present the closed-time path integral in terms of the initial density matrix and thimble. In section 3, we solve the textbook quantum mechanical tunnelling problem of ground state tunnel splitting, utilizing different methods. Finally, section 4 is our summary and conclusions.

Setting up the path integral
The section is devoted to a short, but hopefully sufficiently self-contained, description of the closed-time path integral and thimble methods. For a thorough discussion of the topics, we refer the reader to [35]. Since quantum field theory problems are what we ultimately want to tackle, we keep a spatial index in this section. We have in mind a single real scalar field with a double-well potential, which we will use later in the next section. All the conclusions drawn in the present section, however, hold true for arbitrary potential.

The closed-time path integral
As described briefly above, the Schwinger-Keldysh closed-time path integral [1, 2] appears naturally when we evaluate correlation functions at time t, given an initial density matrix at t 0 , As shown in Fig. 1, the path comprises two branches starting at some initial time t 0 (typically, t 0 = 0), and ending at some final time t m . The upper branch hosts fields φ + j , while the lower one hosts φ − j . In the following, we will think of a finite set of discrete values of time t separated by dt, as appropriate for the numerical simulations to be performed below. Formally, the continuum limit arises from taking dt → 0. These branches appear when one inserts complete sets of states 2 , |φ j , t j , into (1.1), located at times from t 0 to some time t m , later than any time appearing in the operatorÔ. One then finds a number of inner products of the form φ i , t i |φ j , t j , which may be evaluated by using the Hamiltonian to evolve each state to some intermediate time t c . In [35] we took t c to lie at the midpoint between t i and t j . Besides the choice of midpoint, one can also choose the forward or backward difference, for instance, to select t j or t i to be t c . But this comes with a warning. For instance, if one chooses the forward difference on the upper branch, i.e., evaluating φ + j+1 , t j+1 |φ + j , t j at t j+1 , then one should also evaluate the lower branch amplitude φ − j , t j |φ − j+1 , t j+1 at t j+1 , even though this counts as a backward difference. In particular, the expressions on the upper and lower branches should be the complexconjugate of each other. This requirement is obvious when we calculate the expectation value of φ m in Quantum Mechanics, given the wavefunction Ψ = φm Dφ exp i tm t 0 dtL Ψ 0 . As a bonus of eqn. (2.3), the complex-conjugation enforces all φ m terms in the exponent to vanish, except the temporal derivative terms, since the Lagrangian is real. This conclusion also holds true for quantum field theory. That is to say, as long as t c is chosen consistently, the exponent in (2.2) does not contain φ m terms other than through φ m φ ± m−1 . This linear dependence on φ m is crucial in what follows. 3 While φ ± prove useful in the practical numerical calculation, another basis of fields 1 This initial could be Gaussian, or in some other distribution that is easy to sample. 2 These states are eigenstates of the operatorφ(tj). One may reach the same result by using the unitary time-evolution operator. 3 We have inserted only a single φm field at tm, which is manifested in eqn. (2.3). Actually, even if we insert both upper and lower fields at tm, such as φ ± m , the two fields would be the same, since . Therefore, we can always integrate out one of them, and this leads us back to the single field insertion at tm. more easily reveals the separation of variables, and so we introduce the Keldysh basis where we adopt the coefficients of the transformation from [39,40], but the name for the fields from [41]. It is important to realize that the notation (φ cl , φ q ) does not mean φ q is "quantum" and φ cl is "classical", as they are both required for the full quantum description. The naming convention of φ cl follows from the behaviour of this field at the critical point of the action, and its role in the classical-statistical approximation to be discussed later.
Having introduced (φ cl , φ q ), we can now split the partition function into two parts: the initial density matrix and the dynamic part of the path integral. Since it is natural to pair up φ q j and φ cl j+1 in the (φ cl , φ q ) basis, we attribute all terms that only consist of φ cl 0 , φ q 0 and φ cl 1 to the initial density matrix part. Nevertheless, the dynamic part still contains φ cl 0 and φ cl 1 , but not φ q 0 . For the free initial density matrix, we can further integrate out the dummy variable φ q 0 [35], and this integral is in fact a Wigner-Weyl transform of the initial density matrix [42], which for the quantum field theory in momentum space is of a form, , provided that both f (x) and g(x) in the configuration space are real functions. In practice, the splitting proceeds as Dφ q j Dφ cl j+1 e −I , (2.6) where from the first line to the second 4 we arrange the terms in the exponent into two parts, according to whether or not the terms contain φ q 0 , and −I refers to the part without φ q 0 . Obviously, the conjugate momentum field π cl 0 is defined via which is also the initialφ cl 0 , up to the order of O (dt) 2 . For one thing, we find that the Wigner function of the initial density matrix for the free thermal state reads, up to an overall constant factor, (2.8) 4 We specified the second line by imposing the theory to be free at t0, that is, considering only the quadratic terms in the potential. Otherwise we should include extra non-linear terms of φ q 0 . In particular, we have in mind to turn on the interaction gradually or instantly by varying the coupling constants.
where the particle number n p follows the Bose-Einstein distribution, n p = 1/(e ωpβ − 1), and the frequency ω p refers to the frequency of the free theory, namely ω p = p 2 + m 2 . Concretely, the initialization satisfies, according to which we can draw random numbers for practical simulations as the initialization. In particular, we are going to useφ cl 0 andφ cl 1 to denote those random numbers generated by this way, and as a convention of the paper, the fields with tildes always refer to the classical trajectories satisfying the classical equation of motion, except of courseφ cl 0 andφ cl 1 , which serve as the initialization satisfying (2.9). For the dynamic part of the path integral, the integrals over all field variables not belonging to the initial condition, there are a number of different ways to express the exponent, depending on which basis we use. In order to better see the classical-statistical approximation, we can use the basis introduced in (2.4), the (φ cl , φ q ) basis, for which we find Here we have defined in order to highlight the connection with the equation of motion. This expression makes it manifest that the exponent I is an odd function of φ q , and we will call those non-linear terms of φ q the quantum vertices. The classical-statistical approximation then states that we drop the quantum vertices in (2.10), in which case we may perform the Dφ q part of the path integral, which leaves us with a delta function ∼ δ(φ cl j (x) − φ cl j (x)) for j > 1, while we still need to perform the path integral over the initial condition φ cl 0 and φ cl 1 . In practise, this means that φ cl satisfies the classical equations of motion, where the initial data, φ cl 0 and φ cl 1 , is drawn from the Gaussian distribution (2.9). This approximation has been used in many calculations, including early-Universe simulations (for instance [3,4,43,44]), bubble nucleation [42,45,46] and thermalization [8,9,39]. It is clear from these expressions that the exponent I is purely imaginary, and so the exponential term in (2.6) is purely a phase, leading to a highly oscillatory integral. Dealing with this oscillatory behaviour is the main role of Picard-Lefschetz thimbles [17,47], and in this paper we use the generalized thimbles of [28][29][30] to compute the integral.

Lefschetz thimble and Generalized Thimble Method
The central idea of thimble methods is to start with an integration of some function, e −I(X) , over real variables, X j , and then promote the X j to complex variables, Z j , whilst also promoting the integrand to a holomorphic function e −I(Z) [17-31, 35-37, 47, 48]. The problem then corresponds to integrating along the surface R n ⊂ C n given by Im(Z j ) = 0, but this surface may be deformed 5 to some n−dimensional surface M ⊂ C n without affecting the value of the integral, due to Cauchy's theorem for contour integrals in the complex plane. It is then convenient to perform a coordinate transformation from this curve back to the initial real surface, parametrized by the X j The trick to making the integrals converge is to pick an appropriate surface M, and that is done with the aid of the Picard-Lefschetz flow equations, where the second equation can be derived from the first. As initial conditions for the flow, one takes Z j = X j , and sets J jk = ∂Z j /∂X k to be the n × n identity matrix. In this way, each point X on the real surface is flowed into the complex plane. The surface M is called a generalized thimble [28][29][30], and its shape depends on how long one chooses to flow the equations. In the limit of infinite flow time one recovers the Lefschetz thimbles.
In the case of field theory, our φ i in (2.6) are promoted to complex variables, in analogy to the X j → Z j above. 6 Along the flow, Re[I] always increases, so the weight e −Re [I] decreases exponentially, making the integrand of (2.6) take the form of damped oscillations, leading to a well-behaved integral. In practice, it is convenient to combine this with the determinant of the Jacobian and construct a probability distribution function P = e −Re[I]+ln |det(J)| for use in the Monte Carlo evaluation of the integral [27,35]. Now, by generating samples according to this probability distribution function, we are able to compute the expectation value of any observable, which is done by reweighting the observable by e −iIm[I]+iarg(det(J)) , as this is what remains of the e −I after we have accounted for P .
For a single initialization, i.e. one choice ofφ cl 0 andφ cl 1 , we can measure the expectation value due to its single generalized thimble as where the brackets ... P correspond to an average using P as the probability distribution function. To compute the full quantum average, denoted ... , we need to average over an ensemble of initializationsφ cl 0 ,φ cl 1 drawn from a Gaussian as determined by the initial distribution (2.9) [35]. 5 The curve must not cross any poles, and the integrand is assumed to vanish asymptotically on the deformed curved. 6 To avoid confusion: The time variable indexing the fields φ lives in a complex plane on the Keldysh contour. The time index is part of the j label on Xj. But now the domain of integration, the values taken by the field variables φj, is also deformed into the complex plane.

Critical points
Of central importance to our discussion are the critical points (or critical configurations), the fixed points of the flow, ∂I/∂Z = 0. These are the classical trajectories, and so as the manifold M is deformed by the flow, it does so while the classical trajectories stay pinned at their real values. All configurations away from the critical points give an exponentially suppressed contribution (depending on how much we flow), while all the critical points contribute with equal weight. Hence it is essential to identify and include all critical points, and their neighbourhoods in configuration space. By definition, we only need to solve ∂I ∂φ q j (x) = 0, and ∂I ∂φ cl j+1 (x) = 0, for all x and j ≥ 1, (2.14) with the explicit expression of I in (2.10). The linear dependence on φ m , which is treated as a φ cl field, together with the odd dependence on φ q , allows us to solve the equations exactly. Firstly, since φ m (x) only couples to φ q m−1 (x), we obtain φ q m−1 (x) = 0 as the solution of ∂I/∂φ m (x) = 0. Then, given φ q m−1 (x) = 0 for all x, we can derive φ q m−2 (x) = 0 from ∂I/∂φ cl m−1 (x) = 0. In fact, we can continue this process to discover that all φ q (x) = 0 at the critical points. With this in mind, ∂I/∂φ q j (x) = 0 leads to the classical equation of motion for j ≥ 1, where we use tilde to refer to the solution of φ cl , given the initializationφ cl 0 (x) andφ cl 1 (x). Unsurprisingly, we have shown that the stationary points are the classical trajectoriesφ cl . And by sampling the entire ensemble of initial conditionsφ cl 0 (x),φ cl 1 (x) we can ensure that we include all of them.
For a given initial conditionφ cl 0 (x),φ cl 1 (x), we then have a unique classical trajectory, and we proceed to generate a Markov chain of random trajectories starting from there. The resulting trajectories are then, in general, no longer fixed points of the flow, and are therefore evolved into some complex-valued trajectory through applying the flow equation. The thimble corresponding to the original classical trajectory is the collection of all such nearby trajectories flowed into the complex space.
Still, had all the field variables φ j been part of this Monte-Carlo sampling, one could imagine at some stage randomly generating another classical trajectory. This is then the fixed point corresponding to another neighbourhood of trajectories, belonging to another thimble than the original one. This means that when sampling all the field variables simultaneously, one should take care to properly include all the thimbles, in effect finding all the neighbourhoods of classical trajectories.
In our case, we sample the initial conditions separately, and keep them fixed, while the Monte-Carlo sampling of the other variables takes place. Thus as an initial value problem, the existence and uniqueness of the critical point in the complex space are guaranteed. This means that we will never at random end up on a different classical trajectory, because there is only one classical trajectory corresponding to given values ofφ cl 0 (x),φ cl 1 (x) (or indeed the field variables at any two times, we choose to keep fixed). And so in our two-step sampling, we are sure to find all thimbles and that in any Markov chain, only one thimble is present.
We stress that the conclusion here applies to the quantum field theory with any number of spatial dimensions.
In summary, our approach has a number of advantages for real-time simulations using the Keldysh contour: It admits the two-part splitting into an "external" average over realisations of the initial density matrix and the "internal" Monte-Carlo sampling of the rest of the field variables in the path integral. We can generate the initialization according to the initial density matrix, and with each external initialization there exists one and only one critical point/thimble for the internal part. Thirdly, the critical point (φ q = 0, φ cl =φ cl ) lies in the real plane, as long asφ cl 0 andφ cl 1 are real, which is the case in our approach. Thus the Lefschetz thimble and the generalized thimble coincide at the unique critical point. At the moment, it is easier to implement an algorithm to generate samples on the generalized thimble than on the Lefschetz thimble, so we are going to adopt the Generalized Thimble Method in the following calculation. For more details and numerical applications, see [35].

Tunnel splitting of the ground state
Having set up our formalism and algorithm, we wish to apply it to the inherently quantum mechanical process of tunneling. We therefore introduce a system with a double well potential, which corresponds to m = 1, g 3 = g, λ = 12g 2 in (2.1), and we shall take = 1 in the simulations. In the following, we will consider field theory in zero spatial dimensions, i.e. quantum mechanics. The two minima of the potential are situated at φ = 0 and φ = 1/g respectively, and the height of the barrier between the two minima is 1/(32g 2 ). We set the system up with an initial Gaussian vacuum state of the φ = 0 minimum. When g is small, the separation between the two minima is large and the barrier is high, but the wave function will nevertheless oscillate between the two minima. This is a textbook problem [49,50] known as tunnel splitting of the ground state, which is why we use it as a test of the thimble approach.
To get an idea of how the system behaves we may start off by making the approximation of focussing only on the first two energy eigenstates. These are approximately given by the superposition of the two Gaussian vacuum states located around each minimum, where Ψ L is the Gaussian wavefunction associated to the φ = 0 minimum, and Ψ R is the Gaussian wavefunction associated to the φ = 1/g minimum.
Thus if we take the initial state to be Ψ L , we expect the wave function to evolve as, In particular, if we measure the average φ, we would obtain dφ Ψ * (t)φΨ(t) dφ where, for last relation, we utilize the fact that dφΨ L φΨ L = 0 and dφΨ R φΨ R = 1/g. This shows that the oscillation is described by a single parameter ∆E, which is the energy difference between the first two energy eigenstates. In the literature, there exist many ways to compute the energy difference, here we mention three of them: WKB estimate of ∆E: The textbook of Landau and Lifshitz [49] provides a WKB analysis of the problem, leading to the formula, where a and b are turning points, determined by V − E 0 = 0. In theory, E 0 should be the average energy between the first two states, but this may be no easier to calculate than ∆E. However, in practice we may take the free vacuum energy, 0.5, to be an estimate of E 0 . As g becomes larger, g > 0.25 in our example, then the height of the barrier is lower than 0.5, meaning that the WKB formula (3.5) is no longer valid, so we can only trust this approximation for g 0.25.

Instanton estimate of ∆E:
The tunnelling property can also be captured by an instanton calculation, and the energy difference is estimated via a simple formula [50,51], The formula holds when g is small and we note that one does not need to know the average energy E 0 a priori.
Finite Fock space estimate of ∆E: Finally, for this simple problem, we can find ∆E to high precision by approximating the Fock space with a finite set of vectors [52] |n = 1 √ n! (a † ) n |0 , a † |n = √ n + 1|n + 1 , a|n = √ n|n − 1 , (3.7) and write φ = (a † + a)/ √ 2 into a matrix, The lowest two eigenvalues of H = H 0 − gφ 3 + g 2 φ 4 /2 are then the energy of the ground state and the first excited state, respectively. With a 100 × 100 matrix, we can already achieve a high precision, as will be seen in the next section.  Table 1. Both WKB and Instanton approximations to ∆E work best when g is small. In the cases of g = 0.3 and 0.5, E 0 is already larger than the height of potential barrier, so the WKB approach is not valid.
In table 1, we present ∆E for different g using the three different methods.
Having gained an analytic understanding for the behaviour of the wavefunction, specifically (3.4), we may now solve the full Schrödinger equation numerically, and use our analytic approach to interpret the results.

Solving Schrödinger equation
The time dependent Schrödinger equation reads, and we give the wave function an initial condition associated to the Gaussian ground state of the φ = 0 minimum, It is fairly straightforward to simulate this system with such an initial condition, and for an explicit realization with Python, see [53]. In Fig. 2, we show the expectation value of φ in the time range 0 < mt ≤ 300, with g = 0.3, as calculated with the full Schrödinger equation (solid red curve), and the analytic approximation (3.4) based on the ansatz of (3.3). The approximate solution shows good agreement to the full solution, with the oscillations clearly visible at the expected frequency, along with further structure that is not captured by the two-state approximation. The presence of this extra structure is seen by looking at a snapshot of the wavefunction at some later time, as in (1 − cos(∆Et)), with g = 0.3 and ∆E = 0.285763, matching the matrix approximation found in Tab. 1. Also notice there exist oscillations with higher frequencies, which indicate that other eigenstates contribute too, but they are not as dominant as the first two eigenstates.  Given that the two-state approximation yields good agreement with the full solution we may introduce the idea of the tunnelling time, given by t tun = π ∆E , and we see from Tab. 1 that small g leads to longer tunnelling times, as is to be expected given that in this regime the barrier height is increased, and the vacua are further separated. We should note here that one can in fact rescale φ to remove g, at the cost of redefining → g 2 . Here we would interpret the small g limit as a small limit, and again we expect the tunnelling times to be increased in this limit.

Results with the thimble approach
We adopt the algorithm used in [35] to implement the Markov Chain Monte Carlo evaluation of the integral. In practice, we first generate initialφ cl 0 andφ cl 1 according to the Gaussian initial distribution (2.9). Then, for each initialization, we utilize a Metropolis-Hastings algorithm to generate samples on the generalized thimble. Specifically, according to the proposal distribution, the algorithm is classified as a Random Walk Metropolis-Hastings. A robust choice of the proposal function is exp −(dx) T J † n J n dx/δ 2 , given in [27]. Notice the proposal function above is expressed in the real space, and it corresponds to exp −(dz) † dz/δ 2 in the complex space, given that dz = Jdx. Therefore the covariance matrix is proportional to the identity matrix in the complex manifold. Random Walk Metropolis-Hastings algorithm becomes hard to explore the target space, when the number of integration variables is large, and in practice we find that when the total number of variables is larger than ∼20, each Markov chain will require O(100) CPU-hours or more. Thus, in order to guarantee we achieve ergodicity for each chain, we restrict to 18 field variables, and spend more computational time in increasing the number of updates.
The numerical challenge we encounter here also reflects the fact that the pure quantum tunnelling effect is hard to simulate, as it is accompanied with exponential suppression of the evolution across the barrier, for instance as in [54]. To have a sensible signal of tunnelling, a reasonably-sized coupling constant is required. But for the double-well potential, when the quantum tunnelling effect is big, so is the classical-statistical effect, whereby the field has enough energy to classically roll over the barrier for a significant proportion of the random initial conditionsφ cl 0 andφ cl 1 . Therefore, we want the simulations to be able to distinguish the two effects with a reasonable number of field variables. We shall present results for g = 0.3 and 0.5, both of which correspond to g > 0.25, and so the height of potential barrier is close to, but just below, 0.5, i.e. the Gaussian vacuum energy. As we start the system in the Gaussian ground state, this means that the average energy of the particle is close to the barrier height, meaning that the classical-statistical approximation is just about able to probe the second minimum.
In Fig. 4 we show the results for g = 0.3. The plot on the left gives the solution for the full Schrödinger equation (solid curve), and the results for the classical-statistical approximation are shown as the curve with 1σ statistical error bars. For this value of the coupling the barrier is located at φ = 5/3, and it is clear that the classical-statistical approximation starts to break down as the expectation value moves into the second minimum. The plot confirms our expectation that the classical-statistical approximation can not capture what really happens during the quantum tunnelling. The right hand plot of Fig. 4 shows a comparison of the three methods, the path integral and the Schrödinger methods are seen  to agree within error bars, as they should, while the classical-statistical approximation is starting to deviate at later times.  The simulations with a slightly higher coupling, g = 0.5, correspond to a smaller tunnelling time, and so behave more quantum mechanically. The results of this case are shown in Fig. 5, where we see that the Schrödinger and the thimble approaches again agree within error bars, but now the classical-statistical approximation is showing a more pronounced deviation than in the g = 0.3 case. For g = 0.3, we have implemented 400 initializations, each with 20 millions updates, for g = 0.5 we used 200 initializations, each with 50 millions updates.
The key difference between the classical-statistical approximation and the full thimble approach is that the approximation only uses the critical point for each set of initial conditions, i.e. the solution to the classical equations of motion for that initial data. The full thimble, however, is able to explore more of the field space in order to get a better estimate for the path integral.
In our simulations, the thimble is a 16 (real-)dimensional subspace of C 16 , and so it is difficult to visualize. But we can show a comparison the field values at the critical point versus the average over the whole generalized thimble. In Fig. 6 we present the data for a single Monte-Carlo chain, for one choice of initial data. This shows the statistics of φ cl i for each time slice of the path integral in the g = 0.3 case, where the tunnelling rate is suppressed by the larger barrier height. The solid red lines show the solution to the classical equations of motion for the same initial data, which will form one member of the ensemble of the classical-statistical approximation. We see that at early times the generalized thimble values match the red lines, but deviate at later time slices, leading to the discrepancy between the full quantum behaviour and the classical-statistical approximation.

Conclusion
We have provided a calculation of quantum mechanical tunnelling using the generalized thimble approach [28][29][30], and compared it to the full Schrödinger equation computation, as well as the classical-statistical approximation [39]. The thimble method presented by the authors in [35] explicitly breaks up the path integral into initial conditions (external) and a dynamic (internal) part, and so is ideally suited to understanding the classical-statistical approximation, where the path integral is approximated by summing over solutions to the classical equations of motion. We have demonstrated that the thimble calculation of the real-time path integral does indeed reproduce the expected results of the Schrödinger equation, and we have seen how it improves on the classical-statistical approximation by exploring more of the field trajectory space required in calculating the path integral. While the Schrödinger equation is by far the superior method for calculating these tunnelling processes in quantum mechanics, the Schrödinger functional method is not wellsuited to numerical studies of quantum field theory, which is where we expect the thimble approach to be superior in understanding real-time, non-perturbative quantum dynamics of (at least) scalar field theories.
In the present work, we adopt the Generalized Thimble Method, but because of the uniqueness of the thimble for each initialization, it is also possible to sample directly on the Lefschetz thimble. The latter approach might be optimal when there are more integration variables. In either case, it becomes more and more difficult to explore the manifold of the thimble through a Markov chain, as the dimension of the manifold increases. Also, it seems that the tunnelling process in itself also amplifies the difficulty.
The numerical challenge is substantial. Given n integration variables, the most time consuming part of solving the flow equations is to calculate det(J) which alone requires O(n 3 ) computational time (or less, see [27]). Another large factor is the total number of updates of the Markov chain, whose exact form depends on the algorithm used. To explore the manifold of each thimble, we are now using the Random Walk Metropolis-Hastings algorithm, which on the tunnelling problem behaves less efficiently in higher dimensions. The hope is that we can find some efficient Monte Carlo algorithm, which can maintain O(n) or some other polynomial power. This would mirror the efficiency gained for imaginary time lattice simulations through Hybrid Monte Carlo or Metropolis with local update. An interesting option is to implement a machine learning algorithm to better estimate the thimble [55,56]. Another is to perhaps adapt a complex Langevin algorithm [48,57,58] to explore along a thimble rather than in the entire complexified manifold.