1. Spin-orbit coupling

In this tutorial we will perform calculations related to physics in which spin-orbit coupling (SOC) plays an important role. In Fleur SOC typically is treated in so-called second variation. This means that a conventional Hamiltonian matrix, i.e., a Hamiltonian matrix with a scalar-relativistic description in the MT spheres, is set up and diagonalized. After diagonalization the eigenfunctions are used as a new basis for which a new Hamiltonian matrix with SOC is set up and diagonalized again to obtain the final energy eigenvalues and eigenfunctions.

This calculation procedure implies that the number of eigenfunctions that is calculated in the first step now becomes a cutoff parameter defining the basis set size for the second step. This parameter can be set in /calculationSetup/cutoffs/@numbands. A value of zero indicates that a material-dependent default value is chosen by the program. In principle SOC calculations have to be converged with respect to the numbands parameter. In this tutorial we will not do this but you are free to perform convergence tests to get a feeling for the relevance of this parameter. Since some of the investigated quantities are sensitive to very small energy differences a significant quantitative dependence on the cutoff parameters is expected.

1.1. Magnetocrystalline anisotropy of hcp Co

In perfect infinite ferromagnetic crystals the preferred magnetization direction is determined by spin-orbit coupling. In this example we will calculate the energy difference between different magnetization directions. Note that this typically is a very small quantity. Also be aware that magnetocrystalline anisotropy is not the only source of magnetic anisotropy. For example, in finite systems the shape of the sample also contributes to the magnetic anisotropy.

As example system we choose hcp Co. We generate an inp.xml file with the following inpgen input:

hcp Cobalt

&lattice latsys='hdp' a0=1.8897269 a=2.49680 c=4.03081 /

2
27 1.0 1.0 1.0
27 -1.0 -1.0 -1.0
&factor 3.0 3.0 4.0 /
&soc 0.10 0.37 /

Note the last line &soc 0.10 0.37 /. With this line we define a spin quantization axis (SQA) in terms of the two angles $\theta=0.10$ and $\varphi=0.37$ (spherical coordinate system). These are not the angles we will use in the actual calculations. We choose these angles in the inpgen input because the SQA breaks the symmetry of the crystal. Here we want to have the freedom to choose any other SQA after generating the inp.xml file. So we need to define an axis that breaks as many symmetries as possible. Otherwise a later choice of another axis in the Fleur input file may be incompatible to the here detected symmetry operations.

In the inp.xml file the angles appear in /calculationSetup/soc/@theta and /calculationSetup/soc/@phi. The SOC calculation is activated with the XML attribute calculationSetup/soc/@l_soc. If there is a spin quantization axis defined in the inpgen input spin-orbit coupling is activated by default.

Perform SOC calculations for the following SQAs:

1. theta="0.0", phi="0.0" (z axis)
2. theta="Pi/2.0", phi="0.0" (x axis)
3. theta="Pi/2.0", phi="Pi/2.0" (y axis)
4. theta="Pi/4.0", phi="0.0"

Extract the converged total energy for each of the calculations. hcp Co is an uniaxial system, i.e., the magnetocrystalline anisotropy in the x-y plane is negligible. Verify that the total energy difference between calculations 2 and 3 is much smaller than the difference to calculation 1. What is the "easy axis" (energetically most preferable), what is the "hard axis"?

Depending on the unit cell there are different models to describe the magnetocrystalline anisotropy. For uniaxial systems the total energy in dependence of the magnetization direction is expressed by the equation

$$\begin{equation} E_\text{total} = K_0 + K_1 \text{sin}^2(\theta) + K_2 \text{sin}^4(\theta) + \dots \end{equation}$$

We approximate that this expansion in terms of powers of sines of $\theta$ ends after the $K_2$ term. Use the total energies from calculations 1, 2, and 4 to set up a system of linear equations and determine the coefficients $K_0$, $K_1$, and $K_2$.

We want to visualize the anisotropy. For this we ignore the $K_0$ term. Plot it. You can adapt the following gnuplot script for this. (Note that it does not have a file output so you have to start gnuplot with the -p option to keep the image on the screen.)

set parametric
set isosamples 50
set urange [0:2.0*pi]
set vrange [0:pi]

a=1.0 # plug in K_1 here
b=0.0 # plug in K_2 here

k(v) = sin(v)*sin(v)*a
l(v) = sin(v)*sin(v)*sin(v)*sin(v)*b

fx(v,u) = sin(v)*cos(u) * (k(v)+l(v))
fy(v,u) = sin(v)*sin(u) * (k(v)+l(v))
fz(v) = cos(v) * (k(v)+l(v))
splot fx(v,u),fy(v,u),fz(v)

The result of plotting the magnetic anisotropy should be similar to the plot below.

2. Exercises

2.1. Magnetocrystalline anisotropy of fcc Ni

Ni crystallizes in a ccp structure with a lattice constant of $352.4~\text{pm}$. Set up an inpgen input for the determination of the magnetocrystalline anisotropy and generate an inp.xml file. In cubic systems like this one the magnetocrystalline anisotropy is described as

$$\begin{equation} E_\text{total} = K_0 + K_1 (\alpha_1^2 \alpha_2^2 + \alpha_2^2 \alpha_3^2 + \alpha_3^2 \alpha_1^2) + K_2 (\alpha_1^2 \alpha_2^2 \alpha_3^2) + \dots \end{equation}$$

where the $\alpha_i$ are the direction cosines.

Choose good SQAs for the determination of the $K_i$ coefficients up to $K_2$, calculate the total energies, and determine the coefficients. Plot the magnetocrystalline anisotropy. If you want to use the gnuplot script above for this you have to perform some adaptions.

Hint: Multiple equivalent axes are not good. Such a choice will lead to a set of linearly dependent equations.

Results to be delivered: 1. The values of the $K_i$ up to $K_2$. 2. A visualization of the respective plot.

2.2. Rashba effect in BiTeCl

The Rashba effect, similarly to the Dresselhaus effect, is a band splitting due to spin-orbit coupling that appears in systems with broken inversion symmetry. While the Rashba effect appears in low dimensional systems (films, heterostructures) and uniaxial bulk systems, the Dresselhaus effect is observed in cubic bulk systems. In this exercise we calculate the Rashba splitting in bulk BiTeCl. Note that it also appears in other bismuth tellurohalides (BiTeX with X=Cl,Br,I).

BiTeCl features a hexagonal unit cell (use latsys='hP' for inpgen) with 6 atoms and the lattice parameters $a=4.2426$ Angstrom and $c=12.3970$ Angstrom. The symmetry of the unit cell is defined by space group 186 ($P{6_3}mc$). The Bi atoms are at Wyckoff 2b positions with the representative atom being at $(1/3,2/3,0.3606)$ (in internal coordinates). The Te atoms are also at Wyckoff 2b positions. Here the representative atom is at $(1/3,2/3,0.7189)$. The Cl atoms are located at Wyckoff 2a positions with the representative atom at $(0.0,0.0,0.0024)$.

Note: the other atom positions can easily be obtained at the Bilbao Crystallographic Server -> Space-group symmetry -> WYCKPOS.

Set up an inpgen input for this system in which you specify the spin quantization axis with $\theta=0.0$ and $\varphi=0.0$.

We compare calculations with and without spin-orbit coupling.

Set numbands to 150, deactivate the SOC, perform a self-consistency calculation and afterwards a band structure calculation with 200 points along the $\text{K} - \Gamma - \text{M}$ path.

Set numbands to 150, activate the SOC, perform a self-consistency calculation and afterwards a band structure calculation with 200 points along the $\text{K} - \Gamma - \text{M}$ path.

Note: K = (1/3,1/3,0) ; M = (0.0,1/2,0.0)

Compare the two band structures.

Results to be delivered: Plots of the two band structures.