Transport routines within the kubo approach

[apezol]

Dear developers, I'm sharing some routines to calculate transport properties using SIESTA and sisl. I'm not sure if this submit is done the properly nor I'm an actual developer, so I've tried to comment the lines within the codes to highlight how they work.

This routines allow to calculate some transport quantities using an approach base on the Kubo formalism. The velocities are constructed using both derivatives of the Hamiltonian and Overlap matrices. Within this formalism, there are three examples showing how to calculate the spin Hall conductivity in Pt, the Orbital Hall conductivity in MoS2 and the Orbital Rashba Edelstein effect in Cu/O.

Example cases SHC in bulk Pt (reference https://arxiv.org/abs/1810.07637) - routine -> shc.py OHC in 2d MoS2 - routine -> ohc.py Orbital Edelstein effect in oxidized Cu (Cu/O) (used in https://pubs.aip.org/aip/apm/article/12/5/051105/3287821/Quantifying-the-large-contribution-from-orbital) - routine -> job_surface_orbital ( In this case the quantity is calculated as a projection on the layers forming the system)

Comments: routines_sisl_kubo.zip

*The main advantage of using sisl as a post-processing tool to siesta is that we don't need to rely in any wannierization procedure since the orbitals are already localized.

• The SHC and OHC are the intrinsic contributions to the Hall response and there's an issue appearing when Spin is set to SOC+onsite. I still don't know if this is related to the version of Siesta under which the SCF simulation was performed.

• Finally, if the routines aren't well explained through the comments, I'll be more than happy to help with any doubt.

• [x] Author : Armando Pezo

• [x] Which version of sisl it was originally developed and tested with sisl version: 0.0.0.dev+eacb8351ff416ad9f7f93732777c3d7931c5f711 routines_sisl_kubo.zip

[zerothi]

Thanks a lot for sharing your code! I will play a bit with it and see what can be done! :) And I'll probably ask you some questions! :)

[zerothi]

@apezol I am playing with your scripts to get them into the repo... After #790 I can now run your [35]*3 (Pt) BZ calculation on 6 cores (my laptop) in roughly 2 minutes. I'll share them there in a couple of days hopefully. :-)

[apezol]

@zerothi That's perfect, thanks a lot!. Such an enhancement will be very beneficial for studying other spin and/or optical properties that might be integrated over the BZ =D .

[zerothi]

@apezol could you perhaps clarify a bit more details on the L operator in the ohc.py code?

I am no expert here, so dumbing your discussion down (to my level) would be extremely helpful!! ;)

[apezol]

Sure, these operators come after the orbital angular momentum which can be defined in a given basis. In particular, they are defined like this

for the "p" states, and

for "d" states. Taken from https://arxiv.org/pdf/1808.05546.

What I did is, perform a 'rotation' such the basis is given in the order for which SIESTA takes the orbitals, something that we can check from the *.ORB_INDX file for instance. When I have SOC, I only need to perform a tensor product $L \otimes I_2$. The meanvalue of this operator is usually zero due to quenchin, however, when the system is perturbed (e.g. by an electric field), there should be finite orbital current generated. Hopefully this was clear enough.

[zerothi]

Thanks for these explanations.

One thing I am still not understanding is your double fermi-distributions?

In Eq. 4 of https://arxiv.org/pdf/1808.05546 it is written with a single Fermi distribution. In Eq. 1 of https://arxiv.org/abs/1810.07637 it is also only written with a single Fermi distribution. And here the energy dependence $\hbar\omega$ is implicit in the Lorentzian, but that is not how you have done it?

So why does this enter?

[zerothi]

Is it perhaps due to the discussion in https://journals.aps.org/prb/pdf/10.1103/PhysRevB.74.195118 surrounding Eq. 32?

[apezol]

Indeed, the formula containing the Berry (spin,orbital,etc) curvature is usually taken per band and it has the index (n) which ends up in the sum weighted by the FD distribution. I think that also can be done (in fact, it's the way it worked for the AHC in sisl, I guess). However in the routines I shared, it takes the form:

Taken from reference (https://arxiv.org/pdf/2105.14662). So, in that case the sum is running along both n and m indices. Other wise, it should be given the same if the (spin,orbital,etc) Berry curvature is calculated separately and then integrated. This would have the advantage that it can be use to plot the Berry curvature for every momentum $\vec{k}$ (maybe?)

[zerothi]

Sure, these operators come after the orbital angular momentum which can be defined in a given basis. In particular, they are defined like this

for the "p" states, and

I have tried to search around for these, why do they look like this? Is there a way to generalize this?

I.e. I only found this https://quantummechanics.ucsd.edu/ph130a/130_notes/node247.html which doesn't show the same matrices... Probably I am not understanding something ;)

[apezol]

Hi @zerothi, I think this is related to the fact that in a periodic environment like crystals, orbital angular momentum isn't conserved and then you need to consider a 'cubic harmonic' basis:

https://iopscience.iop.org/article/10.1088/1361-648X/ac5867/pdf (arxiv https://arxiv.org/pdf/2202.12475)

From eq. 23, you can see that they are obtained from the 'p' states such that the only non-vanishing matrix elements obey .

Indeed there are some works considering the usual expression of the atomic angular momentum (as the link you shared), but explicitly saying that they consider a spherical approximation:

From reference https://arxiv.org/pdf/cond-mat/0502345.

[zerothi]

So can this be generalized here? I.e. it depends on the basis? I guess not, the best would be to create the recipe for the calculation then?

[apezol]

Well, in SIESTA the order of the orbitals is always the same, I mean, from the file .ORB_INDX, it isalways like: s,px,pz,py,dxy.... and so on, and always on groups {s},{p's},{d's}. What I did is create the arrays L's for both p and d orbitals and from that, we can obtain a block diagonal operator containing these L operators in the order that they appear in .ORB_INDX. This will give a 'Total' angular momentum operator.

What might be tricky is how these need be organized to have the final diagonal form, once this will depend on the chosen basis for the simulations, that is: SP,DZP and so on, as well as if you introduce semi-core states, but isn't certainly hard to do, I guess.

[zerothi]

Well, in SIESTA the order of the orbitals is always the same, I mean, from the file .ORB_INDX, it isalways like: s,px,pz,py,dxy.... and so on, and always on groups {s},{p's},{d's}. What I did is create the arrays L's for both p and d orbitals and from that, we can obtain a block diagonal operator containing these L operators in the order that they appear in .ORB_INDX. This will give a 'Total' angular momentum operator.

What might be tricky is how these need be organized to have the final diagonal form, once this will depend on the chosen basis for the simulations, that is: SP,DZP and so on, as well as if you introduce semi-core states, but isn't certainly hard to do, I guess.

No, I mean, others might introduce a more generic approach to the basis used in sisl, currently it can be anything.

Currently, sisl can quite easily extract the l quantum numbers without problems. So figuring out which orbitals belongs to where shouldn't be a problem, what I am questioning is whether the L_p matrices are always the same, in case a user changes the Siesta NAO to something else?

[apezol]

Ah, ok I see. Yes, as long as SIESTA works with a LCAO basis, this should be the only way to write the operator for p and d states. Extending the same for 'f' states shouldn't be difficult if by any chance SIESTA is used for rare earths for instance.

So, going a little deeper, from my understanding, Wannier90 can give you a basis in terms of hybridized orbitals like sp for instance. SIESTA always has the same basis in terms of atomic orbitals (is that correct?), it's by this argument that I was saying that the operator should be always written the same way, that is, the L_p doesn't change.

[zerothi]

Ah, ok I see. Yes, as long as SIESTA works with a LCAO basis, this should be the only way to write the operator for p and d states. Extending the same for 'f' states shouldn't be difficult if by any chance SIESTA is used for rare earths for instance.

So, going a little deeper, from my understanding, Wannier90 can give you a basis in terms of hybridized orbitals like sp for instance. SIESTA always has the same basis in terms of atomic orbitals (is that correct?), it's by this argument that I was saying that the operator should be always written the same way, that is, the L_p doesn't change.

Ok, I see. So that means we can't implement it, since sisl can also read wannier90 basis information... :)

And yes, by default SIESTA always uses the same basis (although users can supply their own basis, if they want!).

[apezol]

Ah, ok I see. Yes, as long as SIESTA works with a LCAO basis, this should be the only way to write the operator for p and d states. Extending the same for 'f' states shouldn't be difficult if by any chance SIESTA is used for rare earths for instance. So, going a little deeper, from my understanding, Wannier90 can give you a basis in terms of hybridized orbitals like sp for instance. SIESTA always has the same basis in terms of atomic orbitals (is that correct?), it's by this argument that I was saying that the operator should be always written the same way, that is, the L_p doesn't change.

Ok, I see. So that means we can't implement it, since sisl can also read wannier90 basis information... :)

Ahhh, that's true!. In the case of Wannier90, this is something that only depends on the user's choice, sure.