Time propagation: error estimation of the Siegert states expansions

The time_propagation_initial_momentum.ipynb notebook showed the overcompleteness of the Siegert states basis set, where we saw that the Berggren expansion and the one presented in Santra et al., *PRA* **71** (2005) give almost similar results.

The goal of this notebook is to assert the quality of different Siegert states expansion when it comes to evaluate the time-propagation of an initial wavepacket. The same initial wave-packets as in the previous notebook will be used: the interest of the present notebook is only to quantify what was presented in the previous one.

Three different Siegert states expansions will be considered: the Mittag-Leffler expansion (MLE), the Berggren expansion and the exact Siegert states expansion (see notebook time_propagation.ipynb for more details on these expansions). They will be compared to the exact time-propagation using the continuum states.

Initialization

Import useful modules and classes

[ ]:
# Make the notebook aware of some of the SiegPy module classes
from siegpy import SWPBasisSet, Gaussian, SWPotential
# Other imports
import numpy as np

Define functions to estimate and show the errors

The absolute and relative errors of different Siegert states expansions with respect to the reference result using bound and continuum states are computed for various times.

[ ]:
def print_errors(exact, siegert_exp, dx):
    """
    Function printing the absolute and relative error of a given
    Siegert states expansion of the time propagation of a wave-packet
    for different times.

    :param exact: Exact time-propagation of a wavepacket.
    :type exact: 2D numpy.array
    :param siegert_exp: Siegert states expansion of the time-
                        propagation of the same wavepacket.
    :type siegert_exp: 2D numpy.array
    :param dx: Space grid-step.
    :type: float
    """
    # Print lines to initialize the table
    title = "  t    abs. error   rel. error"
    line = "-"*len(title)
    print(line)
    print(title)
    print(line)
    # Get the error estimations for each time
    for i, t in enumerate(time_grid):
        abs_err = np.trapz(np.abs(exact[i] - siegert_exp[i]), dx=dx)
        norm = np.trapz(np.abs(exact[i]), dx=dx)
        #abs_err = np.trapz(np.abs(exact[i] - siegert_exp[i])**2, dx=dx)
        #norm = np.trapz(np.abs(exact[i])**2, dx=dx)
        rel_err = abs_err / norm
        print("{:.2f}   {: .3e}   {: .3e}".format(t, abs_err, rel_err))
    # Print lines to finalize the table
    print(line)
    print()


def SSE_errors(filename, nres, k_max, h_k, nx, test_func, time_grid,
               exact=True, MLE=True, Ber=True):
    """
    Function performing the time-propagation of an initial wave-packet
    using different Siegert states expansions and printing the absolute
    and relative errors when compared to the exact time-propagation of
    this initial wavepacket.

    :param filename: Name of the file where to look for the Siegert states.
    :type filename: str
    :param nres: Number of resonant states to use in the Siegert states expansions.
    :type nres: int
    :param k_max: Wavenumber of the last continuum state involved in the exact calculations.
    :type k_max: float
    :param h_k: Wavenumber grid step of the continuum basis set.
    :type h_k: float
    :param nx: Number of space grid points.
    :type nx: int
    :param test_func: Initial wave-packet.
    :type test_func: Wavefunction
    :param time_grid: Times at which the
    :type time_grid: list or numpy.array
    :param choice: Optional, used to say which Siegert state expansion is desired.
    :type choice: dict
    """
    # Siegert states basis set initialized from a file
    siegerts = SWPBasisSet.from_file(filename, nres=nres)

    # Read the potential from a data file
    potential = siegerts[0].potential
    l = potential.width

    # Define a grid space
    xgrid = np.linspace(-l/2, l/2, nx)
    dx = xgrid[1] - xgrid[0]

    # Discretize the Siegert states over the grid
    for s in siegerts:
        s.grid = xgrid

    # Exact basis set
    cont = SWPBasisSet.find_continuum_states(potential, k_max, h_k, grid=xgrid)
    exact = cont + siegerts.bounds

    # Evaluation of the exact time propagation
    exact_tp = exact.exact_propagation(test_func, time_grid)

    if exact:
        # Evaluation of the exact Siegert expansion of the time propagation
        exact_S_tp = siegerts.exact_Siegert_propagation(test_func, time_grid)
        # Find the error with respect ot the exact time propagation
        print("         Exact Siegert")
        print_errors(exact_tp, exact_S_tp, dx)

    if MLE:
        # Evaluation of the Mittag-Leffler expansion of the time propagation
        MLE_tp = siegerts.MLE_propagation(test_func, time_grid)
        # Find the error with respect ot the exact time propagation
        print("              MLE")
        print_errors(exact_tp, MLE_tp, dx)

    if Ber:
        # Evaluation of the Berggren expansion of the time propagation
        Ber_tp = siegerts.Berggren_propagation(test_func, time_grid)
        # Find the error with respect ot the exact time propagation
        print("            Berggren")
        print_errors(exact_tp, Ber_tp, dx)

Define the filename of the data file

[ ]:
# Siegert states basis set initialized from a file
filename = 'siegerts.dat'
siegerts = SWPBasisSet.from_file(filename)

Case 1: initial wave-packet with a small initial momentum

The first initial wave-packet considered uses a small initial momentum.

Define the test function

[ ]:
l = siegerts.potential.width  # width of the potential
sigma = l/20.  # width of the Gaussian
x_c   = -l/4.  # center of the Gaussian
k_0   = 1.  # initial momentum
gauss_small_momentum = Gaussian(sigma, x_c, k0=k_0)

Estimation of the errors of the Siegert states expansions

As a starting point, let us use parameters allowing for a moderate accuracy of the error measurements. The Siegert states and the exact basis sets have a close extension (\(|k_{res, max}| \approx k_{max}\)) and the space grid is made of a limited number of points.

[5]:
# Definition of the time grid
time_grid = [0.0, 0.25, 0.5, 0.75, 1.0, 2.0, 3.0]

SSE_errors(filename, 50, 40, 0.01, 201, gauss_small_momentum, time_grid)
         Exact Siegert
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    8.722e-06    1.566e-05
0.25    6.219e-07    5.567e-07
0.50    4.425e-07    3.656e-07
0.75    3.630e-07    3.132e-07
1.00    3.156e-07    2.864e-07
2.00    2.260e-07    2.090e-07
3.00    1.863e-07    1.684e-07
------------------------------

              MLE
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    8.722e-06    1.566e-05
0.25    6.329e-01    5.665e-01
0.50    2.953e+03    2.441e+03
0.75    9.029e+08    7.791e+08
1.00    3.302e+14    2.996e+14
2.00    5.546e+36    5.130e+36
3.00    9.162e+58    8.281e+58
------------------------------

            Berggren
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    1.575e-01    2.828e-01
0.25    7.557e-02    6.764e-02
0.50    4.717e-02    3.898e-02
0.75    3.292e-02    2.841e-02
1.00    2.455e-02    2.228e-02
2.00    1.082e-02    1.001e-02
3.00    6.360e-03    5.749e-03
------------------------------

Let us see the influence of using denser wavenumber and space grids, that should give more accurate results:

[6]:
SSE_errors(filename, 50, 40, 0.005, 1601, gauss_small_momentum, time_grid)
         Exact Siegert
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    8.562e-06    1.538e-05
0.25    6.024e-07    5.392e-07
0.50    4.262e-07    3.522e-07
0.75    3.481e-07    3.004e-07
1.00    3.016e-07    2.736e-07
2.00    2.136e-07    1.976e-07
3.00    1.746e-07    1.578e-07
------------------------------

              MLE
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    8.562e-06    1.538e-05
0.25    6.329e-01    5.665e-01
0.50    2.953e+03    2.440e+03
0.75    9.028e+08    7.791e+08
1.00    3.302e+14    2.996e+14
2.00    5.545e+36    5.130e+36
3.00    9.160e+58    8.280e+58
------------------------------

            Berggren
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    1.574e-01    2.828e-01
0.25    7.556e-02    6.764e-02
0.50    4.717e-02    3.897e-02
0.75    3.292e-02    2.840e-02
1.00    2.455e-02    2.228e-02
2.00    1.082e-02    1.001e-02
3.00    6.360e-03    5.749e-03
------------------------------

As you can see, the previous results were almost converged: the evaluation of the erros is only slightly modified by the densification of the grids (especially for the exact Siegert states expansion).

Combining the use of denser wavenumber and space grids grids with an increase of the range of the continuum states wavenumber grid (associated to an increase of the number of resonant couples), the results should be even more accurate:

[7]:
SSE_errors(filename, 78, 60, 0.005, 1601, gauss_small_momentum, time_grid)
         Exact Siegert
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    8.542e-06    1.534e-05
0.25    6.023e-07    5.391e-07
0.50    4.261e-07    3.521e-07
0.75    3.481e-07    3.004e-07
1.00    3.016e-07    2.736e-07
2.00    2.136e-07    1.976e-07
3.00    1.746e-07    1.578e-07
------------------------------

              MLE
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    8.542e-06    1.534e-05
0.25    6.873e+01    6.153e+01
0.50    1.511e+11    1.248e+11
0.75    5.496e+20    4.743e+20
1.00    1.836e+30    1.665e+30
2.00    1.995e+68    1.845e+68
3.00    2.288e+106    2.069e+106
------------------------------

            Berggren
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    1.574e-01    2.828e-01
0.25    7.556e-02    6.764e-02
0.50    4.717e-02    3.897e-02
0.75    3.292e-02    2.840e-02
1.00    2.455e-02    2.228e-02
2.00    1.082e-02    1.001e-02
3.00    6.360e-03    5.749e-03
------------------------------

The errors Berggren and the Siegert states expansion are the same as earlier, showing that the results are converged in terms of the extension of the basis sets. However, the MLE is even more wrong, and diverges even more badly; this is not a surprise, since more divergent terms are added in the sum.

We are now convinced that the error evaluation is sufficiently converged for our purpose. As expected, the results give the same conclusions as in the previous notebook:

  • It is clear that the MLE diverges in the long time limit.
  • The exact Siegert states expansion produces the best result, whatever the time.
  • The errors induced by the use of the Berggren expansion decreases over time, making it a rather good expansion in the long-time limit, but not as good as the exact Siegert states expansion.

Case 2: initial wave-packet with large initial momentum

We saw in the previous notebook that the Berggren expansion was very good for all times in the case of an initial wave-packet with a large initial momentum. What are the relative and absolute errors in that case?

Define the test function

[ ]:
k_0 = 20.  # initial momentum
gauss_large_momentum = Gaussian(sigma, x_c, k0=k_0)

Estimation of the errors of the Siegert states expansions

We start with the same parameters leading to moderate accuracy results:

[9]:
SSE_errors(filename, 50, 40, 0.01, 201, gauss_large_momentum, time_grid)
         Exact Siegert
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    1.441e-04    2.588e-04
0.25    5.216e-06    5.375e-05
0.50    4.181e-06    3.773e-04
0.75    3.731e-06    1.190e-03
1.00    3.338e-06    2.036e-03
2.00    2.074e-06    2.424e-03
3.00    1.822e-06    1.985e-03
------------------------------

              MLE
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    1.441e-04    2.588e-04
0.25    2.464e+02    2.539e+03
0.50    5.505e+06    4.969e+08
0.75    7.641e+11    2.438e+14
1.00    3.037e+17    1.852e+20
2.00    4.225e+39    4.939e+42
3.00    7.028e+61    7.657e+64
------------------------------

            Berggren
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    5.335e-04    9.582e-04
0.25    2.731e-04    2.814e-03
0.50    1.774e-04    1.601e-02
0.75    1.270e-04    4.053e-02
1.00    9.689e-05    5.911e-02
2.00    4.642e-05    5.427e-02
3.00    2.899e-05    3.158e-02
------------------------------

Using denser wavenumber and space grids, we show that the previous results could already be considered as almost converged with respect to those parameters:

[10]:
SSE_errors(filename, 50, 40, 0.005, 1601, gauss_large_momentum, time_grid)
         Exact Siegert
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    1.440e-04    2.586e-04
0.25    5.207e-06    5.365e-05
0.50    4.216e-06    3.806e-04
0.75    3.762e-06    1.201e-03
1.00    3.340e-06    2.037e-03
2.00    2.099e-06    2.454e-03
3.00    1.772e-06    1.931e-03
------------------------------

              MLE
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    1.440e-04    2.586e-04
0.25    2.464e+02    2.539e+03
0.50    5.505e+06    4.969e+08
0.75    7.639e+11    2.438e+14
1.00    3.036e+17    1.852e+20
2.00    4.225e+39    4.939e+42
3.00    7.027e+61    7.656e+64
------------------------------

            Berggren
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    5.332e-04    9.575e-04
0.25    2.731e-04    2.815e-03
0.50    1.774e-04    1.602e-02
0.75    1.270e-04    4.053e-02
1.00    9.690e-05    5.911e-02
2.00    4.643e-05    5.427e-02
3.00    2.899e-05    3.158e-02
------------------------------

The impact of these parameters (wavenumber and space grid densities) is smaller than in the case of an initial wavepacket with a smaller initial momentum. The previous results could already be considered as converged with respect to those parameters.

Finally, increasing the number of states while using the dense wavenumber and space grids lead to the following results:

[11]:
SSE_errors(filename, 78, 60, 0.005, 1601, gauss_large_momentum, time_grid)
         Exact Siegert
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    1.522e-05    2.733e-05
0.25    7.945e-07    8.187e-06
0.50    5.606e-07    5.060e-05
0.75    4.574e-07    1.460e-04
1.00    3.960e-07    2.416e-04
2.00    2.799e-07    3.271e-04
3.00    2.285e-07    2.489e-04
------------------------------

              MLE
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    1.522e-05    2.733e-05
0.25    2.767e+02    2.851e+03
0.50    1.249e+11    1.127e+13
0.75    4.538e+20    1.448e+23
1.00    1.517e+30    9.255e+32
2.00    1.649e+68    1.928e+71
3.00    1.892e+106    2.061e+109
------------------------------

            Berggren
------------------------------
  t    abs. error   rel. error
------------------------------
0.00    4.602e-04    8.264e-04
0.25    2.731e-04    2.815e-03
0.50    1.774e-04    1.602e-02
0.75    1.270e-04    4.052e-02
1.00    9.689e-05    5.910e-02
2.00    4.643e-05    5.428e-02
3.00    2.899e-05    3.158e-02
------------------------------

These parameters impact the quality of the error, especially for short time. This is not suprising because the completeness relation for the intitial wave-packet with a larger initial momentum converges at a larger wavenumber (see the previous notebook): the convergence with respect to the highest wavenumber of the basis set is more difficult to reach.

The results point to the same conclusions as in the first case, mainly showing that even if the Berggren expansion is in very good agreement for all times, the exact Siegert states gives better results.

Conclusion

This concludes our study of the error estimation, showing that the exact Siegert state expansion truly is the most efficient when it comes to reproduce the time propagation of wavepacket, whatever the time and the wavepacket.

For more details on the convergence of the error estimation with respect to the different parameters, you might consider going through the time_propagation_error_estimation_convergence.ipynb notebook.