MATLAB: Faber approximation of the matrix exponential to solve Liouville-von-Neumann equation

Question

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:

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.

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
            break
        else
        M = M+1;
        end
    end

    % 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);
    end

    % 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);
    end
    z = temp;
end

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

Answer 1

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

Concept

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.

Implementation

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.