Ab Initio Simulations of Condensed Matter under Arbitrary Magnetic Field

Ab Initio Simulations of Condensed Matter under Arbitrary Magnetic Field

Alfredo Correa, Eunseok Lee, Giulia Galli and Wei Cai

Abstract

---

[sections will be erased in PRL format, they are here to organize the editing]

Introduction

A plethora of ab initio electronic structure methods has been developed in the latest years. These detailed of implementation of these methods vary widely depending on the system of interest. For condensed systems, density function theory (DFT) within the planewaves (PW) formulation is a widely used approach due to considerations of efficiency and accuracy.

Within the PW method, periodic boundary conditions (PBC) or Bloch-periodic boundary conditions can be naturally applied, by adding k-points an arbitrarily accurate solution of the Kohn-Sham equations can be obtained. The method relies heavily in the so called Bloch theorem for periodic systems. As long as a periodic supercell can be defined (either exactly or as an approximation) for the condensed system, the method can be applied.

Unfortunately, in condensed systems where the Bloch theorem does not hold, the applicability is reduced. To give a realistic example, under external electric or magnetic field.

In the case of electric field the only non-perturbative approach requieres the use of the Berry phase [Resta, etc].

In the case of the magnetic field, these system had been treated only perturbatively or by approximate methods that are not compatible with condensed (extended) systems [Thonhauser]. Alternatively, a non-perturbative but limited solution requires the used of so called magnetic periodic boundary conditions (MPBC) [Cai04, Lee] and it is only valid for integer magnetic flux through the supercell.

The case of magnetic field is particularly interesting because even in the presence of uniform magnetic field in a periodic system the translational symmetry is compromised due to the quantum mechanical formulation. In a uniform magnetic field, all cells are still equivalent; however we can not apply Bloch theorem and straight PW methods fail because it can not impose the right boundary conditions and can not capture the group representation of the existing symmetries.

In this Letter we analyze the problem in a group theoretical framework for the single particle wavefunction that these ab initio methods rely on. This will inspire a complete solution that takes full advantage of the general symmetries. We propose an alternative method to describe the wavefunction that recover the irreducible representations that takes the place of Bloch theorem in presence of a magnetic field. We describe the method and propose certain algorithm for solving the one particle problem in presence of arbitrary uniform magnetic field. We implement the resulting method in a major DFT/PW code, and show possible applications for this new capability.

Problem statement

In the mentioned MPBC method it is required that the flux is an integer multiple of the magnetic flux quantum

$\Phi =n\Phi _{\circ }$

Such imposition is rather artificial and undermines the efficiency of the method since for small or moderated magnetic fields large supercell has to be imposed. The method proposed here does not suffer from the integer-flux limitation of the MPBC method.

The solution developed in Letter enjoys the possibility of imposing fractional magnetic flux conditions, which makes the applied magnetic field arbitrary in the simulation, such that

$\Phi ={\frac {p}{q}}\Phi _{\circ }$

Any magnetic field can be obtained by choosing an appropriate rational approximation. These two integers (numerator and denominator) determine the nature of the mathematical problem.

Background

We will briefly review the case with no magentic field and with magnetic field in order to stress the differences and similarities between the two.

No magnetic field

In the case of periodic potential

${\hat {H}}={\frac {{\hat {p}}^{2}}{2m}}+{\hat {V}}_{\text{periodic}}$

Bloch theorem holds and can be used to classify the solutions by its k-vector. The reason is two-fold. First, there is a set of quantum operators Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \hat T_\mathbf{N} = e^{i \hat\mathbf{p}\cdot \mathbf{N}} } that commute with the Hamiltonian operator $[{\hat {T}}_{\mathbf {N} },{\hat {H}}]=0,$ where $\mathbf {N}$ is an arbitrary lattice vector. Second, these quantum operators that are a representation of the translation group, commute with each other ${\hat {T}}_{\mathbf {N} _{1}}{\hat {T}}_{\mathbf {N} _{2}}={\hat {T}}_{\mathbf {N} _{2}}{\hat {T}}_{\mathbf {N} _{1}}$; this waranties that all irreducible representations are one dimensional

${\hat {T}}_{\mathbf {N} }\phi _{\mathbf {k} }=e^{-i\mathbf {k} \cdot \mathbf {N} }\phi _{\mathbf {k} }$

which is at the core of the applicability of the Bloch theorem. In particular this implies the following boundary conditions in real space for a cubic system:

$\phi _{\mathbf {k} }(x-a,y)=e^{-ik_{x}a}\phi _{\mathbf {k} }(x,y)$

$\phi _{\mathbf {k} }(x,y-a)=e^{-ik_{y}a}\phi _{\mathbf {k} }(x,y)$

In presence of magnetic field

In the presence of a magnetic field, the Hamiltonian changes to

${\hat {\mathcal {H}}}={\frac {({\hat {\mathbf {p} }}-e\mathbf {A({\hat {r}})} )^{2}}{2m}}+{\hat {V}}_{\text{periodic}}$

and the situation changes dramatically. The translation operators as they where defined do not commute with the magnetic Hamiltonian. In fact the transformation applied on the Hamiltonian operator

${\hat {T}}(\mathbf {N} ){\hat {\mathcal {H}}}{\hat {T}}^{\dagger }(\mathbf {N} )={\frac {({\hat {\mathbf {p} }}-e\mathbf {A({\hat {r}}-N)} )^{2}}{2m}}+{\hat {V}}_{\text{periodic}}$

induces a gauge transformation characterized by the gauge function $\xi$:

$\nabla \xi (\mathbf {x} )=\mathbf {A(r)-A(r-N)}$

A new set of magnetic translation has to be defined in order to compensate for a gauge transformation that the normal translation operators generate

${\hat {\mathcal {T}}}_{\mathbf {N} }=e^{ie\xi (\mathbf {N} )/\hbar }{\hat {T}}_{\mathbf {N} }$

These operators depend on the magnetic field and on the original gauge $\mathbf {A}$. They form a projective representation [Hammermesh] of the translation group with Schur factor system $\omega (\mathbf {N} _{1},\mathbf {N} _{2})=e^{-ie\xi _{\mathbf {N} _{2}}(\mathbf {N} _{1})/\hbar }$. They commute with the Hamiltonian $[{\hat {\mathcal {T}}}_{\mathbf {N} },{\hat {\mathcal {H}}}]\neq 0$ -- However Bloch theorem can not be applied since these operators do not commute with each other

${\hat {\mathcal {T}}}_{\mathbf {N} _{1}}{\hat {\mathcal {T}}}_{\mathbf {N} _{2}}=e^{ie(\mathbf {N} _{1}\times \mathbf {N} _{2})\cdot \mathbf {B} /\hbar }{\hat {\mathcal {T}}}_{\mathbf {N} _{2}}{\hat {\mathcal {T}}}_{\mathbf {N} _{1}},$

where we used the gauge-independent property $\xi _{\mathbf {N} _{1}}(\mathbf {N} _{2})-\xi _{\mathbf {N} _{2}}(\mathbf {N} _{1})=[\mathbf {N} _{1}\times \mathbf {N} _{2}]\cdot \mathbf {B}$.

The operators do not commute except in the precise case where the flux is an integer multiple of the quantum of flux. This exceptional case, which will be called integer magnetic flux or integer case, is exactly the case treated in the MPBC method, the particularity is that Bloch theorem can be again applied (but this time to the magentic translations for integer flux)

${\hat {\mathcal {T}}}_{\mathbf {N} }\phi _{\mathbf {k} }=e^{-i\mathbf {k} \cdot \mathbf {N} }\phi _{\mathbf {k} }\quad {\text{(integer case)}}$

In the case of $\mathbf {k} =\mathbf {0}$, this implies MPBC [Cai], for $\mathbf {k} \neq \mathbf {0}$ it implies Bloch-MPBC. In its simplest version, two dimensions and square cell and for the Landau gauge, equation (XX) implies

$\left.{\begin{array}{l}\phi _{\mathbf {k} }(x-a,y)=e^{-ik_{x}a}\phi _{\mathbf {k} }(x,y)e^{-ieBay/\hbar }\\\phi _{\mathbf {k} }(x,y-a)=e^{-ik_{y}a}\phi _{\mathbf {k} }(x,y)\end{array}}\right\}\quad {\text{(integer case)}}$

These equations boundary conditions make sense only when $Ba^{2}=n\Phi _{0}$ for the wavefunction to be singlevalued in the cell. (In this case $eBay/\hbar =n2\pi y/a$).

Solution

A generalization proposed here for the fractional case is obtained inspired in the fact that, since Bloch theorem is not valid we could in principle mix different k-point to make the boundary conditions consistent with a fractional flux. For example in the case of $\Phi =\Phi _{0}/2$ we could impose a different set of boundary conditions that involves a particular pair of k-points:

$\left.{\begin{array}{l}\phi _{(k_{x},k_{y})}(x-a,y)=e^{-ik_{x}a}\phi _{(k_{x},k_{y}+\pi /a)}(x,y)e^{-ieBay/\hbar }\\\phi _{(k_{x},k_{y})}(x,y-a)=e^{-ik_{y}a}\phi _{(k_{x},k_{y})}(x,y)\end{array}}\right\}\quad {\text{(half integer case)}}$

This shows that only half of the k-points are independent now. This makes sense since the Bloch theorem (and the k-ppoint classification) is valid in the doubled cell. The functions are singlevalued provided that $Ba^{2}=p\Phi _{0}/2$ (half-integer condition, i.e. $q=2$) and that we identify the points $k_{x}+2\pi /a=k_{x}$. (Note that in this case $eBay/\hbar =n\pi y/a$.)

To cope with this ambiguity we can change the notation and work with only half (in the y-direction) of the BZ and introduce another index. Instead of working with two k-points we can introduce a new index to the wavefunction.

${\begin{array}{l}\phi _{\mathbf {\kappa } }^{(0)}(x-a,y)=+e^{-i\kappa _{x}a}\phi _{\kappa }^{(1)}(x,y)e^{-ieBay/\hbar }\\\phi _{\mathbf {\kappa } }^{(1)}(x-a,y)=+e^{-i\kappa _{x}a}\phi _{\kappa }^{(0)}(x,y)e^{-ieBay/\hbar }\\\phi _{\mathbf {\kappa } }^{(0)}(x,y-a)=+e^{-i\kappa _{y}a}\phi _{\kappa }^{(0)}(x,y)\\\phi _{\mathbf {\kappa } }^{(1)}(x,y-a)=-e^{-i\kappa _{y}a}\phi _{\kappa }^{(1)}(x,y)\end{array}}$

In matrix form we can see a pattern emerge

${\begin{array}{l}\left[{\begin{array}{c}\phi _{\mathbf {\kappa } }^{(0)}\\\phi _{\mathbf {\kappa } }^{(1)}\end{array}}\right](x-a,y)=e^{-i\kappa _{x}a}\left[{\begin{array}{cc}0&+1\\+1&0\end{array}}\right]\left[{\begin{array}{c}\phi _{\kappa }^{(0)}\\\phi _{\kappa }^{(1)}\end{array}}\right](x,y)e^{-ieBay/\hbar }\end{array}}$

${\begin{array}{l}\left[{\begin{array}{c}\phi _{\mathbf {\kappa } }^{(0)}\\\phi _{\mathbf {\kappa } }^{(1)}\end{array}}\right](x,y-a)=e^{-i\kappa _{y}a}\left[{\begin{array}{cc}+1&0\\0&-1\end{array}}\right]\left[{\begin{array}{c}\phi _{\kappa }^{(0)}\\\phi _{\kappa }^{(1)}\end{array}}\right](x,y)\end{array}}$

We can reconstruct the magnetic translation operation in this particular representation:

${\begin{array}{l}{\hat {\mathcal {T}}}(a,0)\left[{\begin{array}{c}\phi _{\mathbf {\kappa } }^{(0)}\\\phi _{\mathbf {\kappa } }^{(1)}\end{array}}\right]=e^{-i\kappa _{x}a}\left[{\begin{array}{cc}0&+1\\+1&0\end{array}}\right]\left[{\begin{array}{c}\phi _{\kappa }^{(0)}\\\phi _{\kappa }^{(1)}\end{array}}\right]\end{array}}$

${\begin{array}{l}{\hat {\mathcal {T}}}(0,a)\left[{\begin{array}{c}\phi _{\mathbf {\kappa } }^{(0)}\\\phi _{\mathbf {\kappa } }^{(1)}\end{array}}\right]=e^{-i\kappa _{y}a}\left[{\begin{array}{cc}+1&0\\0&-1\end{array}}\right]\left[{\begin{array}{c}\phi _{\kappa }^{(0)}\\\phi _{\kappa }^{(1)}\end{array}}\right]\end{array}}$

The representation is obviously irreducible.

We are at this point ready to generalize to the case of higher denominators. For a flux defined by $p/q$, $q$ k-points have to be entangled via boundary conditions in the same way as before. In matrix form we obtain

[in the LaTeX file the notation will capture that phi is a vector (a kind of spinor)]

${\hat {\mathcal {T}}}(a,0)\phi _{\kappa }=e^{-i\kappa _{x}a}\Sigma _{1}^{\dagger }\phi _{\kappa }$

${\hat {\mathcal {T}}}(0,a)\phi _{\kappa }=e^{-i\kappa _{y}a}\Sigma _{3}^{\dagger }\phi _{\kappa }$

Where the matrices are generalizations of the matrices in Eq.(). The first being a (cyclic) shift matrix

$\Sigma _{1}=\left[{\begin{array}{cccc}0&1&&\\&0&1&\\&&\ddots &\ddots \\1&&&\end{array}}\right]$

And the second, a diagonal matrix with the qth roots of unity

$\Sigma _{3}=\left[{\begin{array}{cccc}1&&&\\&e^{i2\pi /q}&&\\&&\ddots &\\&&&e^{i2\pi (q-1)/q}\end{array}}\right]$

This is as close as we can get to replace the Bloch theorem for the case of fractional magnetic field -- a sort of fractional-magnetic-Bloch theorem. Note that the dimension of the representation is given by the denominator q of the magnetic flux in units of the quantum of flux $\Phi =p/q\Phi _{0}$ and that the continuous index $\kappa$ lives in a fraction (the qth part) of the original BZ.

A great simplification has been obtained by using the Landau gauge and a cubic cell. The same formulation can be obtained for non cubic cells by using a crystalline-Landau gauge.

Algorithm

As in the case of no-magnetic field Bloch-type envelope functions defined in the unit cell can be defined in the same way

$u_{\kappa }^{(n)}(\mathbf {r} )=e^{-i\kappa \cdot \mathbf {r} }\phi _{\kappa }^{(n)}(\mathbf {r} )\quad (n=0,\dots ,q-1)$

By replacing this definitions in Eq. (), () or () we obtain the boundary conditions for $u$. These boundary conditions are difficult to visualize in real space, however in an intermediate reciprocal space (Fourier transform applied in only one direction) the interpretation is straight forward.

The functions are normalized such that $\int _{\Omega }\sum _{i=0}^{q-1}|u_{\kappa }^{(i)}(\mathbf {r} )|^{2}d\mathbf {r} =1$.

Figure 1 illustrates the interpretation of the boundary conditions by connecting the two transformed functions at the ends of the cell for the $\Phi =1/2\Phi _{0}$ case.

Figure 1. (color online)