# Computational parameters

Discussion of the computational parameters in a TEMPO or PT-TEMPO computation and establishing convergence of results

launch binder (runs in browser),

read through the text below and code along.

## Contents

The following packages will be required

```
import sys
sys.path.insert(0,'..')
import oqupy
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 14.0, 'lines.linewidth':2.50, 'figure.figsize':(8,6)})
```

The OQuPy version should be `>=0.5.0`

```
oqupy.__version__
```

```
'0.5.0'
```

## Introduction - numerical exactness and computational parameters

The TEMPO and PT-TEMPO methods are numerically exact meaning no approximations are required in their derivation. Instead error only arises in their numerical implementation, and is controlled by a set of computational parameters. The error can, in principle (at least up to machine precision), be made as small as desired by tuning those numerical parameters. In this tutorial we discuss how this is done to derive accurate results with manageable computational costs.

As introduced in the Quickstart tutorial a TEMPO or PT-TEMPO calculation has three main computational parameters:

A memory cut-off

`tcut`

, which must be long enough to capture non-Markovian effects of the environmentA timestep length

`dt`

, which must be short enough to avoid Trotter error and provide a sufficient resolution of the system dynamicsA precision

`epsrel`

, which must be small enough such that the numerical compression (singular value truncation) does not incur physical error

In order to verify the accuracy of a calculation, convergence should be
established under all three parameters, under increases of `tcut`

and
decreases `dt`

and `epsrel`

. The challenge is that these parameters
cannot necessarily be considered in isolation, and a balance must be
struck between accuracy and computational cost. The strategy we take is
to firstly determine a suitable `tcut`

(set physically by properties
of the environment) with rough values of `dt`

and `epsrel`

, then
determine convergence under `dt->0,epsrel->0`

.

We illustrate convergence using the TEMPO method, but the ideas
discussed will also generally apply to a PT-TEMPO computation where one
first calculates a process tensor - fixing `tcut`

, `dt`

, `epsrel`

- before calculating the system dynamics (see PT-TEMPO).
Note some of the calculations in this tutorial may not be suitable to
run in a Binder instance. If you want to run them on your own device,
you can either copy the code as you go along or download the .ipynb
file
to run in a local jupyter notebook session. Example results for all
calculations are embedded in the notebook already, so this is not
strictly required.

## Choosing tcut

### Example - memory effects in a spin boson model

We firstly define a spin-boson model similar to that in the Quickstart tutorial, but with a finite temperature environment and a small additional incoherent dissipation of the spin.

```
sigma_x = oqupy.operators.sigma('x')
sigma_y = oqupy.operators.sigma('y')
sigma_z = oqupy.operators.sigma('z')
sigma_m = oqupy.operators.sigma('-')
omega_cutoff = 2.5
alpha = 0.8
T = 0.2
correlations = oqupy.PowerLawSD(alpha=alpha,
zeta=1,
cutoff=omega_cutoff,
cutoff_type='exponential',
temperature=T)
bath = oqupy.Bath(0.5 * sigma_z, correlations)
Omega = 2.0
Gamma = 0.02
system = oqupy.System(0.5 * Omega * sigma_x,
gammas=[Gamma],
lindblad_operators=[sigma_m], # incoherent dissipation
)
t_start = 0.0
t_end = 5.0
```

To determine a suitable set of computational parameters for
`t_start<=t<=t_end`

, a good place to start is with a call to the
`guess_tempo_parameters`

function:

```
guessed_paramsA = oqupy.guess_tempo_parameters(bath=bath,
start_time=t_start,
end_time=t_end,
tolerance=0.01)
print(guessed_paramsA)
```

```
../oqupy/tempo.py:865: UserWarning: Estimating TEMPO parameters. No guarantee subsequent dynamics calculations are converged. Please refer to the TEMPO documentation and check convergence by varying the parameters manually.
warnings.warn(GUESS_WARNING_MSG, UserWarning)
```

```
----------------------------------------------
TempoParameters object: Roughly estimated parameters
Estimated with 'guess_tempo_parameters()' based on bath correlations.
dt = 0.125
tcut [dkmax] = 2.5 [20]
epsrel = 6.903e-05
add_correlation_time = None
```

As indicated in the description of this object, the parameters were
estimated by analysing the correlations of `bath`

, which are discussed
further below.

From the suggested parameters, we focus on `tcut`

first, assuming the
values of `dt`

and `epsrel`

are reasonable to work with. To do so we
compare results at the recommend `tcut`

to those calculated at a
smaller (`1.25`

) and larger (`5.0`

) values of this parameter,
starting from the spin-up state:

```
initial_state = oqupy.operators.spin_dm('z+')
for tcut in [1.25,2.5,5.0]:
# Create TempoParameters object matching those guessed above, except possibly for tcut
params = oqupy.TempoParameters(dt=0.125, epsrel=6.9e-06, tcut=tcut)
dynamics = oqupy.tempo_compute(system=system,
bath=bath,
initial_state=initial_state,
start_time=t_start,
end_time=t_end,
parameters=params)
t, s_z = dynamics.expectations(sigma_z, real=True)
plt.plot(t, s_z, label=r'${}$'.format(tcut))
plt.xlabel(r'$t$')
plt.ylabel(r'$\langle\sigma_z\rangle$')
plt.legend(title=r'$t_{cut}$')
```

```
--> TEMPO computation:
100.0% 40 of 40 [########################################] 00:00:00
Elapsed time: 0.8s
--> TEMPO computation:
100.0% 40 of 40 [########################################] 00:00:01
Elapsed time: 1.6s
--> TEMPO computation:
100.0% 40 of 40 [########################################] 00:00:01
Elapsed time: 1.9s
```

```
<matplotlib.legend.Legend at 0x76ef12d393d0>
```

We see that `tcut=2.5`

(orange) does very well, matching `tcut=5.0`

(green) until essentially the end of the simulation (the precision
`epsrel`

could well be causing the small discrepancy). We know
`tcut=5.0`

should capture the actual result, because
`tcut=5.0=t_end`

means no memory cutoff was made! In general it is not
always necessary to make a finite memory approximation. For example,
perhaps one is interested in short-time dynamics only. The memory cutoff
can be disable by setting `tcut=None`

; be aware computation to long
times (i.e. many hundreds of timesteps) may then be infeasible.

The `tcut=1.25`

result matches the other two exactly until `t=1.25`

(no memory approximation is made before this time), but deviates shorlty
after. On the other hand, the cost of using the larger `tcut=2.5`

was
a longer computation: 1.6s vs 0.8s above. This was a trivial example,
but in many real calculations the runtimes will be far longer
e.g. minutes or hours. It may be that an intermediary value
`1.25<=tcut<=2.5`

provides a satisfactory approximation - depending on
the desired precision - with a more favourable cost: a TEMPO (or
PT-TEMPO) computation scales **linearly** with the number of steps
included in the memory cutoff.

#### A word of warning

`guess_tempo_parameters`

provides a reasonable starting point for many
cases, but it is only a guess. You should always verify results using a
larger `tcut`

, whilst also not discounting smaller `tcut`

to reduce
the computational requirements. Similar will apply to checking
convergence under `dt`

and `epsrel`

.

Also, note we only inspected the expectations
\(\langle \sigma_z \rangle\). To be most thorough all unique
components of the state matrix should be checked, or at least the
expectations of observables you are intending to study. So, if you were
interested in the coherences as well as the populations, you would want
to add calls to calculate \(\langle \sigma_x \rangle\),
\(\langle \sigma_y \rangle\) above (you can check `tcut=2.5`

is
still good for the above example).

### Discussion - environment correlations

So what influences the required `tcut`

? The physically relevant
timescale is that for the decay of correlations in the environment.
These can be inspected using
`oqupy.helpers.plot_correlations_with_parameters`

:

```
fig, ax = plt.subplots()
params = oqupy.TempoParameters(dt=0.125, epsrel=1, tcut=2.5) # N.B. epsrel not used by helper, and tcut only to set plot t-limits
oqupy.helpers.plot_correlations_with_parameters(bath.correlations, params, ax=ax)
```

```
<AxesSubplot:>
```

This shows the real and imaginary parts of the bath autocorrelation
function, with markers indicating samples of spacing `dt`

. We see that
correlations have not fully decayed by `t=1.25`

, but have - at least
by eye - by `t=2.5`

. It seems like `tcut`

around this value would
indeed be a good choice.

The autocorrelation function depends on the properties of the bath: the
form the spectral density, the cutoff, and the temperature. These are
accounted for by the `guess_tempo_parameters`

function, which is
really analysing the error in performing integrals of this function. The
`tolerance`

parameter specifies the maximum absolute error permitted,
with an inbuilt default value of `3.9e-3`

- passing `tolerance=0.01`

made for slightly ‘easier’ parameters.

Note, however, what is observed in the *system dynamics* also depends
the bath coupling operator and strength (`alpha`

), and that these are
*not* taken into account by the guessing function. More generally, the
nature of the intrinsic system dynamics (see below) and initial state
preparation also has to be considered.

Finally, the guessing function uses specified `start_time`

and
`end_time`

to come up with parameters providing a manageable
computation time over a timescale `end_time-start_time`

, so make sure
to set these to reflect those you actually intend to use in
calculations.

## Choosing dt and epsrel

### Example - convergence for a spin boson model

Continuing with the previous example, we now investigate changing `dt`

at our chosen `tcut=2.5`

.

```
plt.figure(figsize=(10,8))
for dt in [0.0625, 0.125, 0.25]:
params = oqupy.TempoParameters(dt=dt, epsrel=6.9e-05, tcut=2.5)
dynamics = oqupy.tempo_compute(system=system,
bath=bath,
initial_state=initial_state,
start_time=t_start,
end_time=t_end,
parameters=params)
t, s_z = dynamics.expectations(sigma_z, real=True)
plt.plot(t, s_z, label=r'${}$'.format(dt))
plt.xlabel(r'$t$')
plt.ylabel(r'$\langle\sigma_z\rangle$')
plt.legend(title=r'$dt$')
```

```
--> TEMPO computation:
100.0% 80 of 80 [########################################] 00:00:03
Elapsed time: 3.0s
--> TEMPO computation:
100.0% 40 of 40 [########################################] 00:00:00
Elapsed time: 0.9s
--> TEMPO computation:
100.0% 20 of 20 [########################################] 00:00:00
Elapsed time: 0.3s
```

```
<matplotlib.legend.Legend at 0x76ef12cdf190>
```

That doesn’t look good! If we had just checked `dt=0.25`

and
`dt=0.125`

we may have been happy with the convergence, but a halving
of the timestep gave very different results (you can check `dt=0.0625`

is even worse).

The catch here is that we used the same precision `epsrel=6.9e-05`

for
all runs, but `dt=0.125`

requires a smaller `epsrel`

: halving the
timestep *doubles* the number of steps `dkmax`

for which singular
value truncations are made in the bath’s memory `tcut=dt*dkmax`

.

Let’s repeat the calculation with a smaller `epsrel`

at `dt=0.125`

:

```
for dt, epsrel in zip([0.0625,0.125, 0.25],[6.9e-06,6.9e-05,6.9e-05]):
params = oqupy.TempoParameters(dt=dt, epsrel=epsrel, tcut=2.5)
dynamics = oqupy.tempo_compute(system=system,
bath=bath,
initial_state=initial_state,
start_time=t_start,
end_time=t_end,
parameters=params)
t, s_z = dynamics.expectations(sigma_z, real=True)
plt.plot(t, s_z, label=r'${}$, ${:.2g}$'.format(dt,epsrel))
plt.xlabel(r'$t$')
plt.ylabel(r'$\langle\sigma_z\rangle$')
plt.legend(title=r'$dt, \epsilon_{rel}$')
```

```
--> TEMPO computation:
100.0% 80 of 80 [########################################] 00:00:04
Elapsed time: 5.0s
--> TEMPO computation:
100.0% 40 of 40 [########################################] 00:00:00
Elapsed time: 0.9s
--> TEMPO computation:
100.0% 20 of 20 [########################################] 00:00:00
Elapsed time: 0.2s
```

```
<matplotlib.legend.Legend at 0x76ef12a6b1d0>
```

That looks far better. The lesson here is that one cannot expect to be
able to decrease `dt`

at fixed `tcut`

without also decreasing
`epsrel`

. A heuristic used by `guess_tempo_parameters`

, which you
may find useful, is

where tol is a target tolerance (e.g. `tolerance=0.01`

above) and
`dkmax`

the number of steps such that `tcut=dt*dkmax`

.

Note `TempoParameters`

allows the memory cutoff to be specified as the
integer `dkmax`

rather than float `tcut`

, meaning this estimation of
`epsrel`

doesn’t change with `dt`

. However, the author prefers
working at a constant `tcut`

which is set physically by the decay of
correlations in the environment; then one only has to worry about the
simultaneous convergence of `dt`

and `epsrel`

.

Comparing the simulation times at `dt=0.0625`

between the previous two
sets of results, we see the cost of a smaller `epsrel`

is a longer
computation (5 vs. 3 seconds). The time complexity of the singular value
decompositions in the TEMPO tensor network scales with the **third
power** of the internal bond dimension, which is directly controlled by
the precision, so be aware that decreasing `epsrel`

may lead to rapid
increase in computation times.

The last results suggest that we may well already have convergence w.r.t
`epsrel`

at `dt=0.125`

. This should be checked:

```
for epsrel in [6.9e-06,6.9e-05,6.9e-04]:
params = oqupy.TempoParameters(dt=dt, epsrel=epsrel, tcut=2.5)
dynamics = oqupy.tempo_compute(system=system,
bath=bath,
initial_state=initial_state,
start_time=t_start,
end_time=t_end,
parameters=params)
t, s_z = dynamics.expectations(sigma_z, real=True)
plt.plot(t, s_z, label=r'${:.2g}$'.format(epsrel))
plt.xlabel(r'$t$')
plt.ylabel(r'$\langle\sigma_z\rangle$')
plt.legend(title=r'$\epsilon_{rel}$')
```

```
--> TEMPO computation:
100.0% 20 of 20 [########################################] 00:00:00
Elapsed time: 0.3s
--> TEMPO computation:
100.0% 20 of 20 [########################################] 00:00:00
Elapsed time: 0.2s
--> TEMPO computation:
100.0% 20 of 20 [########################################] 00:00:00
Elapsed time: 0.2s
```

```
<matplotlib.legend.Legend at 0x76ef12aec9d0>
```

In summary, we may well be happy with the parameters `dt=0.125`

,
`epsrel=6.9e-05`

, `tcut=2.5`

for this model (we could probably use a
larger `epsrel`

, but the computation is so inexpensive in this example
it is hardly necessary).

So far we have discussed mainly how the environment - namely the memory length - dictates the parameters. We now look at what influence the system can have.

### Resolving fast system dynamics

In the above you may have noticed that the results at `dt=0.125`

,
while converged, were slightly undersampled. This becomes more
noticeable if the scale of the system energies is increased (here by a
factor of 4):

```
Omega = 8.0 # From 2.0
Gamma = 0.08 # From 0.02
system = oqupy.System(0.5 * Omega * sigma_x,
gammas=[Gamma],
lindblad_operators=[sigma_m])
params = oqupy.guess_tempo_parameters(bath=bath,
start_time=t_start,
end_time=t_end,
tolerance=0.01)
print(params)
dynamics = oqupy.tempo_compute(system=system,
bath=bath,
initial_state=initial_state,
start_time=t_start,
end_time=t_end,
parameters=params)
t, s_z = dynamics.expectations(sigma_z, real=True)
plt.plot(t, s_z)
plt.xlabel(r'$t$')
plt.ylabel(r'$\langle\sigma_z\rangle$')
plt.show()
```

```
../oqupy/tempo.py:865: UserWarning: Estimating TEMPO parameters. No guarantee subsequent dynamics calculations are converged. Please refer to the TEMPO documentation and check convergence by varying the parameters manually.
warnings.warn(GUESS_WARNING_MSG, UserWarning)
```

```
----------------------------------------------
TempoParameters object: Roughly estimated parameters
Estimated with 'guess_tempo_parameters()' based on bath correlations.
dt = 0.125
tcut [dkmax] = 2.5 [20]
epsrel = 6.903e-05
add_correlation_time = None
--> TEMPO computation:
100.0% 40 of 40 [########################################] 00:00:03
Elapsed time: 3.5s
```

The call to `guess_tempo_parameters`

returned the same set
`dt=0.125`

, `epsrel=6.9e-05`

, `tcut=2.5`

as before, because it did
not use any information of the system. We can change this, and hopefully
resolve the system dynamics on a more appropriate grid, by providing the
system as an optional argument:

[Warning: long computation]

```
Omega = 8.0 # From 2.0
Gamma = 0.08 # From 0.02
system = oqupy.System(0.5 * Omega * sigma_x,
gammas=[Gamma],
lindblad_operators=[sigma_m])
params = oqupy.guess_tempo_parameters(system=system, # new system argument (optional)
bath=bath,
start_time=t_start,
end_time=t_end,
tolerance=0.01)
print(params)
dynamics = oqupy.tempo_compute(system=system,
bath=bath,
initial_state=initial_state,
start_time=t_start,
end_time=t_end,
parameters=params)
t, s_z = dynamics.expectations(sigma_z, real=True)
plt.plot(t, s_z)
plt.xlabel(r'$t$')
plt.ylabel(r'$\langle\sigma_z\rangle$')
```

```
../oqupy/tempo.py:865: UserWarning: Estimating TEMPO parameters. No guarantee subsequent dynamics calculations are converged. Please refer to the TEMPO documentation and check convergence by varying the parameters manually.
warnings.warn(GUESS_WARNING_MSG, UserWarning)
```

```
----------------------------------------------
TempoParameters object: Roughly estimated parameters
Estimated with 'guess_tempo_parameters()' based on bath correlations and system frequencies (limiting).
dt = 0.03125
tcut [dkmax] = 2.5 [80]
epsrel = 6.903e-06
add_correlation_time = None
--> TEMPO computation:
100.0% 160 of 160 [########################################] 00:01:09
Elapsed time: 69.5s
```

Text(0, 0.5, '$\langle\sigma_z\rangle$')

As both `dkmax`

increased and `epsrel`

decreased to accommodate the
smaller `dt=0.03125`

, the computation took far longer - over a minute
compared to a few seconds at `dt=0.125`

(it may now be worth
investigating whether a larger `epsrel`

can be used).

With a `system`

argument, `guess_tempo_parameters`

uses the matrix
norm of the system Hamiltonian and any Lindblad operators/rates to
estimate a suitable timestep on which to resolve the system dynamics.
This is compared to the `dt`

required to meet the tolerance on error
for the bath correlations, and the smaller of the two is returned. The
description of the `TempoParameters`

object indicates which part was
‘limiting’ i.e. required the smaller `dt`

.

Often it is not necessary to calculate the system dynamics on such a
fine grid. For example, if one only needs to calculate the steady-state
polarisation. Moreover, the undersampling is easy to spot and adjust by
eye. Hence you may choose to not pass a `system`

object to
`guess_tempo_parameters`

. However, note there are cases where not
accounting for system frequencies can lead to more physical features
being missed, namely when the Hamiltonian or Lindblad operators/rates
are (rapidly) *time-dependent.*

#### What sets dt, really?

The main error associated with `dt`

is that from the Trotter splitting
of the system propagators. In a simple (non-symmetrised) splitting, a
basic requirement is

In words: error arises from non-commutation between the system and bath coupling operator. This simply reflects the fact that in the discretisation of the path integral the splitting is made between the system and environment Hamiltonians. In cases where \(H_S\) commutes with \(H_E\), such as the independent boson model, \(dt\) can be arbitrarily large without physical error.

Note `guess_tempo_parameters`

does *not* attempt to estimate the
Trotter error, even when both `system`

and `bath`

objects are
specified - another reason to be cautious when using the estimate
produced by this function.

## Further considerations

### Additional TempoParameters arguments

For completeness, there are a few other parameters that can be passed to
the `TempoParameters`

constructor: - `subdiv_limit`

and
`liouvillian_epsrel`

. These control the maximum number of subdivisions
and relative error tolerance when integrating a time-dependent system
Liouvillian to construct the system propagators. It is unlikely you will
need to change them from their default values (`265`

, `2**(-26)`

) -
`add_correlation_time`

. This allows one to include correlations
*beyond* `tcut`

in the final bath tensor at `dkmax`

(i.e., have your
finite memory cutoff cake and eat it). Doing so may provide better
approximations in cases where `tcut`

is limited due to computational
difficultly. See
[Strathearn2017] for
details.

### Bath coupling degeneracies

The bath tensors in the TEMPO network are nominally \(d^2\times d^2\) matrices, where \(d\) is the system Hilbert space dimension. If there are degeneracies in the sums or differences of the eigenvalues of the system operator coupling to the environment, it is possible for the dimension of these tensors to be reduced.

Specifying `unique=True`

as an argument to `oqupy.tempo_compute`

,
degeneracy checking is enabled: if a dimensional reduction of the bath
tensors is possible, the lower dimensional tensors will be used. We
expect this to provide in many cases a significant speed-up without any
significant loss of accuracy, although this has not been extensively
tested (new in `v0.5.0`

).

### Mean-field systems

For calculating mean-field dynamics, there is an additional requirement
on `dt`

being small enough so not as to introduce error when
integrating the field equation of motion between timesteps (2nd order
Runge-Kutta). Common experience is that this is generally satisfied if
`dt`

is sufficiently small to avoid Trotter error. Still, you will
want to at least inspect the field dynamics in addition to the system
observables when checking convergence.

Note that, for the purposes of estimating the characteristic frequency
of a `SystemWithField`

object, `guess_tempo_parameters`

uses an
arbitrary complex value \(\exp(i \pi/4)\) for the field variable
when evaluating the system Hamiltonian. This may give a poor estimation
for situations where the field variable is not of order \(1\) in the
dynamics.

## PT-TEMPO

The above considerations for `tcut`

, `dt`

and `epsrel`

all apply
to a PT-TEMPO computation, with the following caveats:

Convergence for a TEMPO computation does not necessarily imply convergence for a PT-TEMPO computation. This is important as it is often convenient to perform several one-off TEMPO computations to determine a good set of computational parameters to save having to construct many large process tensors. You can still take this approach, but must be sure to check for convergence in the PT-TEMPO computation when you have derived a reasonable set.

Similar to 1., the best possible accuracy of a TEMPO and PT-TEMPO computation may be different. In particular, there may be a trade-off of accuracy for overall reduced computation time when using the PT approach. We note that the error in PT-TEMPO has also been observed ([FowlerWright2022]) to become unstable at very high precisions (small

`epsrel`

), i.e., there may be a higher floor for how small you can make`epsrel`

.The computational difficultly of constructing the PT may not be monotonic with memory cutoff

`dkmax`

(or`tcut`

). In particular, the computational time may diverge*below*a certain`dkmax`

, a.k.a, the ‘dkmax anomaly’. We highlight this counter-intuitive behaviour, which seems to occur at relatively high precisions (small`epsrel`

) with short timesteps, because it may lead one to falsely believe the computation of a process tensor is not feasible. See below for a demonstration and the Supplementary Material of [FowlerWright2022] for further discussion.

### The dkmax anomaly

We consider constructing a process tensor of 250 timesteps for the
harmonic environment discussed in the Mean-Field
Dynamics
tutorial, but with a smaller timestep `dt=1e-3`

ps and `epsrel=1e-8`

than considered there. Note that the following computations are very
demanding.

```
alpha = 0.25 # Doesn't affect PT computation
nu_c = 227.9
T = 39.3
start_time = 0.0
dt = 1e-3
epsrel = 1e-8
end_time = 250 * dt # 251 steps starting from t=0.0
correlations = oqupy.PowerLawSD(alpha=alpha,
zeta=1,
cutoff=nu_c,
cutoff_type='gaussian')
bath = oqupy.Bath(oqupy.operators.sigma("z")/2.0, correlations)
```

We firstly set `dkmax=250`

(or `None`

), i.e., make no memory
approximation:

```
params = oqupy.TempoParameters(dt=dt,
epsrel=epsrel,
dkmax=250)
process_tensor = oqupy.pt_tempo_compute(bath=bath,
start_time=start_time,
end_time=end_time,
parameters=params)
```

```
--> PT-TEMPO computation:
100.0% 250 of 250 [########################################] 00:01:37
Elapsed time: 97.3s
```

Including the full memory didn’t take too long, just over one and a half minutes on a modern desktop (AMD 6-Core processor @4.7GHz, python3.6).

What about if we now make a memory approximation, say only after
`dkmax=225`

timesteps:

```
params = oqupy.TempoParameters(dt=dt,
epsrel=epsrel,
dkmax=225)
process_tensor = oqupy.pt_tempo_compute(bath=bath,
start_time=start_time,
end_time=end_time,
parameters=params)
```

```
--> PT-TEMPO computation:
100.0% 250 of 250 [########################################] 00:08:04
Elapsed time: 484.6s
```

That was far slower (8 mins)! You can try `dkmax=200`

- on my hardware
the computation took half an hour; it may not be possible to complete
the calculation `dkmax`

much below this.

In general, there may exist some `dkmax`

value (here close to 250)
below which the computational time grows quickly. On the other hand, for
longer computations, e.g. a 500 step process tensor, increases of
`dkmax`

will eventually lead to increasing compute times again,
although the dynamics will surely be converged with respect to this
parameter well before then.

The take-home message is to not discount that a stalling PT computation
may in fact be possible with an increase in the memory length. In these
situations one approach is to start with `dkmax=None`

and work
backwards to find the `dkmax`

offering the minimum compute time
(before the rapid increase).