Some questions

[C-ZWang]

Hi, Dr Dalum, I have some problems regarding the calculation procedure, could you help me out?

1. About data generation. cols = gen_bands(Helicene(N=7), xsymbol=:δz, bounds=(0.5, 1.0), nsamples=100); gen_near_bands!(cols, f=TheoryOfCISS.calc_data1, bounds=(-5, 0), nsamples=100);

Based on my understanding. is the function 'gen_bands' to create the related Hamiltonian and coupling for the molecule and the function 'gen_near_bands' to calculate Green's function and transport quantities?

2. Spin orbit coupling term. There are not spin orbit coupling terms in the function f=TheoryOfCISS.calc_data1? But there are in calc_data3 and cal_data6. What are the differences between these calc_data1-calc_data6 ? Where in the code the SOC is introduced?

3. Are all the codes for photoemission calculation? Can we calculate the magnetic lead case using the existing code?

[dalum]

Thanks for the questions, let me try to answer them:

1. Your understanding is correct. gen_bands essentially sets up all the parts that depend only on the x parameter. So, as you suggest, it calculates all the non-energy dependent quantities such as the Hamiltonian, couplings, etc. These are contained within the LiteSimulation struct, and stored in the metadata field of DataColumn structs. gen_near_bands! then calculates the energy-dependent quantities such as Green's functions and transport properties.

2. calc_data1 calculates the spin and particle currents to first order in the SOC, with spin degrees traced out, as shown in the the paper. calc_data2 does the same, but gives the contribution from each individual carbon atom. calc_data3 calculates the same as calc_data1, but also stores both the full Green's function (with SOC) and the Green's function without, projected into the subspace of the two nearest eigenstates. I can't remember why I made calc_data4, it looks like it can be removed. calc_data5 is a poor man's way of calculating the time dependence of the polarization in a single lead, for an electron starting out on the molecule. This was for my PhD thesis and will be elaborated on in a future paper. calc_data6 is also related to that paper. Also notice that the spin-orbit coupling strength, λ, in the case of calc_data1, is not used until the spintransmission or polarization functions are called. I had a reason for doing it that way once, but it doesn't make a lot of sense anymore, so I will probably change that.

3. The code does not treat the magnetic lead case, as described in the Theory of CISS paper. In fact, we have recently discovered, by some free energy considerations, that the induced magnetization of the non-magnetic lead must be totally symmetric in the reversal of the magnetic field, which invalidates some of our arguments about the magnetic lead case in the paper. We are going to publish a correction to this.

[C-ZWang]

Thanks so much for your response, Dr Dalum. I now have a better understanding of your excellent work, and I have read your PhD thesis! It is so well written too! Hope, I can understand as much details as I can about your work :)

[dalum]

I am happy to help, and thank you for the kind words! :blush:

[C-ZWang]

Hi, Dr Dalum :) Sorry to bother you again... Do you know how to give a plot with the x-axis being the transport energy E and y-axis being the spin polarization, the twist angle is fixed to a certain value, that is equivalent to plot the two dimension view for your main figure with fixed twist angle theta (or psi in your paper). Thank you!!!

[dalum]

There are a few ways. To keep it similar to the heatmap plot, you can use gen_bands with nsamples = 1. To get a progress bar for the data generation, you then have to use the gen_col! function on cols[1]. If you don't care about that, you can still generate the data with gen_near_bands!, which accepts the same arguments as gen_col!, but the progress bar only updates every time a column has finished calculating.

So, to calculate the transport polarization along the z-axis for θ = pi/2, you would do:

cols = gen_bands(Polyacetylene(N=24, l_cc=1.4, l_ch=1.4), correct_overlaps = false, xsymbol = :θ, bounds=(pi/2, pi/2), nsamples=1);
gen_col!(cols[1], ysymbol = :E, bounds = (0, 3), nsamples=2000)
plot(cols[1].ys, map(percentage∘polarization(3), cols[1].data), xguide="E [eV]", yguide="polarization (%)")

[C-ZWang]

Got it! Another question is that we can change the chirality by changing the bounds (pi/2, pi/2) to (-pi/2,-pi/2)? Why is that if I take the bound as (0, 0) with chirality being 0, the spin polarization still having non zero values?

[dalum]

It took me a while to come to terms with it too, but it is related to the fact that the couplings between the leads and the molecule are, by default, randomly generated, positive definite matrices. Such a coupling is, in general, introducing chirality into the whole system of molecule + leads. Thus, you will in general get non-zero polarizations even for achiral molecules with such a coupling. You can change the seed that generates the coupling, by passing an integer to gen_bands using the keyword argument seed.

The effect should disappear if you set rng = 1.0I in gen_bands, which will force the couplings to be diagonal (though still only at the end atoms) and hence achiral. So,

cols = gen_bands(Polyacetylene(N=24, l_cc=1.4, l_ch=1.4), correct_overlaps = false, rng=1.0I, xsymbol = :θ, bounds=(pi/2, pi/2), nsamples=1);

When you play around with the couplings, you should find that the magnitude of the polarization, including its sign, will vary, but generally have its largest values when two states are close in the spectrum.

I just pushed some changes to both SlaterKoster.jl, Molecules.jl and this repository, that allows you to limit the couplings to some that are slightly more realistic, by calculating the coupling using a fictitious atom with manually specified Slater-Koster parameters. It is currently hard-coded in the LiteSimulation constructor in the lite.jl file. To enable it, you have to pass coupling_type = :atom. Additionally, with this update, if you want the diagonal coupling, in addition to specifying rng = 1.0I, you must also pass coupling_type = :identity.

Hope this helps!

[C-ZWang]

Thank you so much! Dr. Dalum. I have tried that, it seems the polarization is around 0 for E less than 2eV. But still each time, the curves are different. And I want to get the 0 regardless of energy range when twist angle is 0. Besides, I want to get opposite spin polarization for right and left handedness. Right now, I still cannot get what I want.. Maybe, I need to change the random number to some fixed number in the makecoupling function part? Another strange thing is that when I pass rng = 1.0I (use the command you give me above), it shows ERROR: LoadError: UndefVarError: I not defined. But in other code, they are OK. Do I need to use some package? I am so sorry to bother you so much. You can just point the way to solve the problems, because it is time consuming. I can read your thesis again to solve it with your advice.

[dalum]

Sorry, yes. You need to using LinearAlgebra for I to work.

[C-ZWang]

OK,Thank you!

[C-ZWang]

Hi, Dr Dalum, I want to ask two quick questions....

1. How do I pass the coupling_type= :identity into the code? I couldn't find where the variable is.

2. How could I get opposite spin polarization for opposite chirality? For example, if I set chirality to be left-handed, then at twist angle pi/2, I can get positive spin polarization; If I set chirality to be right-handed, then at twist angle -pi/2, I can get negative spin polarization.

Thank you!!!

[dalum]

1. In gen_bands (which relegates them to LiteSimulation). But it was added in the latest version of the code, a few days ago. Before that, you didn't have to set coupling_type.

2. Sorry I didn't answer this before. If you look at the implementation for Polyacetylene, you will find that negative twist angles, actually also ends up reversing the left and right lead, because the molecule now "grows" in the wrong direction, so you end up with the current running the other way as well, and so the polarization is the same. To fix it, you have to set handedness (defaults to 1) to -1, if you want to change the handedness: Polyacetylene(N=14, handedness=-1). Let me know if it doesn't change sign in that case, for non-chiral couplings, because it has every time I've tested it.

[C-ZWang]

Thank you so much, Dr Dalum! I deleted the old file and installed them, then I get perfect opposite spin polarization for different chirality. Last time, I updated the file, maybe it fails. Now, the spin polarization for theta=pi/2 is around 5%-30%. But, there is still spin polarization for theta=0, although generally, it is around 1%-15%. Do you have any idea about this? I expect it should be at least one order of magnitude smaller.

[dalum]

Is this with rng = 1.0I? If you are using calc_data1, it should be analytically zero for all energies in that case.

[C-ZWang]

This is what I use:

using TheoryOfCISS using LinearAlgebra cols = gen_bands(Polyacetylene(N=24,handedness=1, l_cc=1.4, l_ch=1.4), correct_overlaps = false, rng=1.0I, coupling_type = :identity, xsymbol = :θ, bounds=(0, 0), nsamples=1); gen_col!(cols[1], ysymbol = :E, bounds = (0, 3), nsamples=2000)

using Plots plot(cols[1].ys, map(percentage∘polarization(3), cols[1].data), xguide="E [eV]", yguide="polarization (%)")

And this is the figure I obtained:

[dalum]

Ah! That is expected. If you show a 3D model of polyacetylene at θ = 0 (do scatter(Polyacetylene(N=24, θ=0)()), and you should get a 3D scatter plot of the molecule), you will see that it is actually a highly unphysical model of a molecule. Essentially, it is a helix that has been squeezed flat, i.e., it has zero pitch. But we still model the leads as attaching to the ends of the molecule, so the electron will still have to do a circular motion to go from one end to the other. In other words, the presence of the leads break mirror symmetry of the system in the planes orthogonal to the z-axis, and so symmetry constraints only demand that the polarization along the x- and y-axes are zero. If you plot them, you should find that they indeed vanish.

To get actual achiral polyacetylene, you need to set θ = π. There's an oddity with range in Julia not accepting irrational values, so you need to convert it to a float first. So try:

cols = gen_bands(Polyacetylene(N=24,handedness=1, l_cc=1.4, l_ch=1.4), correct_overlaps = false, rng=1.0I, coupling_type = :identity, xsymbol = :θ, bounds=(float(π), float(π)), nsamples=1);

instead, and you should find that it has zero spin polarization along the z-axis. For this molecule, symmetry considerations now allow polarization along the y-axis though. But this is generally not measurable by experiments, and if the molecule is allowed to rotate around the z-axis, the "y-axis" is not well-defined.

So, I think the issue is that the name "chiral" is misleading in the name Chiral-Induced Spin Selectivity. If you consider the symmetry properties of an axial vector, chiral symmetry leaves it invariant, so you can't derive any properties from this. The relevant symmetries are mirror symmetries, and not purely chiral symmetry or lack thereof.

[C-ZWang]

Yes! Problem Solved! Still need time for me to fully understand the meaning of chirality! Thank you so much for you quick and thorough response, Dr Dalum!

[C-ZWang]

Hi, Dr Dalum. Could I ask a small question that is why the values of transmission T sometimes are larger than 1? Usually T is less or equal to 1.

[dalum]

Can you show an example?

[C-ZWang]

using TheoryOfCISS using LinearAlgebra

cols = gen_bands(Polyacetylene(N=24, handedness=1, l_cc=1.4, l_ch=1.4), correct_overlaps = false, rng=1.0I, coupling_type = :identity, xsymbol = :θ, bounds=(0, 0), nsamples=1);

cols = gen_bands(Polyacetylene(N=24,handedness=1, l_cc=1.4, l_ch=1.4), correct_overlaps = false, rng=1.0I, coupling_type = :identity, xsymbol = :θ, bounds=(float(π0.25), float(π0.25)), nsamples=1);

gen_near_bands!(cols, f=TheoryOfCISS.calc_data1, bounds=(0, 3), nsamples=2000);

gen_col!(cols[1], ysymbol = :E, bounds = (0, 3), nsamples=2000)

ns = 2000;

col_1 = cols[1].data; sx1 = []; sx2 = []; sx3 = []; t = []; for (idx, col) in enumerate(col_1) push!(sx1,real(col_1[idx][4][1])); push!(sx2,real(col_1[idx][4][2])); push!(sx3,real(col_1[idx][4][3])); push!(t,abs(col_1[idx][5])); end

L0 = (1:1:2000)/20003; sp1 = 6e-1sx1./t; sp2 = 6e-1sx2./t; sp3 = 6e-1sx3./t;

using Plots

plot(L0, t, xguide="E [eV]", yguide="Spin Polarization") # Spin Polarization Spin component Transmission

plot(L0, sx3, xguide="E [eV]", yguide="Spin component") # Spin Polarization Spin component Transmission

plot(L0, t, xguide="E [eV]", yguide="Transmission")

[C-ZWang]

Oh no! The code is roughly like this, but multiply symbol does not show in the copy process. π0.25, L0 = (1:1:2000)/2000*3

[dalum]

Uh-huh. I think it is just a normalisation issue, but I thought I had fixed all of those. I'll take a look at it later. It shouldn't affect the polarisation, though, since the normalisation of the transmission and spin transmission cancel.

[dalum]

Just had a look. If you look in the calc_data1 function in lite_datagen.jl, you will see that the transmission, t, is defined as:

X = G*sim.γL*G'
t = 2real(dot(X, sim.ΓR))

where dot(X, sim.ΓR) is just a slightly faster way to do tr(X'sim.ΓR). dot here is the Frobenius inner product (https://en.wikipedia.org/wiki/Frobenius_inner_product). The factor of 2 comes from the derivation as the trace over spin, so the transmission is normalised to two instead of unity.

The normalisation may be a bit worse for the randomly generated matrices. It would probably be nice to fix it.

[C-ZWang]

Gotcha! Thank you so much!!!

[C-ZWang]

Dear Dr. Dalum, thanks so much for your help! Can I get a connection with you via Email? My Email address is: cwang239@asu.edu