Thermodynamic quantities of small systems fluctuate and have to be characterised by their appropriate probability distributions. Within the two-point energy measurement prescription, the distribution of work in a quantum system can be derived from the transition probability from initial to final energy. We consider a simple yet representative model system starting in thermodynamic equilibrium and driven by an external force, and compare two different numerical techniques to determine this transition probability with respect to accuracy and numerical effort. In addition, we perform a semi-classical analysis of the process using the WKB approximation. The results agree well with the numerically exact values if Airy-tails modelling the tunnelling into classically forbidden regions of phase space are properly taken into account.
The ongoing miniaturisation of technical devices like engines and frigistors as well as the rapid development of new experimental techniques in biophysics aiming at the analysis of single molecules, molecular motors, and ion pumps require new concepts for the theoretical modelling of energy conversion at the nanoscale. Within the last 20 years thermodynamics has been extended to small systems with energy turnover comparable to their energy fluctuations. The key concept is to describe thermodynamics quantities like work, heat, and entropy by their corresponding probability distributions. The new field of stochastic thermodynamics (for introductory reviews, see , , ) has meanwhile proven to provide the appropriate framework for analysing the efficiency of nanomachines , , , estimating free-energy differences in single-molecule experiments , and highlighting the role of information as thermodynamic resource .
Focussing on small systems invariably brings quantum effects into play. While stochastic thermodynamics has reached a fairly mature state for classical systems, the thermodynamic description of fluctuating quantum systems is much less complete. One reason is that central thermodynamic quantities like work and heat are no state variables, and their proper quantum definition remains controversial , , , . Work performed on a classical system, e.g. depends on the complete trajectory the system undergoes during the process. However, no entire quantum analogue of the system trajectory is available. One possible workaround restricts the analysis to isolated systems for which the work performed or delivered must be equal to the energy difference of the system as required by the first law of thermodynamics . Then, measuring the systems energy twice, first before the driving and second after its end, should yield a first estimate for the work. Although being simple and operative, this prescription has its own deficiencies because it involves two projective measurements that destroy existing quantum correlations. As a result, the system dynamics is merely replaced by a classical stochastic process with just the transition rates derived from quantum mechanics.
After the first energy measurement, the system is in an eigenstate of the initial Hamiltonian. Being isolated from its surroundings, it then follows a unitary evolution specified by the process under consideration. The final state of this evolution determines the probability for the results of the second energy measurement. In the simplest setting, the system starts in equilibrium at some temperature T, so the probability for the first energy value is given by the Gibbs measure. The central quantity of interest determining the histogram of work values is then the transition probability from the initial state to the one corresponding to the value of the energy at the second measurement.
In the present paper, we study this transition probability for a simple but representative and not analytically solvable model system introduced in . We compare two different numerical methods for its determination, characterise their accuracy as a function of the various parameters, and recommend suitable choices for these. Finally, we check our results against an approximate analytic treatment within a semi-classical framework.
The paper is organised as follows. Section 2 contains the basic equations and fixes the notation. Section 3 describes the two numerical methods we use for our comparison. In section 4 we compare the results of these methods with each other and discuss in detail the influence of the various parameters on accuracy and computing time. Section 5 is devoted to the semi-classical analysis. Finally, section 6 contains our conclusions.
We consider the one-dimensional motion of a quantum particle of mass M in a potential
in the time interval with the parameter changing from at t = 0 to at final time tf . The dynamics is described by the time-dependent Schrödinger equation for the wave function of the particle:
By the first energy measurement at t = 0, the system is prepared in an energy eigenstate of the initial Hamiltonian satisfying
with denoting the corresponding energy eigenvalue. For t > 0, the state then follows the unitary time evolution described by
where 𝒯 denotes time ordering. The driving of the system is specified by the explicit time dependence of the Hamiltonian . At the final time , a second energy measurement is performed that projects the state on an eigenstate of the final Hamiltonian with corresponding eigenvalue .
Our central quantity of interest is the transition probability
for the state to be projected onto at when starting in at t = 0. We will only consider protocols linear in time, i.e.
as units for length, time, and energy, respectively, and measuring λ in multiples of λ0, we arrive at the dimensionless form of the Schrödinger equation (2)
with starting at λ0 = 1.
For the classical oscillator, the dependence of the oscillation period T on the amplitude A is given by
Replacing A by the energy E of the particle, this corresponds in dimensionless units to
On the other hand, we get for the energy of the m-th quantum level in the potential from Bohr–Sommerfeld quantisation the approximate expression 
for the classical oscillation period corresponding to the m-th quantum level of energy. We will mainly be interested in values of λ between 1 and 5 and energy levels up to m = 200. The relevant dimensionless oscillation periods then lie in the range 0.2…1. In order to find appreciable probabilities for off-diagonal transitions, , the driving process parametrised by should not be slow on these time scales. Accordingly, for our numerical investigations, we choose and .
We use two different procedures to solve the Schrödinger equation (8) numerically. The first one is the well-known Crank–Nicolson method  that provides a direct prescription to compute the wave function at times t > 0 given its initial form at t = 0. The second method employs a special ansatz for the wave function involving time-dependent coefficients . It transforms the Schrödinger equation (8) into a system of ordinary differential equations for these coefficients that is then solved numerically. In the following, we shortly characterise both methods.
The central point of the Crank–Nicolson method is the discretisation of the time evolution operator (4) in Cayley form
where denotes the temporal step size. The main virtue of this replacement is that the first term on the r.h.s. of (13) is unitary and, therefore, norm preserving. The discretised Schrödinger equation acquires the form
as given by (8).
It is convenient to rewrite (14) in terms of a new function
and using the customary notations , , and where we find
with the abbreviations
The second-order differential equation (17) has to be solved numerically at every time step. To this end, we also discretise the space coordinate with the spatial step size and and use the notation , etc.
From its structure, (17) is amenable to the Numerov method that allows to improve on the accuracy of the solution with very little additional effort. The method is particularly efficient if a recursive procedure is used to perform the necessary matrix inversion (see , ). The main idea is the replacement
as implied by (17) in the approximation
for the second derivative with respect to x. In terms of the auxiliary quantities
one then gets a three-term recursion relation of the form
that approximates the original Schrödinger equation (8) to order and .
The procedure is then as follows. With given from the k-th time step, is determined using (19). The potential fixes according to (18) and, subsequently, by using (21). Therefore, the recursion (23) allows to determine all from two initial values, and . From , we find using (22), which via (16) gives .., the wave function at the next time step:
The final dodge of the algorithm concerns the point that (17) has to be solved as a boundary value problem with (and therefore also ) given for j = 0 and rather than as an initial value problem fixing for j = 0 and j = 1. As detailed in , this problem can be dealt with by introducing yet another set of variables and according to
which transforms the recursion (23) into
The initial conditions and finite ensure and allow to successively calculate all using (26) and (27). From these all, can be determined using (25) by starting with = 0 and decreasing j in each step down to j = 1.
The second method to solve the Schrödinger equation (8) builds on an ansatz for the wave function . We introduce the eigenstates and corresponding eigenvalues of the Hamiltonian for a fixed value λt of the parameter ,
Next, we expand the state at time t on the basis of these instantaneous eigenstates according to
The maximal value N of n in this expansion is a crucial parameter for the accuracy of this method (cf. Section 4). Plugging (29) into the Schrödinger equation (8), it transforms into a set of linear ordinary differential equations for the time evolution of the expansion coefficients :
Here, the second line follows from (28) and the Hellmann–Feynman-theorem.
To solve this system of coupled ordinary differential equations, we need the matrix elements and the instantaneous eigenvalues for each value of λ. Since the stationary Schrödinger equation is not analytically solvable, this can be a time-consuming numerical task. However, if we rescale the space coordinate x in the initial, i.e. , Schrödinger equation (28)
we are left with
Adapting the normalisation of the eigenfunctions to the rescaled coordinate we therefore find the simple mappings
These in turn allow to rewrite the system (31) in the form
Hence, we have to solve the stationary Schrödinger equation only once, namely, for , which can be done, e.g. with the matrix Numerov method . Plugging the results for and into (37) a completely determined system of coupled linear differential equations with time-dependent coefficients for the expansion parameters cn have to be solved numerically. For a better comparison between the two methods, we use for this task a second-order Runge–Kutta algorithm, which has the same accuracy as in the Crank–Nicolson method.
In the present section, we compare the two numerical methods described above with respect to their accuracy as a function of the various parameters. To this end, we study four representative examples for transition probabilities corresponding to the combinations (n, m) = (46, 50), (50, 50), (58, 50) and (50, 46). Their respective numerically exact values obtained with parameter combinations detailed below are compiled in Table 1. We then change the parameters such as, e.g. to reduce the necessary computation time and monitor the resulting deviations from these reference values.
The numerically exact values of the second column were determined by the Crank–Nicolson method with and . The expansion method yields exactly the same values for the digits shown. The semi-classical results displayed in the third column were obtained using classical trajectories. The last column shows the Crank–Nicolson values for a slower process with .
Let us first look at the influence of the temporal and spatial resolution on the results. Figure 1 shows the dependence of the relative accuracy on the spatial step size . Irrespective of the specific values of m and n, both methods yield accurate results for already. This is due to the high convergence rate of the Numerov method discussed in Section 3.
The dependence on the time step is qualitatively similar, but now, step sizes of the order of or smaller are necessary to obtain reliable results. This reflects the fact that the chosen time integration methods are only accurate to order . Nevertheless, both methods are accurate and reliable for sufficiently small increments in x and t.
The plateau values to which converge in Figures 1 and 2 for small step sizes depend on the value of the complementary parameter and , respectively, and for the expansion method also on N. As an example for the interplay of these parameters, we show in Figure 3 examples of heat maps of the accuracy in the - plane. The qualitative behaviour is as expected: both and have to be sufficiently small to obtain high accuracy. But the figure also offers quantitative information. It is again clearly seen that has to be significantly smaller than for getting a desired accuracy. Moreover, comparison of the figures with each other shows that for higher values of m and n, a finer resolution in x and t is necessary to get the same accuracy. This is in line with intuition since higher values of m and n imply shorter wavelengths and higher frequencies of the participating states.
When applying the expansion method, the maximal value N of the expansion parameter is an additional important parameter affecting the accuracy of the results. Although it is obvious that N has to exceed both m and n, it is not so clear to which extent, since during the driving of the system between t = 0 and , a proper approximation of its state may well require eigenstates in the superposition (29) with values of l significantly larger than those of n and m. Figure 4 gives an impression on how the relative accuracies of the considered transition probabilities depend on N. It shows that in the present context, it is sufficient to include roughly 20 more modes than given by the maximum of n and m. Note the large error in if N remains below this threshold.
Guided by the results displayed above, we have chosen the parameters for the numerically exact reference values as given in the caption of Table 1. In the last column of this table, we have included the values of the transition probabilities for instead of , i.e. for five times slower driving. These values show the approach of the transition probabilities to as required by the adiabatic theorem of quantum mechanics. These values, therefore, may serve as an additional independent check of our numerical code.
In addition to the numerical methods discussed above, we give in the following a short report on a semi-classical treatment of the problem. Quite generally, semi-classical methods yield approximations to quantum mechanical results by building on classical trajectories . For the problem at hand, the classical equations of motion have again to be solved numerically.
To obtain a semi-classical approximation for the transition probabilities (5), we need semi-classical expressions for both the eigenfunctions of the final Hamiltonian and the state to which the initial state evolved at the end of the driving. Following , we implement WKB wave functions of the form
where and denote the classical phase space densities corresponding to the states and , and are their classical actions, and and denote the Maslov indices labelling the different possible branches b of classical transitions . For notational simplicity, we suppressed the dependence on n and m in the r.h.s. (39) and (40), respectively.
The scalar product in (5) then acquires the form
In a semi-classical analysis, one is interested in the leading behaviour for ; the integral in (41) may therefore be performed by the stationary phase approximation. Because of
where stands for the classical momentum as function of position, the condition
for the stationary point xs translates to leading order into
The stationary points xs are therefore nothing but the crossing points between the final classical orbit in phase space and the curve . Here Hcl denotes the classical Hamilton function.
Expanding the densities and as well as to second order in , we may replace the sum over the branches b and b′ by a sum over the crossing points and find to leading order in ℏ
Here, b and b′ denote the branches corresponding to the crossing point xs.
Finally, the transition probabilities are given by
The amplitudes as and phases θs are fully determined by the properties of the classical trajectories at the intersection points xs. Integrating the classical equations of motion numerically for a certain number of different initial conditions corresponding to the same initial energy, , we determine the values of ϕs and κs and produce reliable estimates for the densities and . Table 1 shows results for the transition probabilities obtained in this way using 105 initial conditions and compares them to the findings of Section 3.
For and also for , the results from the numerical solution of the Schrödinger equation and those coming from the semi-classical analysis agree rather well. For , there is a larger discrepancy. It is due to the fact that, classically, only a finite interval of energy transfers can be accomplished and the difference lies just outside this energy window. Building solely on classical trajectories, the semi-classical method as described above must then fail. This is a standard problem of the WKB method that can be remedied by a higher-order expansion of resulting in Airy-tails of the semi-classical wave function extending into the classically forbidden region. Proceeding as described in detail in Section V of , we get in this way the refined semi-classical value that agrees much better with the corresponding result of Section 3.
In some circumstances, the semi-classical limit is found to be tantamount to the adiabatic limit . This seems to make a WKB treatment of non-diagonal transition probabilities questionable. However, ℏ cannot be compared with a time scale as such, but only with a product of a time scale and a characteristic energy. By considering rather large energy levels m ≃ 50, we chose this energy scale large enough to keep the small ℏ limit compatible with a comparatively short driving time tf.
In the present paper, we have studied two numerical methods to determine the transition amplitude between energy eigenstates of a simple quantum system driven by a time-dependent external protocol. The corresponding transition probability is the pivotal quantity to compile the work statistics of the system, which in turn determines many of its thermodynamic properties. Our motivation was two-fold: on the one hand, we systematically investigated the accuracy of both methods as function of their parameters, in particular of the temporal and spatial resolution used in the discretisation. On the other hand, we compared the two methods with each other with respect to the numerical effort involved.
Technically speaking, the problem consists of solving the time-dependent Schrödinger equation for the one-dimensional model system. No analytical solution is available. Our first procedure implements the standard Crank–Nicolson method for partial differential equations. The second one builds on an expansion of the quantum state into a superposition of instantaneous eigenstates of the Hamiltonian.
Implementation of the Crank–Nicolson method gives rise to a rather universal solver for the initial value problem of the Schrödinger equation with Dirichlet boundary conditions. It is easily implemented, reliable, and yields accurate results. However, it is a purely mathematical prescription to solve a partial differential equation. There is no way to improve its performance by using additional physical insight into the dynamics of the system.
Expanding the searched-for wave function in the instantaneous eigenfunctions of the system, on the other hand, allows to calculate part of its time dependence analytically. Moreover, studying the dependence of the results on the number of states included valuable information on which states contribute substantially to a specific transition and which do not. On the down side, this method requires to solve the stationary Schrödinger equation of the system for each value of the protocol parameter separately. As far as numerical effort is concerned, the method is hence only competitive with the Crank–Nicolson method if a way is found to circumvent the repeated numerical solution of the eigenvalue problem of the Hamiltonian. For our system, this was possible by a scaling transformation as explained in Section 3.
In both methods, we found the Numerov prescription for the spatial discretisation extremely valuable. With only little additional effort, the accuracy is considerably improved. As a result, rather large spatial step sizes could be used, which reduced the computation time substantially.
Given appropriate values of its parameters, both methods gave consistent results for the transition probabilities considered. In order to additionally verify their correctness, we compared them with the outcome of a semi-classical analysis using the WKB method. Again, very good agreement was obtained if Airy-tails extending into the classically forbidden region were taken into account.
This work has been funded by the Deutsche Forschungsgemeinschaft (DFG; German Research Foundation) – 397082825. We would like to thank Vincent Preut, Tim Utermöhlen, and the members of the DFG Research Unit FOR2692 for fruitful discussions.
 C. Jarzynski, Annu. Rev. Condens. Matter Phys. 2, 329 (2011). Search in Google Scholar
 U. Seifert, Rep. Prog. Phys. 75, 126001 (2012). Search in Google Scholar
 K. Sekimoto, Stochastic Energetics, Springer, Berlin 2010. Search in Google Scholar
 T. Schmiedl and U. Seifert, EPL 81, 20003 (2008). Search in Google Scholar
 M. Esposito, R. Kawai, K. Lindenberg, and C. van den Broeck, Phys. Rev. E 81, 041106 (2010). Search in Google Scholar
 P. Chvosta, M. Einax, V. Holubec, A. Ryabov, and P. Maass, J. Stat. Mech. Theor. Exp. 2010, P03002 (2010). http://dx.doi.org/10.1088/1742-5468/2010/03/P03002. Search in Google Scholar
 C. Chipot and A. Pohorille, Free Energy Calculations: Theory and Applications in Chemistry and Biology, Springer, Berlin 2007. Search in Google Scholar
 J. M. P. Parrondo, J. M. Horowitz, and T. Sagawa, Nat. Phys. 11, 131 (2015). Search in Google Scholar
 S. Yukawa, J. Phys. Soc. Jpn. 69, 2367 (2000). Search in Google Scholar
 A. Engel and R. Nolte, EPL 79, 10003 (2007). Search in Google Scholar
 J. Gemmer, M. Michel, and G. Mahler, Quantum Thermodynamics, Springer, Berlin 2009. Search in Google Scholar
 M. Campisi, P. Hänggi, and P. Talkner, Rev. Mod. Phys. 83, 771 (2011). Search in Google Scholar
 C. Jarzynski, H. T. Quan, and S. Rahav, Phys. Rev. X 5, 031038 (2015). Search in Google Scholar
 L. D. Landau and E. M. Lifshitz, Course of Theoretical Physics III: Quantum Mechanics, Non-Relativistic Theory, Butterworth-Heinemann, Oxford 2005. Search in Google Scholar
 W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes, Cambridge U. P., Cambridge 2007. Search in Google Scholar
 C. A. Moyer, Am. J. Phys. 72, 351 (2004). Search in Google Scholar
 A. Goldberg, H. M. Schey, and J. L. Schwartz, Am. J. Phys. 35, 177 (1967). Search in Google Scholar
 M. Pillai, J. Goglio, and T. G. Walker, Am. J. Phys. 80, 1017 (2012). Search in Google Scholar
 R. G. Littlejohn, J. Stat. Phys. 68, 7 (1992). Search in Google Scholar
© 2020 Walter de Gruyter GmbH, Berlin/Boston