FermiFab User's Guide (Matlab version)

FermiFab

Christian B. Mendl, 2008-2015

Contents

1.1 What is FermiFab?

FermiFab stands for "fermion laboratory" and is a quantum physics toolbox for small numbers of fermions. It is mainly concerned with the representation of (symbolic) fermionic wavefunctions and the calculation of their reduced density matrices (RDMs). The toolbox transparently handles the inherent antisymmetrization of wavefunctions and incorporates the creation and annihilation formalism.

FermiFab is available at http://sourceforge.net/projects/fermifab.

1.2 Fermi states

A fundamental building block of multi-fermion quantum mechanics are Slater determinants, which you can think of as a collection of "orbitals" (or slots), some of which are occupied by a fermionic particle (e.g., an electron):

In this illustration, each circle represents an orbital, and the filled circles denote occupied orbitals. In mathematical terms, the available number of orbitals 'orbs' is the dimension of the underlying single-particle Hilbert space and the number of occupied orbitals is the particle number . Thus there are altogether

Slater determinants. Each quantum mechanical Fermi state is a linear combination of these Slater determinants with complex coefficients. Formally, the quantum Fermi state is an element of the -particle Fock space .

In FermiFab, you can define, say, a particle state (or wavefunction) and 6 available orbitals with the fermistate command:

orbs = 6; N = 4;
x = zeros(nchoosek(orbs,N),1); x(1)=1/sqrt(2); x(2)=1i/sqrt(2);
psi = fermistate(orbs,N,x)
psi = 
	Fermi State (orbs == 6, N == 4)
	(0.70711)|1234> + (0+0.70711i)|1235>

The vector contains the Slater determinant coefficients of in lexicographical order. Let's assign names to the orbitals of :

psi = set(psi,'orbnames',{'1s' '1s~' '2s' '2s~' '2p' '2p~'})
psi = 
	Fermi State (orbs == 6, N == 4)
	(0.70711)|1s 1s~ 2s 2s~> + (0+0.70711i)|1s 1s~ 2s 2p>

In this example, the names denote electronic subshells of atoms, with the tilde ~ representing spin-down, otherwise spin-up.

The rank-one projector in bra-ket notation

or "density matrix" of can be calculated intuitively by

psi*psi'
ans = 
	Fermi Operator wedge^4 H -> wedge^4 H (orbs == 6)

Matrix representation w.r.t. ordered Slater basis
(|1s 1s~ 2s 2s~>, |1s 1s~ 2s 2p>, ... |2s 2s~ 2p 2p~>) -> (|1s 1s~ 2s 2s~>, |1s 1s~ 2s 2p>, ... |2s 2s~ 2p 2p~>):

  Columns 1 through 4

   0.5000 + 0.0000i   0.0000 - 0.5000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.5000i   0.5000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
...

Note that the result is now an operator acting on the 4-particle antisymmetrized Hilbert space !

1.3 Reduced density matrices (RDMs)

The core feature of the toolbox is the efficient calculation of reduced density matrices. For example, the two-body reduced density matrix defined by

can be calculated by the rdm command as

rdm(psi,2)
ans = 
	Fermi Operator wedge^2 H -> wedge^2 H (orbs == 6)

Matrix representation w.r.t. ordered Slater basis
(|1s 1s~>, |1s 2s>, ... |2p 2p~>) -> (|1s 1s~>, |1s 2s>, ... |2p 2p~>):

  Columns 1 through 4

   1.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   1.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   1.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.5000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.5000i
...

In the above formula, and denote creation and annihilation operators, respectively.

More generally, for wavefunctions and orthonormal basis sets of with and , define the reduced density matrix

by

The significance of this definition can be seen as follows. Any linear map with matrix representation may be "lifted" to an operator by

implemented as p2N in FermiFab. A prominent example is the Coulomb operator with , which describes the pairwise interaction between charged particles. Now, the expectation value with respect to equals

In other words, this equation switches from to for . For many applications, this is the only possibility to avoid the "curse of dimensionality" induced by the large , -particle systems.

1.4 Symbolic computations

FermiFab is compatible with the "Symbolic Math Toolbox". So you may insert symbolic variables as well in the above examples (assuming the Symbolic Math Toolbox is installed)

syms a b c
y = sym(zeros(1,nchoosek(orbs,N)));
y(1) = a; y(3) = 1i*b^2; y(4) = 1/c;
psi = set(psi,'data',y)
rdm(psi,2)
psi = 
	Fermi State (orbs == 6, N == 4)
	(a)|1s 1s~ 2s 2s~> + (b^2*i)|1s 1s~ 2s~ 2p> + (1/c)|1s 2s 2s~ 2p>


ans = 
	Fermi Operator wedge^2 H -> wedge^2 H (orbs == 6)

Matrix representation w.r.t. ordered Slater basis
(|1s 1s~>, |1s 2s>, ... |2p 2p~>) -> (|1s 1s~>, |1s 2s>, ... |2p 2p~>):

[ abs(a)^2 + abs(b)^4,       (b^2*i)/conj(c),              0,                                0,                   0,                     0,             a/conj(c),                0,               0,                     0, 0, 0, 0, 0, 0]
[    -(conj(b)^2*i)/c, abs(a)^2 + 1/abs(c)^2,              0,                                0,                   0,                     0,         a*conj(b)^2*i,                0,               0,                     0, 0, 0, 0, 0, 0]
[                   0,                     0,       abs(a)^2,                                0,                   0,                     0,                     0,    a*conj(b)^2*i,      -a/conj(c),                     0, 0, 0, 0, 0, 0]
[                   0,                     0,              0, abs(a)^2 + abs(b)^4 + 1/abs(c)^2,                   0,                     0,                     0,                0,               0,                     0, 0, 0, 0, 0, 0]
...

1.5 Tensor products of operators

For any linear operator (i.e., matrix) , it holds that

That is, we obtain a matrix representation of acting on . The tensor_op command implements precisely this operation. The following code lines are taken from the "natural orbitals" example in test/norbs.m:

orbs = 6; N = 4;
psi = fermistate(orbs,N,crand(nchoosek(orbs,N),1));
[U,D] = eig(rdm(psi,1));

crand generates pseudorandom complex numbers (similar to rand), and eig computes eigenvalues and -vectors. Thus, the eigenvectors of the one-body reduced density matrix of are stored in . Performing a corresponding base change on using these eigenvectors should result in a diagonal one-body RDM [Löwdin 1955]:

psi = tensor_op(U,N)'*psi;
G = get(rdm(psi,1),'data');
err = norm(G-diag(diag(G)))
err =

   1.0232e-15

In many physical applications, one can take advantage of unitary base changes on such that subsequent computations are simplified. The above code shows how to implement the according base change on .

1.6 State configurations

For performance and memory efficiency reasons, the concept of "configurations" are built into FermiFab, i.e. you can partition the orbitals into several groups, each containing a fixed number of particles. Let's say your system involves 3 particles in 9 orbitals, with exactly 2 particles in the first 5 orbitals and 1 particle in the remaining 4 orbitals. Then a fermistate reflecting this configuration is specified by

orbs = [5,4]; N = [2,1];
psi = fermistate(orbs,N)
psi = 
	Fermi State (orbs == 9, N == 3)
	|126>

Note that is the lexicographically first basis vector respecting the configuration constraints, and that requires only

rather than complex entries:

length(psi)
ans =

    40

The command rdm works transparently for all configurations, so behaves like a standard 9-orbital 3-particle state.

What happens if you add two fermistate wavefunctions with different but compatible configurations (i.e., the total number of orbitals and particles agrees)?

orbs = [2,7]; N = [1,2];
phi = fermistate(orbs,N)
phi = 
	Fermi State (orbs == 9, N == 3)
	|134>

length(phi)
ans =

    42

Adding and gives

chi = psi+phi
chi = 
	Fermi State (orbs == 9, N == 3)
	|134> + |126>

as expected - so how is this accomplished? FermiFab has detected that it needs to merge the two configurations, resulting in the full-fledged 9-orbital 3-particle state, which is reflected by

length(chi)
ans =

    84

1.7 References