We study systems of linear ordinary differential equations of Schrödinger type
with a time-dependent Hermitian matrix . The exact flow of (1) is denoted by in the following.
We focus on a system that is described by a Hubbard model of electron-hopping between discrete sites. For the numerical experiments we use the parameters corresponding to a simple approximation of a Mott transistor . The time dependence is introduced by coupling of the electronic system to a source–drain potential which is switched on rapidly. The corresponding electromagnetic field is treated classically and enters the equation as a scalar potential and modifies hoppings via Peierls’ substitution .
The coefficients , are determined from the order conditions such that the method attains convergence order . In this study, we will use the methods referred to as CF2 (exponential midpoint rule), CF4oH and CF6n in , the two latter constructed by us with the aim to optimize the error for a given computational effort. For comparison, we also show the embedded Runge–Kutta method by Dormand & Prince DoPri45. The computational challenge results on the one hand from the high dimension of the underlying system. Indeed, for a model with discrete locations, the state space has dimension . Thus, for a model with the claim of physical relevance, the problem quickly reaches the limitations of modern supercomputers. On the other hand, the modeling of the switching process in a Mott transistor, features a quickly attenuating electric field in a small time interval. Thus the smoothness of the solution varies strongly in the course of the time propagation, which suggests to use adaptive choice of the time step to increase the efficiency and reliability of a numerical computation. To this end, we use classical step-size choice based on asymptotically correct defect-based estimators to equidistribute the local error, which were constructed and analyzed in .
2 Model of a Mott Transistor
The Mott transistor, which is based on the Mott-insulator-to-metal transition driven by the change of the gate voltage , can be described within the Hubbard model , the paradigm model for the description of strongly interacting electrons [7, 8]. We resort here to the second-quantization formalism. It describes the electron occupation on a given number of sites , corresponding to Wannier discretization. Only a single orbital per site is considered with nearest-neighbor hopping between the sites. Due to the Pauli exclusion principle, there are only four states per site allowed (no electrons, one electron with spin-up or -down, two electrons with opposite spins). The electrons interact via Coulomb interaction . In the considered model the Hamiltonian in (1) has the form
The static part for arbitrary hopping reads
where sum over all sites and the spins are either up or down. The notation describes the “hopping” of an electron from site to with creation and annihilation operators and . The hopping amplitudes with
give the probability (rate) of such an electron hopping;is the occupation number operator, which counts the number of electrons with spin at site . For details on the notation in (3) refer to [6, 7, 8]. In the following we take only the local part of the Coulomb interaction, i.e. .
The dynamic part of the Hamiltonian is in general given by
The real matrices and have the following properties: , is symmetric,
is skew-symmetric,is diagonal. Here, we restrict ourselves to the Coulomb gauge where For the performance analysis we choose a rapidly attenuating potential
3 Error Analysis
In this section, we investigate the error structure of commutator-free Magnus-type time integrators for the model of a Mott transistor. To this end, we recapitulate general convergence results given in :
The local time-stepping error can be expressed in terms of an iterated defect , where
For the exponential midpoint rule, we thereby obtain for the local error
It follows from the arguments in  that more generally, the leading local error term for a higher-order CFM of order depends on the difference .
For the exponential midpoint rule of order two, this means in particular that the following
error estimate holds:
[3, Proposition 4.1] Consider the solution of (1) by the exponential midpoint rule. If , then the local error satisfies
For the Hubbard model with the potential (5), the convergence analysis above has the following implications.
Differentiation of shows the following asymptotics:
When considering the asymptotics of the error as a function of , clearly the local error is dominated by the highest derivative of , where the time-dependence occurs in . For the exponential midpoint rule, for example, (6) shows that this is . Analogously, for a method of order the local error is dominated by and is thus proportional to .
4 Numerical Results
As a numerical illustration, we show the adaptive solution of our model of a Mott transistor by Magnus-type integrators. The geometry of the model is illustrated in Figure 1 (left). The dynamics of the transistor is approximated by the dynamics of an 8-site chain with 8 electrons. The sites are then split into three sections which correspond to an electron reservoir (left and right) which is coupled to a source–drain potential . The scalar potential added on the mid-sites ( in the left part of Figure 1) models the gate of the transistor. It controls whether the transistor is in a conducting or an insulating state. For the numerical analysis shown here we chose (conducting state). For the source–drain potential we used ; . The potential is switched on at and is chosen as the (unique) ground state of the model. Hubbard interaction . The error displayed is , where the reference solution was computed with adaptive CF4oH with a high accuracy of . The tolerance of the Lanczos iteration was throughout. In the right plot of Figure 1, the index of corresponds to the sites labeled in the left part of the figure, and . It shows the occupation numbers (proportional to the local charge) computed with a fine error tolerance by the Dormand & Prince Runge–Kutta method (solid lines), and the adaptive Magnus-type method CF4oH (crosses) with the same error tolerance of . We observe that the solution computed with adaptive time stepping chooses appropriately small steps where the solution varies strongly, and very large steps where the solution is smooth. Nonetheless, the solutions correspond very well at the grid points. We conclude that adaptive commutator-free Magnus-type methods serve very well to accurately approximate the solution both in regions where dense grids are required and where the dynamics allows for large time steps. The adaptive choice of the steps thus implies a gain in efficiency in contrast to uniform time steps. Indeed, in Figure 2
, we show the error as a function of the number of matrix–vector multiplications, which constitute most of the computational effort required for the integrators. The left plot gives the results for adaptive time stepping, on the right-hand side, the results for equidistant time steps are given for comparison. We observe that the error is smaller for a comparable computational effort when the adaptive strategy is used. This results from the unsmooth solution dynamics associated with the rapid changes in the electric field. The specially constructed CFM integrators are most efficient, adaptive Runge–Kutta methods show consistent convergence behavior, but prohibitively large computational effort.
This work was supported by the Austrian Science Fund (FWF) under grant P 30819-N32.
- Zhong et al.  Z. Zhong, M. Wallerberger, J. M. Tomczak, C. Taranto, N. Parragh, A. Toschi, G. Sangiovanni, and K. Held, Phys. Rev. Lett. 114, p. 246401Jun (2015).
- Li et al.  J. Li, D. Golez, G. Mazza, A. J. Millis, A. Georges, and M. Eckstein, Phys. Rev. B 101, p. 205140 (2020).
- Auzinger et al.  W. Auzinger, H. Hofstätter, O. Koch, M. Quell, and M. Thalhammer, M2AN – Math. Model. Numer. Anal. 53, 197–218 (2019).
- Jawecki, Auzinger, and Koch  T. Jawecki, W. Auzinger, and O. Koch, BIT 60, 157–197 (2020).
-  W. Auzinger, J. Dubois, K. Held, H. Hofstätter, T. Jawecki, A. Kauch, O. Koch, K. Kropielnicka, P. Singh, and C. Watzenböck, “Efficient Magnus-type integrators for Hubbard models of solar cells,” Submitted.
- Hubbard  J. Hubbard, Proc. Roy. Soc. London A 276, 238–257 (1963).
- Mahan  G. Mahan, Many-particle physics, 2nd ed., Physics of solids and liquids (Plenum Press, New York, 1993).
- Pavarini et al.  E. Pavarini, E. Koch, J. van den Brink, and G. Sawatzky, Quantum Materials: Experiments and Theory, Modeling and Simulation, Vol. 6 (Forschungszentrum Jülich, Jülich, 2016) p. 420 p.