EDIT: I moved the full code to my github page so the post can be read more easily.

I am writing a script to take the Faber approximation approach outlined in Hassan Fahs paper (free access) and apply it to the Liouville-von-Neumann equation to propagate the density matrix $\rho$ with some inspiration from the scripts found here for general structure.

Currently the approach is stuck with the approximation of the exponential. The idea behind the propagation scheme is to compute the update step as: $$\rho_{n+1}(\tau) = exp({\mathcal{L}\tau})\rho_n = \sum_{m=0}^{\infty}c_m(\tau) F_m(\mathcal{L})\rho_n.$$ As we can see the Faber polynomials are matrix valued, and perhaps -if you have taken a look into the script or the paper-you will see that the the spectrum of $\mathcal{L}$ has to be scaled with a factor I will simply refer to as $c$ in order to preserve stability, which in turn means the time step has to be adapted to $c*\tau =\tilde{\tau}$ So when we approximate the exponential we can see that our propagation scheme becomes: $$\rho_{n+1}(\tau*c)=\rho_{n+1}(\tilde{\tau}) = exp({\frac{\mathcal{L}}{c}\tau c}\rho_n) = \sum_{m=0}^{\infty}c_m(\tilde{\tau}) F_m(\frac{\mathcal{L}}{c})\rho_n.$$

The matrix valued Faber polynomial recurrence relation for an elliptic domain in my case is: $$F_{m+1}(\mathcal{L})\rho(0) = (\mathcal{L}-b_0\mathcal{I})F_m(\mathcal{L})\rho(0) - b_1 F_{m-1}(\mathcal{L})\rho(0).$$ $\mathcal{L}$ denotes the Liouville superoperator. Obviously, $\mathcal{L}$ is a function of the density matrix.

After initialising a basic density matrix and computing the factors as in the paper above I tried to implement the algorithm as outlined by Dr . Fahs:

enter image description here

But before trying to compute the full scheme I simply wanted to see if the approximation of the exponential $$exp(\mathcal{L}_{sc}\tilde{\tau})= \sum_{m=0}^{\infty}c_m(\tilde{\tau}) F_m(\mathcal{L}_{sc})$$ was correct. In order to do so tried to compare my solution for 100 iterations with the standard diagonalisation approach.

The resulting picture is created by plotting the third line of a matrix whose columns are filled with the main diagonal of the respective matrix. The blue line is the diagonalisation approach, the orange one mine.

enter image description here

Oddly, my approximation stays arround 1 for the entirety of the execution of the script.

I've traced the issue to the polynomial recurrence relation: $F_1...F_m$ all depend on the factors of the Laurent expansion $b_0, b_1$ and the scaled Liouvillian $\mathcal{L}_sc$. They barely change, making me wonder if the algorithm for the approximation mentioned above is correct.

Here is the code snippet for the computation of the Faber coefficients and Matrix valued polynomial recurrence relation, as well as the summing up of all the factors:

% Faber approximation of the matrix exponential
function z = faber(L, dt, N)

    % Faber coefficients for ellipic domain
    c_m = @(dt_tilde, m) (((-1i/sqrt(b_1))^m)*exp(dt_tilde*b_0)*...
    besselj(m, 2*dt_tilde*sqrt(-b_1)));

    % establish polynomial truncation order so that M >e*sf*dt
    M = 1;
    while true
        orderinterator(M) = double(c_m(dt_tilde, M));
        if abs(orderinterator(M))<10e-15
        M = M+1;

    % Compute time dependent Faber coefficients of order M for elliptic domain
    % initialize CM
    CM = zeros(M+1,1);
    for m = 0:M
        CM(m+1) = c_m(dt_tilde, m);

    % Compute matrix valued polynomials with initial value density matrix from
    % previous iteration
    I = ones(size(L_sc));
    P = zeros((M+1)*N, N);

    % Compute matrix valued Faber polynomial recurrence relation
    % Compute initial value polynomials P_0, P_1, P_2:
    % F_0
    P(1:N, 1:N) = I;

    % F_1
    P((N+1):(2*N), 1:N) = (L_sc-b_0*I);

    % F_2
    P((2*N+1):(3*N), 1:N) = (L_sc-b_0*I)-2*b_1*I;

    % c0*F_0+c1*F_1+c_2*F_2
    temp = CM(1)*P(1:N, 1:N)+CM(2)*P(N+1:2*N, 1:N)+CM(3)*P(2*N+1:3*N, 1:N); 

    % F_3 ... F_M = c2*F_2...cm*F_m
    for i = 3:M
        % F_{m+1} = L_sc(rho_fab)*F_m - b_0*F_m - b_1*F_{m-1}
        P((i*N)+1:(i+1)*N, 1:N) = L_sc*P((i-1)*N+1:(i)*N, 1:N)...
            -b_0*P((i-1)*N+1:(i)*N, 1:N)...
            -b_1*P((i-1-1)*N+1:(i-1)*N, 1:N);

        temp = temp + CM(i+1)*P((i*N)+1:(i+1)*N, 1:N);
    z = temp;

Can someone perhaps elaborate on where it went wrong? Have I perhaps misunderstood the approach ?

  • $\begingroup$ Welcome to SciComp.SE. It is not clear to me what your question is. $\endgroup$
    – nicoguaro
    Commented Mar 12, 2019 at 14:01
  • $\begingroup$ Dear nicoguaro, I edited my question. The error is apparently somehow related to indexing/how I computed the recurrence relation, specifically why only the first column is being computed and not the rest. $\endgroup$
    – Tony_V
    Commented Mar 12, 2019 at 14:48
  • $\begingroup$ I just edited again. The issue truly makes me wonder if this algorithm works at all. Yet, it has been utilized by almost half a dozen authors. Any input on the topic is welcome ! $\endgroup$
    – Tony_V
    Commented Mar 13, 2019 at 15:56
  • $\begingroup$ I think you mixed up Liouvillian and Hamiltonian. The algorithm calculates the exponential of the Hamiltonian, your code uses the Liouvillian of the initial density matrix. What do you want to calculate? $\endgroup$ Commented Mar 15, 2019 at 9:33
  • $\begingroup$ I don't see a sparse matrix in your code?! $\endgroup$ Commented Mar 15, 2019 at 14:28

1 Answer 1


Not sure whether this provides a full answer, but at least it provides some thoughts and hints.


I think the expression $\exp(\mathcal{Lt})\rho$ needs clarification. This is not necessarily a matrix exponential, but rather an abstract notation for the solution operator. For ODEs $\partial_t x = A x$, where $x$ is a vector and $A$ is a matrix (such as in the Fahs reference), the solution operator is indeed a matrix exponential. On the other hand, for the Liouville-von Neumann equation the solution operator is (roughly and without any relaxation terms) $\exp(\mathcal{L}t)\rho = U \rho U^\dagger$, where $U = \exp(-\mathrm{i}\hbar^{-1} H t)$ and $H$ is the Hamiltonian. In your code, $U$ is calculated using MATLABs built-in expm function and may serve as reference.

Now I am not sure what your plan is, but you could try to calculate $U$ using your Faber approximation code and check whether it gives the same result as the expm method.


I tried to do this using your code, but the code returned errors. If you write a function that performs a certain calculation, immediately think about ways to test it. Even during early stages of development this will be helpful to detect mistakes.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.