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

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 ?

• Welcome to SciComp.SE. It is not clear to me what your question is. – nicoguaro Mar 12 at 14:01
• 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. – Tony_V Mar 12 at 14:48
• 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 ! – Tony_V Mar 13 at 15:56
• 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? – carlosvalderrama Mar 15 at 9:33
• I don't see a sparse matrix in your code?! – carlosvalderrama Mar 15 at 14:28

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.