Efficient Adaptive Computation of the Dynamics of Mott Transistors

1 Introduction

We study systems of linear ordinary differential equations of Schrödinger type

ψ^{'} (t) = - i H (t) ψ (t) =: A (t) ψ (t), t \in [t_{0}, t_{e n d}], ψ (t_{0}) = ψ_{0} given,

(1)

with a time-dependent Hermitian matrix $H : R \to C^{d \times d}$ . The exact flow of (1) is denoted by $E (t; t_{0}) ψ_{0}$ 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 [1]. 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 [2].

We consider the numerical solution of the system by exponential-based Magnus-type time integrators in conjunction with an adaptive Lanczos method [3, 4],

u_{n + 1} = S (τ; t_{n}) u_{n}, where

	$S (τ; t_{n}) = S_{J} (τ) \dots S_{1} (τ) = e^{Ω_{J} (τ)} \dots e^{Ω_{1} (τ)},$
	$Ω_{j} (τ) = τ B_{j} (τ), j = 1, \dots, J, B_{j} (τ) = K \sum k = 1 a_{j k} A_{k} (τ), A_{k} (τ) = A (t_{n} + c_{k} τ) .$

The coefficients $a_{j k}$ , $c_{k}$ are determined from the order conditions such that the method attains convergence order $p$ . In this study, we will use the methods referred to as CF2 (exponential midpoint rule), CF4oH and CF6n in [5], 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 $N$ discrete locations, the state space has dimension $4^{N}$ . 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 [3].

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 $V_{g}$ , can be described within the Hubbard model [6], 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 $N$ , 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 $U_{i j σ σ^{'}}$ . In the considered model the Hamiltonian in (1) has the form

H (t) = H_{s t a t} + H_{d y n} (t) .

(2)

The static part for arbitrary hopping $v_{i j}$ reads

H_{s t a t} = \sum i, j, σ v_{i j} {^c}_{j σ}^{†} {^c}_{i σ} + \frac{1}{2} \sum i, j, σ, σ^{'} U_{i j σ σ^{'}} {^n}_{i σ} {^n}_{j σ^{'}}

(3)

where $i, j$ sum over all $N$ sites and the spins $σ, σ^{'} \in {↑, ↓}$ are either up or down. The notation ${^c}_{j σ}^{†} {^c}_{i σ}$ describes the “hopping” of an electron from site $i$ to $j$ with creation and annihilation operators ${^c}_{j σ}^{†}$ and ${^c}_{i σ}$ . The hopping amplitudes $v_{i j}$ with $i, j = 1, \dots, N$

give the probability (rate) of such an electron hopping;

{^n}_{j σ} = {^c}_{j σ}^{†} {^c}_{j σ}

is the occupation number operator, which counts the number of electrons with spin

σ

at site

j

. 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.

U_{i j σ σ^{'}} = U δ_{i j} (1 - δ_{σ σ^{'}})

The dynamic part of the Hamiltonian is in general given by

H_{d y n} (t) = N_{α} \sum α = 1 (R (f_{α} (t)) H_{s y m m}^{α} + i I (f_{α} (t)) H_{a n t i}^{α}) + N_{β} \sum β = 1 g_{β} (t) H_{p o t}^{β} .

(4)

The real matrices $H_{s y m m}^{α}, H_{a n t i}^{α},$ and $H_{p o t}^{α}$ have the following properties: $g_{α} \in R$ , $H_{s y m m}^{α}$ is symmetric, $H_{a n t i}^{α}$

is skew-symmetric,

H_{p o t}^{β}

is diagonal. Here, we restrict ourselves to the Coulomb gauge where

f_{α} (t) = 0 \forall α,

N_{β} = 1.

For the performance analysis we choose a rapidly attenuating potential

g_{1} (t) \equiv g (t) = 1 - \frac{1}{e^{(t - t_{0}) / T} + 1}, t_{0}, T \in R^{+}, H_{p o t}^{1} = \sum i σ V_{i} {^n}_{i σ} =: H_{p o t}, V_{i} \in R .

(5)

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 [3]:

The local time-stepping error can be expressed in terms of an iterated defect $D_{i} (τ; t_{0})$ , where

D_{0} = D = ψ^{'} (t) - A (t) ψ (t), D_{i + 1} = D_{i}^{'} - A D_{i}, i = 0, 1, \dots .

For the exponential midpoint rule, we thereby obtain for the local error $L (τ; t_{0}) ψ_{0} := (S (τ; t_{0}) - E (τ; t_{0})) ψ_{0}$

	$L (τ; t_{0})$	$=$	$\int_{0}^{τ} Π (τ, σ_{1}) D (σ_{1}; t_{0}) d σ_{1} = \int_{0}^{τ} Π (τ, σ_{1}) \int_{0}^{σ_{1}} Π (σ_{1}, σ_{2}) \int_{0}^{σ_{2}} Π (σ_{2}, σ_{3}) d σ_{3} d σ_{2} d σ_{1} \cdot D_{2} (0; t_{0}) + O (τ^{4})$
		$=:$	$I_{3} (τ)      = O (τ^{3}) \cdot (Γ^{''} (0) - A^{''} (t_{0})) + O (τ^{4}), % where \frac{d}{d τ} e^{τ A (τ / 2)} = Γ (τ) e^{τ A (τ / 2)} .$

It follows from the arguments in [3] that more generally, the leading local error term for a higher-order CFM of order $p$ depends on the difference $Γ^{(p)} (0) - A^{(p)} (0)$ .

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 $A \in C^{3}$ , then the local error satisfies

∥ L (τ; t_{0}) ∥_{2} \leq \frac{1}{12} τ^{3} ∥ ∥ [A (t_{0})), A^{'} (t_{0})] - \frac{1}{2} A^{''} (t_{0}) {∥ ∥}_{2} + O (τ^{4}) .

(6)

For the Hubbard model with the potential (5), the convergence analysis above has the following implications.

Differentiation of $g (t)$ shows the following asymptotics:

$g (t)$	$=$	$1 - \frac{1}{e^{(t - t_{0}) / T} + 1} \in (0, 1),$	(7)
$g^{'} (t)$	$=$	$\frac{1}{T} \frac{e^{(t - t_{0}) / T}}{(e^{(t - t_{0}) / T} + 1)^{2}} \in (0, 1 / T), g^{''} (t) = \frac{1}{T^{2}} (\frac{e^{(t - t_{0}) / T}}{(e^{(t - t_{0}) / T} + 1)^{2}} - \frac{2 e^{(t - t_{0}) / T}}{(e^{(t - t_{0}) / T} + 1)^{3}}) \in (0, 1 / T^{2}),$
	$⋮$
$g^{(p)} (t)$	$=$	$O (T^{- p}) .$	(9)

When considering the asymptotics of the error as a function of $T$ , clearly the local error is dominated by the highest derivative of $A (t)$ , where the time-dependence occurs in $g$ . For the exponential midpoint rule, for example, (6) shows that this is $A^{''}$ . Analogously, for a method of order $p$ the local error is dominated by $Γ^{(p)} (0) - A^{(p)} (0)$ and is thus proportional to $τ^{p + 1} T^{- p}$ .

4 Numerical Results

Figure 1: *Left*: Schematic depiction of the Mott transistor with the source–drain voltage applied to the outer sites [described by the potential $V_{s d} (t)$ ] and gate voltage ( $V_{g}$ ). *Right*: Time dependence of the charge at the sites $n_{i} (t)$ and of the source–drain potential $g (t)$ . Computed with DoPri45 (solid lines) and CF4oH (crosses).

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 $V_{s d}$ . The scalar potential added on the mid-sites ( $V_{g}$ 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 $V_{g} = 0 = V_{i}, i \in 2, 3, 4, 5, 6, 7$ (conducting state). For the source–drain potential we used $V_{1} = - V_{8} = V_{s d} / 2 = 1.04 U$ ; $T = 2^{- 5}$ . The potential is switched on at $t_{0} = 10$ and $ψ_{0}$ is chosen as the (unique) ground state of the model. Hubbard interaction $U = 10 | v_{12} |$ . The error displayed is $| | ψ (t = 50) - ψ^{r e f} (t = 50) | |_{2}$ , where the reference solution $ψ^{r e f}$ was computed with adaptive CF4oH with a high accuracy of $tol = 10^{- 11}$ . The tolerance of the Lanczos iteration was $10^{- 12}$ throughout. In the right plot of Figure 1, the index of $n_{i} = n_{i ↑} + n_{i ↓}$ corresponds to the sites labeled in the left part of the figure, and $n_{i σ} (t) = ψ^{†} (t) {^n}_{i σ} ψ (t)$ . 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 $t o l = 10^{- 11}$ . 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.

Figure 2: Comparison of the error as a function of the number of matrix–vector multiplications for adaptive time stepping (left) and equidistant time steps (right).

Acknowledgments

This work was supported by the Austrian Science Fund (FWF) under grant P 30819-N32.

References

Zhong et al. [2015] 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. [2020] J. Li, D. Golez, G. Mazza, A. J. Millis, A. Georges, and M. Eckstein, Phys. Rev. B 101, p. 205140 (2020).
Auzinger et al. [2019] 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 [2020] T. Jawecki, W. Auzinger, and O. Koch, BIT 60, 157–197 (2020).
[5] 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 [1963] J. Hubbard, Proc. Roy. Soc. London A 276, 238–257 (1963).
Mahan [1993] G. Mahan, Many-particle physics, 2nd ed., Physics of solids and liquids (Plenum Press, New York, 1993).
Pavarini et al. [2016] 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.

Efficient Adaptive Computation of the Dynamics of Mott Transistors

Authors

Mesh adaptivity for quasi-static phase-field fractures based on a residual-type a posteriori error estimator

ISALT: Inference-based schemes adaptive to large time-stepping for locally Lipschitz ergodic systems

Efficient Magnus-type integrators for solar energy conversion in Hubbard models

Learning to Reason With Adaptive Computation

A Lepskiĭ-type stopping rule for the covariance estimation of multi-dimensional Lévy processes

A local Gaussian filter and adaptive morphology as tools for completing partially discontinuous curves

Online Selection of CMA-ES Variants

1 Introduction

2 Model of a Mott Transistor

3 Error Analysis

4 Numerical Results

Acknowledgments

References

Efficient Adaptive Computation of the Dynamics of Mott Transistors

Authors

Related Research

Mesh adaptivity for quasi-static phase-field fractures based on a residual-type a posteriori error estimator

ISALT: Inference-based schemes adaptive to large time-stepping for locally Lipschitz ergodic systems

Efficient Magnus-type integrators for solar energy conversion in Hubbard models

Learning to Reason With Adaptive Computation

A Lepskiĭ-type stopping rule for the covariance estimation of multi-dimensional Lévy processes

A local Gaussian filter and adaptive morphology as tools for completing partially discontinuous curves

Online Selection of CMA-ES Variants

This week in AI

1 Introduction

2 Model of a Mott Transistor

3 Error Analysis

4 Numerical Results

Acknowledgments

References