# Kwant KPM Module
Here we will learn how to use the KPM module in Kwant.
The outline in this section is:
<div class="lower_greek">
0. The math behind it
1. KPM expansion of a function
1. Spectral densities
1. Density of States
2. Spectral density of a general operator
2. Add moments
1. Choosing energy resolution or number of moments?
3. Integrate
1. What is the effect of temperature?
4. Vector factories
1. Random vectors and the stochastic trace
2. Local vectors
3. Choosing projections over sites and orbitals
5. Conductivity
1. Kubo-Bastin expansion of the conductivity tensor
2. Spin-Hall conductivity
6. Exercises
1. (easy) staggered potential: a simplified model of hBN is a honeycomb lattice with a large staggering potential.
* Compute the DoS for each sublattice (B and N).
* What can you say about the population of the valence and conduction band?
2. (hard)Spin-Hall conductivity: for typical systems[1] the spin-Hall conductivity is equal to the difference between spin-up and spin-down Hall conductivity.
* Construct a Kane-Mele system (spin-full sites with s_z times the Haldane hoppings).
* Get the difference between spin-up and spin-down conductivity.
* What do you think about the results?
* Is this quantization robust? What term can you add to the Hamiltonian to break it?
</div>
<!-- #region -->
## 0. The math behind it
The KPM expansion [\[1\]][1] is based on the expansion of a function that depends on the Hamiltonian in the basis of
Chebyshev polynomials.
This expansion is particularly stable, and very efficient, thanks to the recursion relation of Chebyshev polynomials. For details see Ref. [\[1\]][1].
\[1\] <https://arxiv.org/abs/cond-mat/0504627> "Rev. Mod. Phys., Vol. 78, No. 1 (2006)"
[1]: <https://arxiv.org/abs/cond-mat/0504627> "Rev. Mod. Phys., Vol. 78, No. 1 (2006)"
<!-- #endregion -->
<!-- #region -->
### 0.A. KPM expansion of a function
The Chebyshev polynomials are defined as
and they hold the recursion relation
\begin{align}
T_0(x)&=1\\
T_1(x)&=x\\
T_{n+1}(x) &= 2xT_n(x) - T_{n-1}(x).
\end{align}
These polynomials form a complete basis on the interval $x \in (-1,1)$, and any function $f(x)$ on that interval can be expanded as
$$ f(x) = \dfrac{1}{\pi\sqrt{1-x^2}}\left(\mu_0 + 2\sum_{n=1}^\infty {\mu_n T_n(x)} \right),
$$
where the coefficients $\mu_n$, or *moments* are
$$\mu_n = \int_{-1}^1f(x) T_n(x) \mathrm{d}x.
$$
<!-- #endregion -->
### Total spectral density
The KPM expansion can be applied to operators, such as the Hamiltonian $H$.
To obtain the density of states, we will expand the *delta function* of $\delta(E-H)$.
Let's start from a Hamiltonian $H$ with a complete set of eigenvectors $\lvert k\rangle$, such that
$$ H\lvert k\rangle=E_k \lvert k\rangle.$$
The density of states is defined by
$$\rho(E) = \sum_{k} \delta(E-E_k)=\mathrm{Tr}[\delta(E-H)].
$$
The expansion of the *delta functions* in terms of Chebyshev polynomials is determined by the moments
\begin{align}
\mu_n = &\, \int_{-1}^1 \rho(E) T_n(E) \mathrm{d}E \\
\mu_n = &\, \int_{-1}^1 \sum_{k} \delta(E-E_k) T_n(E) \mathrm{d}E \\
\mu_n = &\, \sum_{k} T_n(E_k)
\end{align}
Now we can use the fact that $T_n(H)\lvert k\rangle = T_n(E_k) \lvert k\rangle$
\begin{align}
\mu_n = &\, \sum_{k} \langle k\lvert T_n(H) \lvert k\rangle\\
\mu_n = &\, \mathrm{Tr}\left[T_n(H)\right]
\end{align}
Applying the recursion relation for Chebyshev polynomials, we can see that the reason behind the efficiency of the method is due to the fact that
\begin{eqnarray}
&\lvert \alpha_n \rangle =& T_n(H) \lvert k\rangle\\
\Rightarrow \,\,&
\lvert \alpha_{n+1}\rangle =& 2 H \lvert \alpha_n \rangle - \lvert \alpha_{n-1} \rangle,
\end{eqnarray}
where the moments are
$$\mu_n = \langle \alpha_0 \lvert \alpha_n\rangle,
$$
therefore, we need only one matrix-vector multiplication for every new moment in the expansion of the density of states.
## 1. Spectral densities
```python
import numpy as np
import kwant
from kwant import kpm
from kwant import HoppingKind
def lattice_and_neighbors(norbs=1, name=('a', 'b')):
honeycomb = kwant.lattice.honeycomb(1, norbs=norbs, name=name)
honeycomb_a, honeycomb_b = honeycomb.sublattices
# neighbors for honeycomb lattice
first_neighbors = [
HoppingKind((0, 1), honeycomb_a, honeycomb_b),
HoppingKind((0, 0), honeycomb_a, honeycomb_b),
HoppingKind((-1, 1), honeycomb_a, honeycomb_b)
]
second_neighbors = [
HoppingKind(( 0, 1), honeycomb_a, honeycomb_a),
HoppingKind(( 1,-1), honeycomb_a, honeycomb_a),
HoppingKind((-1, 0), honeycomb_a, honeycomb_a),
HoppingKind(( 1, 0), honeycomb_b, honeycomb_b),
HoppingKind(( 0,-1), honeycomb_b, honeycomb_b),
HoppingKind((-1, 1), honeycomb_b, honeycomb_b),
]
return honeycomb, first_neighbors, second_neighbors
```
```python
# construct a Hamiltonian of a finite system
def graphene(radius=20, norbs=1, t=-1, t2=0):
"""Construct by default a simple graphene model.
If 't2' is imaginary, different from zero, a Haldane model.
"""
honeycomb, first_neighbors, second_neighbors = lattice_and_neighbors(norbs=norbs)
syst = kwant.Builder()
def disk(pos):
return np.linalg.norm(pos) <= radius
syst[honeycomb.shape(disk, (0, 0))] = 0 * np.eye(norbs)
syst[first_neighbors] = t * np.eye(norbs)
syst.eradicate_dangling()
if t2:
syst[second_neighbors] = t2 * np.eye(norbs)
fsyst = syst.finalized()
return fsyst
```
```python
import holoviews as hv
hv.notebook_extension()
def plot_spectrum(spectrum):
"""Plot the energies vs. densities. If the spectrum is complex,
plot the imaginary part too."""
e, d = spectrum()
d = np.real_if_close(d)
the_plot = hv.Curve((e, d.real))
if type(d) is complex:
the_plot = the_plot * hv.Curve((e, d.imag))
return the_plot.redim(x='$e$', y='$\\rho$')
```
### 1.A. Density of States
```python
radius = 20
fsyst = graphene(radius=radius)
kwant.plot(fsyst);
```
```python
%%time
spectrum = kpm.SpectralDensity(fsyst)
plot_spectrum(spectrum)
```
### 1.B. Spectral density of a general operator
```python
```
## 2. Add moments
The KPM module uses the information of the SpectralDensity instance to
generate more orders in the expansion
### 2.A. Choosing energy resolution or number of moments?
The energy resolution is roughly the width of the expansion of a delta peak.
This resolution scales as the inverse number of moments:
`energy_resolution <= (e_max - e_min) / (pi * m)`
```python
%%time
print('num_moments before:', spectrum.num_moments)
spectrum.add_moments(300)
print('num_moments after:', spectrum.num_moments)
```
```python
plot_spectrum(spectrum)
```
## 3. Integrate
The integral of a KPM expansion has almost no computational cost, and provides
useful information on the system. For example, the number of filled states
below a certain Fermi energy.
$$n = \int_{-\infty}^{e_F} \rho(E) f(E) dE$$
Where $f(E)$ is the zero temperature Fermi distribution function.
```python
e_F = 0.
print(spectrum.integrate(distribution_function=lambda e: 1*(e<e_F)).real)
```
```python
# compare with half the number of sites in the system
print(fsyst.graph.num_nodes / 2)
```
### 3.A. What is the effect of temperature?
Before we assumed zero temperature, so the Fermi distribution function was a
step function with values `1` before `e_F` and `0` after `e_F`.
When there is a finite temperature, the distribution function is smooth and
slowly changes from `1` to `0` over a range of energies similar to the temperature.
We can get the number of filled states at finite temperature `temp` by calling
```python
temp = 0.1 # an arbitrary value as an example
fermi_array = kpm.fermi_distribution(spectrum.energies, mu=e_F, temperature=temp)
hv.Curve((spectrum.energies, fermi_array))
```
```python
def distribution_function(e):
return kpm.fermi_distribution(energy=e, mu=e_F, temperature=temp)
print(spectrum.integrate(distribution_function=distribution_function).real)
```
## 4. Vector factories
The KPM module provides two vector factories that largely influence the output
of the `SpectralDensity` function:
* RandomVectors
* LocalVectors
Besides that, you can define your own vector factory, which must be any python
iterable of vectors: a list, a tuple, a numpy ndarray, etc.
### 4.A. Random vectors and the stochastic trace
So far, we totally overlooked the fact that the spectrum looks a bit noisy.
This is due to the stochastic trace approximation that is used by default.
In this approximation, the trace is not computed explicitly over the Hilbert
space, but approximated with a number of random vectors:
$\mu_n = \mathrm{Tr}\left(T_n(H)\right) \approx \dfrac{1}{R} \sum_{r} \langle r\lvert T_n(H) \lvert r\rangle$.
The default vector factory is `kpm.random_vectors`.
You can improve the approximation by using more random vectors, and the average
densities will be the new average that includes those vectors.
```python
%%time
spectrum.add_vectors(10)
plot_spectrum(spectrum)
```
### 4.B. Local vectors
The `kpm.LocalVectors` is a vector factory that does not rely on the stochastic trace,
but samples the whole Hilbert space exactly: $\mu_n = \, \sum_{k} \langle k\lvert T_n(H) \lvert k\rangle$.
Therefore, it is in general more computationally demanding.
The number of total vectors that cover the Hilbert space of the system is `num_sites * norbs` (assuming all sites with the same number of orbitals):
```python
print(f'Total number of orbitals: {len(fsyst.sites) * fsyst.sites[0].family.norbs}')
```
Each of the vectors created by this factory is a vector on the basis
of your Hamiltonian, with `1` on one site/orbital, and zeroes everywhere else.
### 4.C. Choosing projections over sites and orbitals
Choosing RandomVectors or LocalVectors greatly depends on your system and on
what you want to compute. There is another degree of freedom to create either
of these vector factories, which is the `where` argument.
`where` can be a list of sites, or indices, or a function that returns `True`
for the sites that you want to select. On the basis of those sites, you can
generate random vectors or local vectors.
For example, we can choose to compute the DoS only at the bulk of our sample,
so we can define
```python
# average over a disk half the radius
def center_disk(site):
return np.linalg.norm(site.pos) < (radius * 0.5)
center_rand_vecs = kpm.RandomVectors(fsyst, where=center_disk)
one_rand_vec = next(center_rand_vecs) # gets one vector from the generator
kwant.plot(fsyst, site_color=np.abs(one_rand_vec), cmap='Blues');
```
```python
spectrum = kpm.SpectralDensity(fsyst, vector_factory=center_rand_vecs)
plot_spectrum(spectrum)
```
```python
# sample only the center unitcell, representative of pristine graphene
def center_uc(site):
return site.tag == (0,0) # A and B sites at tag (0, 0)
center_local_vecs = kpm.LocalVectors(fsyst, where=center_uc)
print(f'This factory has {len(center_local_vecs)} vectors.')
spectrum = kpm.SpectralDensity(fsyst, vector_factory=center_local_vecs,
num_vectors=None)
# passing a larger 'num_vectors' than available will raise an error
plot_spectrum(spectrum)
```
<!-- #region -->
## 5. Conductivity
### 5.A. Kubo-Bastin expansion of the conductivity tensor
The implementation in Kwant follows Ref. [\[2\]][2], where the Kubo-Bastin
conductivity tensor is approximated with two KPM expansions:
$$
\sigma_{\alpha \beta}(\mu, T) =
\frac{1}{V} \int_{-\infty}^{\infty}{\mathrm{d}E} f(\mu-E, T)
\left({j_\alpha \rho(E) j_\beta \frac{\mathrm{d}G^{+}}{\mathrm{d}E}} -
{j_\alpha \frac{\mathrm{d}G^{-}}{\mathrm{d}E} j_\beta \rho(E)}\right),
$$
where $\mu$ is the chemical potential, $T$ the temperature, and $j_{\alpha}$ $(j_{\beta})$ are
the current operators in direction $\alpha$ $(\beta)$.
The prefactor $1/V$ is the normalization with the volume (area in 2D) of the system.
\[2\] <https://arxiv.org/abs/1410.8140> "Real-space calculation of the conductivity tensor for disordered topological matter. Phys. Rev. Lett. 114, 116602 (2015)"
[2]: <https://arxiv.org/abs/1410.8140> "Real-space calculation of the conductivity tensor for disordered topological matter. Phys. Rev. Lett. 114, 116602 (2015)"
<!-- #endregion -->
```python
%%time
# first, let's make a topological system, a spinful Haldane model
radius = 20
norbs = 2 # spinful
t2 = 0.5j # SOI of the Haldane model
fsyst = graphene(radius=radius, norbs=norbs, t=-1, t2=t2)
center_local_vecs = kpm.LocalVectors(fsyst, where=center_uc)
# these vectors will cover 2 sites (A and B), and two orbitals on each
print(f'There are {len(center_local_vecs)} local vectors')
cond_xy = kpm.conductivity(
fsyst, alpha='x', beta='y',
vector_factory=center_local_vecs,
num_vectors=None,
mean=True, # return the average over the vectors
)
```
```python
# the conductivity inside the gap should be quantized
print('\nConductivity raw output:\n', cond_xy(0).real)
```
Here, we must remember about the normalization prefactor. In this case, the area should be the area covered by each of vectors.
Finally, since we output the `mean` over these vectors, and not the `sum`, we multiply by the number of vectors.
```python
# We must rescale with the unitcell area per site and the number of orbitals
honeycomb = kwant.lattice.honeycomb(norbs=norbs)
area_uc = np.cross(*honeycomb.prim_vecs)
area_uc_per_site = area_uc / 2
print('\nConductivity scaled with the area of the unit cell:\n',
cond_xy(0).real / area_uc_per_site * norbs)
```
```python
%%time
cond_xy_array = np.array([cond_xy(mu=e, temperature=0.) for e in cond_xy.energies])
# we must normalize with the unitcell area
cond_xy_array = cond_xy_array / area_uc_per_site * norbs
```
```python
hv.Curve((cond_xy.energies, cond_xy_array.real),
kdims='$e\,[t]$', vdims='$\sigma_{xy}\,[e^2/h]$')
```
<!-- #region -->
## 6. Exercises
### 6.A. (easy) Staggering potential: a simplified model of hBN is a honeycomb lattice with a large staggering potential.
* Add a staggering potential to the graphene system 'simple_graphene'
- Hint: modify `graphene()` to pass two onsite energies for the `A` and `B` sublattices with
```
honeycomb_a, honeycomb_b = honeycomb.sublattices
syst[honeycomb_a.shape(disk, (0, 0)) = e_A
syst[honeycomb_b.shape(disk, (0, 0)) = e_B
```
* Compute the DoS for each sublattice (`A` and `B` represent boron and nitrogen).
* What can you say about the population of the valence and conduction band?
### 6.B. (hard) Spin-Hall conductivity: for typical systems\* the spin-Hall conductivity is equal to half the difference between spin-up and spin-down Hall conductivity.
\* typical means, systems where the Hamiltonian commutes with the spin sigma_z operator, which means that the `z` spin projection of the electrons is conserved. See [\[3\]][3] for details.
* Modify the `graphene()` function to construct a Kane-Mele system.
- Hint: It will be similar to `graphene(norbs=2)`, with spin-full sites, but should have the Pauli matrix `s_z` in the Haldane second neighbors hopping.
* Get the (total) Hall conductivity. How is it different from the Haldane model?
* Get the difference (multiplied by `1/2`) between spin-up and spin-down conductivity. What do you think about the results?
* Hint: The local vector factory produces one vector per orbital (spin projection), in order. Then,
```
center_local_vec_factory = kpm.LocalVectors(fsyst, where=center_uc)
center_local_vecs = list(center_local_vec_factory) # get all the vectors from this factory
center_local_vecs_spin_down = center_local_vecs[0::2] # select every second vector
center_local_vecs_spin_up = center_local_vecs[1::2]
```
* Is this quantization robust? What term can you add to the Hamiltonian to break it? See [\[4\]][4].
* Think about what happens if the spin is conserved, but in a different direction than `z`.
Do you think it is still possible to obtain a quantized spin-Hall conductivity? See [\[5\]][5] for a material where this happens.
[4]: <https://arxiv.org/abs/cond-mat/0701293> " Theory of conserved spin current and its application to a two-dimensional hole gas. Phys. Rev. B 77, 075304 (2008)"
[3]: <https://arxiv.org/abs/2006.12071> "Spin Hall conductivity in insulators with non-conserved spin. Phys. Rev. B 102, 125138 (2020)"
[5]: <https://arxiv.org/abs/2007.05626> "Canted Spin Texture and Quantum Spin Hall Effect in WTe2"
<!-- #endregion -->
## References
\[1\] <https://arxiv.org/abs/cond-mat/0504627> "The Kernel Polynomial Method. Rev. Mod. Phys., Vol. 78, No. 1 (2006)"
[1]: <https://arxiv.org/abs/cond-mat/0504627> "The Kernel Polynomial Method. Rev. Mod. Phys., Vol. 78, No. 1 (2006)"
\[2\] <https://arxiv.org/abs/1410.8140> "Real-space calculation of the conductivity tensor for disordered topological matter. Phys. Rev. Lett. 114, 116602 (2015)"
[2]: <https://arxiv.org/abs/1410.8140> "Real-space calculation of the conductivity tensor for disordered topological matter. Phys. Rev. Lett. 114, 116602 (2015)"
\[3\] <https://arxiv.org/abs/cond-mat/0701293> "Theory of conserved spin current and its application to a two-dimensional hole gas. Phys. Rev. B 77, 075304 (2008)"
[3]: <https://arxiv.org/abs/cond-mat/0701293> "Theory of conserved spin current and its application to a two-dimensional hole gas. Phys. Rev. B 77, 075304 (2008)"
\[4\] <https://arxiv.org/abs/2006.12071> "Spin Hall conductivity in insulators with non-conserved spin. Phys. Rev. B 102, 125138 (2020)"
[4]: <https://arxiv.org/abs/2006.12071> "Spin Hall conductivity in insulators with non-conserved spin. Phys. Rev. B 102, 125138 (2020)"
\[5\] <https://arxiv.org/abs/2007.05626> "Canted Spin Texture and Quantum Spin Hall Effect in WTe2"
[5]: <https://arxiv.org/abs/2007.05626> "Canted Spin Texture and Quantum Spin Hall Effect in WTe2"
## 7. Solutions
### 7.A. (easy) Staggering potential: a simplified model of hBN is a honeycomb lattice with a large staggering potential.
* Add a staggering potential to the graphene system 'simple_graphene'
- Hint: modify `graphene()` to pass two onsite energies for the `A` and `B` sublattices with
```
honeycomb_a, honeycomb_b = honeycomb.sublattices
syst[honeycomb_a.shape(disk, (0, 0)) = e_A
syst[honeycomb_b.shape(disk, (0, 0)) = e_B
```
* Compute the DoS for each sublattice (`A` and `B` represent boron and nitrogen).
* What can you say about the population of the valence and conduction band?
```python
# construct a Hamiltonian of a finite system
def graphene_staggered(radius=20, norbs=1, staggering=0.5, t=-1, t2=0):
"""Construct by default a simple graphene model.
If 't2' is imaginary, different from zero, a Haldane model.
"""
honeycomb, first_neighbors, second_neighbors = lattice_and_neighbors(norbs=norbs)
honeycomb_a, honeycomb_b = honeycomb.sublattices
syst = kwant.Builder()
def disk(pos):
return np.linalg.norm(pos) <= radius
syst[honeycomb_a.shape(disk, (0, 0))] = staggering * np.eye(norbs)
syst[honeycomb_b.shape(disk, (0, 0))] = -staggering * np.eye(norbs)
syst[first_neighbors] = t * np.eye(norbs)
syst.eradicate_dangling()
if t2:
syst[second_neighbors] = t2 * np.eye(norbs)
fsyst = syst.finalized()
return fsyst
```
```python
fsyst = graphene_staggered()
# sample only the center unitcell, representative of pristine graphene
def center_uc_a(site):
return site.tag == (0,0) and site.family.name == 'a'
def center_uc_b(site):
return site.tag == (0,0) and site.family.name == 'b'
center_local_vecs_a = kpm.LocalVectors(fsyst, where=center_uc_a)
center_local_vecs_b = kpm.LocalVectors(fsyst, where=center_uc_b)
spectrum_a = kpm.SpectralDensity(fsyst, vector_factory=center_local_vecs_a,
num_vectors=None)
spectrum_b = kpm.SpectralDensity(fsyst, vector_factory=center_local_vecs_b,
num_vectors=None)
# Plot both spectrums
plot_spectrum(spectrum_a).relabel('A') * plot_spectrum(spectrum_b).relabel('B')
```
The population of the valence band near the gap is concentrated on the `B` sites, while the population of the conductance band is concentrated on the `A` sites.
### 7.B. (hard) Spin-Hall conductivity: for typical systems\* the spin-Hall conductivity is equal to the difference between spin-up and spin-down Hall conductivity.
* Modify `graphene()` to construct a Kane-Mele system
```python
# construct a Hamiltonian of a finite system
def graphene_Kane_Mele(radius=20, norbs=2, t=-1, t2=0.5j):
"""Construct by default a simple graphene model.
If 't2' is imaginary, different from zero, a Haldane model.
"""
s_z = np.array([[1, 0], [0, -1]])
honeycomb, first_neighbors, second_neighbors = lattice_and_neighbors(norbs=norbs)
syst = kwant.Builder()
def disk(pos):
return np.linalg.norm(pos) <= radius
syst[honeycomb.shape(disk, (0, 0))] = 0 * np.eye(norbs)
syst[first_neighbors] = t * np.eye(norbs)
syst.eradicate_dangling()
syst[second_neighbors] = t2 * s_z
fsyst = syst.finalized()
return fsyst
# first, let's make a topological system, a spinful Haldane model
radius = 20
norbs = 2 # spinful
t2 = 0.5j # SOI of the Haldane model
fsyst = graphene_Kane_Mele(radius=radius, norbs=norbs, t=-1, t2=t2)
```
* Get the total conductivity
```python
# sample only the center unitcell, representative of pristine graphene
def center_uc(site):
return site.tag == (0,0)
center_local_vecs = kpm.LocalVectors(fsyst, where=center_uc)
# these vectors will cover 2 sites (A and B), and two orbitals on each
cond_xy = kpm.conductivity(
fsyst, alpha='x', beta='y',
vector_factory=center_local_vecs,
num_vectors=None,
mean=True, # return the average over the vectors
)
# We must rescale with the unitcell area per site and the number of orbitals
honeycomb = kwant.lattice.honeycomb(norbs=norbs)
area_uc = np.cross(*honeycomb.prim_vecs)
area_uc_per_site = area_uc / 2
cond_xy_array = np.array([cond_xy(mu=e, temperature=0.) for e in cond_xy.energies])
# we must normalize with the unitcell area
cond_xy_array = cond_xy_array / area_uc_per_site * norbs
```
```python
hv.Curve((cond_xy.energies, cond_xy_array.real),
kdims='$e\,[t]$', vdims='$\sigma_{xy}\,[e^2/h]$')
```
We see zero everywhere!
* Get the difference between spin-up and spin-down conductivity
```python
# sample only the center unitcell, representative of pristine graphene
def center_uc(site):
return site.tag == (0,0)
center_local_vec_factory = kpm.LocalVectors(fsyst, where=center_uc)
center_local_vecs = list(center_local_vec_factory) # get all the vectors from this factory
center_local_vecs_spin_down = center_local_vecs[0::2] # select every second vector
center_local_vecs_spin_up = center_local_vecs[1::2]
cond_xy_up = kpm.conductivity(
fsyst, alpha='x', beta='y',
vector_factory=center_local_vecs_spin_up,
num_vectors=None,
mean=True, # return the average over the vectors
)
cond_xy_down = kpm.conductivity(
fsyst, alpha='x', beta='y',
vector_factory=center_local_vecs_spin_down,
num_vectors=None,
mean=True, # return the average over the vectors
)
# We must rescale with the unitcell area per site and the number of orbitals
honeycomb = kwant.lattice.honeycomb(norbs=norbs)
area_uc = np.cross(*honeycomb.prim_vecs)
area_uc_per_site = area_uc / 2
cond_xy_up_array = np.array([cond_xy_up(mu=e, temperature=0.) for e in cond_xy.energies])
cond_xy_down_array = np.array([cond_xy_down(mu=e, temperature=0.) for e in cond_xy.energies])
# we must normalize with the unitcell area
cond_xy_up_array = cond_xy_up_array / area_uc_per_site * norbs
cond_xy_down_array = cond_xy_down_array / area_uc_per_site * norbs
```
```python
(
hv.Curve((cond_xy_down.energies, cond_xy_down_array.real),
kdims='$e\,[t]$', vdims='$\sigma_{xy}\,[e^2/h]$', label='down') *
hv.Curve((cond_xy_up.energies, cond_xy_up_array.real),
kdims='$e\,[t]$', vdims='$\sigma_{xy}\,[e^2/h]$', label='up')
).relabel('Separated spin down and up')
```
```python
cond_xy_spin_hall = (cond_xy_down_array.real-cond_xy_up_array.real) / 2
hv.Curve((cond_xy_down.energies, cond_xy_spin_hall),
kdims='$e\,[t]$', vdims='$\sigma^z_{xy}\,[e^2/h]$', label='Spin-Hall conductivity')
```