Thermodynamic quantities of small systems fluctuate and have to be characterised by their appropriate probability distributions. Within the twopoint 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 semiclassical analysis of the process using the WKB approximation. The results agree well with the numerically exact values if Airytails 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 [1], [2], [3]) has meanwhile proven to provide the appropriate framework for analysing the efficiency of nanomachines [4], [5], [6], estimating freeenergy differences in singlemolecule experiments [7], and highlighting the role of information as thermodynamic resource [8].
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 [9], [10], [11], [12]. 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 [12]. 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 [13]. 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 semiclassical 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 semiclassical analysis. Finally, section 6 contains our conclusions.
We consider the onedimensional motion of a quantum particle of mass M in a potential
in the time interval
By the first energy measurement at t = 0, the system is prepared in an energy eigenstate
with
where 𝒯 denotes time ordering. The driving of the system is specified by the explicit time dependence of the Hamiltonian
Our central quantity of interest is the transition probability
for the state to be projected onto
Introducing
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
For the classical
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 mth quantum level in the potential
Combining (10) and (11), we find
for the classical oscillation period corresponding to the mth 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
We use two different procedures to solve the Schrödinger equation (8) numerically. The first one is the wellknown Crank–Nicolson method [15] that provides a direct prescription to compute the wave function
The central point of the Crank–Nicolson method is the discretisation of the time evolution operator (4) in Cayley form
where
with
as given by (8).
It is convenient to rewrite (14) in terms of a new function
and using the customary notations
with the abbreviations
and
The secondorder differential equation (17) has to be solved numerically at every time step. To this end, we also discretise the space coordinate
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 [16], [17]). 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
and
one then gets a threeterm recursion relation of the form
that approximates the original Schrödinger equation (8) to order
The procedure is then as follows. With
The final dodge of the algorithm concerns the point that (17) has to be solved as a boundary value problem with
which transforms the recursion (23) into
The initial conditions
The second method to solve the Schrödinger equation (8) builds on an ansatz for the wave function [13]. We introduce the eigenstates
Next, we expand the state
where
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–Feynmantheorem.
To solve this system of coupled ordinary differential equations, we need the matrix elements
according to
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
with
Hence, we have to solve the stationary Schrödinger equation only once, namely, for
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
n  m 




46  50  0.2695  0.2626  0.0023 
50  50  0.1090  0.1089  0.8222 
58  50  0.0449  0.0194 

50  46  0.2010  0.2009  0.0026 
The numerically exact values of the second column were determined by the Crank–Nicolson method 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
The dependence on the time step
The plateau values to which
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
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
In addition to the numerical methods discussed above, we give in the following a short report on a semiclassical treatment of the problem. Quite generally, semiclassical methods yield approximations to quantum mechanical results by building on classical trajectories [19]. For the problem at hand, the classical equations of motion have again to be solved numerically.
To obtain a semiclassical approximation for the transition probabilities (5), we need semiclassical expressions for both the eigenfunctions
where
The scalar product in (5) then acquires the form
with
In a semiclassical analysis, one is interested in the leading behaviour for
where
for the stationary point x_{s} translates to leading order into
The stationary points x_{s} are therefore nothing but the crossing points between the final classical orbit in phase space and the curve
Expanding the densities
with
Here, b and b′ denote the branches corresponding to the crossing point x_{s}.
Finally, the transition probabilities are given by
with
The amplitudes a_{s} and phases θ_{s} are fully determined by the properties of the classical trajectories at the intersection points x_{s}. Integrating the classical equations of motion numerically for a certain number of different initial conditions corresponding to the same initial energy,
For
In some circumstances, the semiclassical limit
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 timedependent 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 twofold: 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 timedependent Schrödinger equation for the onedimensional 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 searchedfor 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
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
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 semiclassical analysis using the WKB method. Again, very good agreement was obtained if Airytails 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.
[1] C. Jarzynski, Annu. Rev. Condens. Matter Phys. 2, 329 (2011). Search in Google Scholar
[2] U. Seifert, Rep. Prog. Phys. 75, 126001 (2012). Search in Google Scholar
[3] K. Sekimoto, Stochastic Energetics, Springer, Berlin 2010. Search in Google Scholar
[4] T. Schmiedl and U. Seifert, EPL 81, 20003 (2008). Search in Google Scholar
[5] M. Esposito, R. Kawai, K. Lindenberg, and C. van den Broeck, Phys. Rev. E 81, 041106 (2010). Search in Google Scholar
[6] 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/17425468/2010/03/P03002. Search in Google Scholar
[7] C. Chipot and A. Pohorille, Free Energy Calculations: Theory and Applications in Chemistry and Biology, Springer, Berlin 2007. Search in Google Scholar
[8] J. M. P. Parrondo, J. M. Horowitz, and T. Sagawa, Nat. Phys. 11, 131 (2015). Search in Google Scholar
[9] S. Yukawa, J. Phys. Soc. Jpn. 69, 2367 (2000). Search in Google Scholar
[10] A. Engel and R. Nolte, EPL 79, 10003 (2007). Search in Google Scholar
[11] J. Gemmer, M. Michel, and G. Mahler, Quantum Thermodynamics, Springer, Berlin 2009. Search in Google Scholar
[12] M. Campisi, P. Hänggi, and P. Talkner, Rev. Mod. Phys. 83, 771 (2011). Search in Google Scholar
[13] C. Jarzynski, H. T. Quan, and S. Rahav, Phys. Rev. X 5, 031038 (2015). Search in Google Scholar
[14] L. D. Landau and E. M. Lifshitz, Course of Theoretical Physics III: Quantum Mechanics, NonRelativistic Theory, ButterworthHeinemann, Oxford 2005. Search in Google Scholar
[15] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes, Cambridge U. P., Cambridge 2007. Search in Google Scholar
[16] C. A. Moyer, Am. J. Phys. 72, 351 (2004). Search in Google Scholar
[17] A. Goldberg, H. M. Schey, and J. L. Schwartz, Am. J. Phys. 35, 177 (1967). Search in Google Scholar
[18] M. Pillai, J. Goglio, and T. G. Walker, Am. J. Phys. 80, 1017 (2012). Search in Google Scholar
[19] R. G. Littlejohn, J. Stat. Phys. 68, 7 (1992). Search in Google Scholar
© 2020 Walter de Gruyter GmbH, Berlin/Boston