Trying to replicate a figure

I am trying to replicate a figure from the paper (Dyrdał, A., and J. Barnaś. "Anomalous, spin, and valley Hall effects in graphene deposited on ferromagnetic substrates." 2D Materials 4.3 (2017): 034003). I have attached the required figure in terms of figure-1, figure-2 and figure-3. I want to replicate figure-6. I have attached the required equations and my code. I am not getting where I am getting the error.

% Defining the parameters   
lambdar= 0.010;        % Rashba parameter (10 meV)
hbar=6.582e-16;    %Reduced Plancks constant in meV-s
vf= 0.812e-6;       % Fermi velocity in m/s .
v=hbar*vf; 
e=1; 

mu1= linspace(-0.001, 0.001, 5);    %Fermi potential
lambdaex1= linspace(0, 0.020, 5);  %Fixed Exchange parameter of 20meV

[mu,lambdaex]= meshgrid(mu1, lambdaex1);

% Defining two notations
k3p=(1/v).*sqrt(lambdaex.^2+mu.^2+2.*(sqrt(mu.^2.*(lambdaex.^2+lambdar.^2))-(lambdaex.^2*lambdar.^2)));
k3n=(1/v).*sqrt(lambdaex.^2+mu.^2-2.*(sqrt(mu.^2.*(lambdaex.^2+lambdar.^2))-(lambdaex.^2*lambdar.^2)));
zetap=sqrt(lambdar.^4+k3p.^2.*v^2.*(lambdar.^2+lambdaex.^2));
zetan=sqrt(lambdar.^4+k3n.^2.*v^2.*(lambdar.^2+lambdaex.^2));

%k3p=abs(k3p1); k3n=abs(k3n1); zetap=abs(zetap1); zetan=abs(zetan1);

%Intermediate terms
t1= (lambdar.^2*v.^2)./(lambdar.^2+lambdaex.^2);
t2= ((lambdar.^2*lambdaex.^2)*(2*lambdar.^2+lambdaex.^2))./((lambdar.^2+lambdaex.^2).^2);

%Defining the regions
cshe1p= (e/(8*pi)).*t1.*((k3p.^2./zetap)-(k3n.^2./zetan))- (e/(4*pi)).*t2.*((1./zetap)-(1./zetan));
cshe1n= -(e/(8*pi)).*t1.*((k3p.^2./zetap)-(k3n.^2./zetan))+ (e/(4*pi)).*t2.*((1./zetap)-(1./zetan));
cshe2p= (e/(8*pi)).*t1.*(k3p.^2./zetap.^2)- (e/(4*pi)).*t2.*((1./zetap)-(1./lambdar.^2));
cshe2n= -(e/(8*pi)).*t1.*(k3p.^2./zetap.^2)+(e/(4*pi)).*t2.*((1./zetap)-(1./(lambdar.^2)));
cshe3p= (e/(8*pi)).*t1.*((k3p.^2./zetap)-(k3n.^2./zetan))- (e/(4*pi)).*t2.*((1./zetap)-(1./zetan));
cshe3n= -(e/(8*pi)).*t1.*((k3p.^2./zetap)-(k3n.^2./zetan))+(e/(4*pi)).*t2.*((1./zetap)-(1./zetan));
cshe4=0;
ax2=axes;

%Final equation
%cshe= cshe1p+cshe1n+cshe2p+cshe2n+cshe3p+cshe3n+cshe4;
cshe1=cshe1p+cshe2p+cshe3p+cshe4;
cshe=cshe1n+cshe2n+cshe3n+cshe4;

figure(1); %Allows working with multiple figures simultaneously
%surf(ax2,mu,lambdaex,((cshe)*(2))/(8*pi)) %mesh plot
surf( mu,lambdaex, cshe/(4*pi))
hold on
%surf(ax2,mu,lambdaex, ((cshe1)*(2))/(8*pi))
surf( mu,lambdaex, cshe1/(4*pi))
colorbar %add colorbar in plot
shading interp
colormap('turbo')
%set(gca,'zlim',[-1.5 1.5])
box on

[Dyuman Joshi]

There is a mistake in the definition of k3p and k3n, the parenthesis are not ordered correctly. Also, you should use element wise power, similar to the element wise multiplication you have done, not sqrt.

Also, you should define the cshe according to the values of mu, as shown in fig 2.

You are simply calculating cshe using all formulae for all values of mu and lambdaex at once, which is not what is done in the reference.

[SWASTIK SAHOO]

I have modified the code as per your comment. But still I am not getting the desired figures. Please let me know, where I am going wrong. I have attached the code here.

lambdar= 0.010;        % Rashba parameter (10 meV)
hbar=6.582e-16;    %Reduced Plancks constant in meV-s
vf= 0.812e-6;       % Fermi velocity in m/s .
v=hbar*vf; 
e=1; 

mu1= linspace(-0.001, 0.001, 5);    %Fermi potential
lambdaex1= linspace(0, 0.020, 5);  %Fixed Exchange parameter of 20meV

[mu,lambdaex]= meshgrid(mu1, lambdaex1);

% Defining two notations
k3p=(1/v)*sqrt((lambdaex.^2+mu.^2)+2*(sqrt(mu.^2.*(lambdaex.^2+lambdar.^2))-(lambdaex.^2*lambdar.^2)));
k3n=(1/v)*sqrt((lambdaex.^2+mu.^2)-2*(sqrt(mu.^2.*(lambdaex.^2+lambdar.^2))-(lambdaex.^2*lambdar.^2)));
zetap=sqrt(lambdar.^4+k3p.^2.*v^2.*(lambdar.^2+lambdaex.^2));
zetan=sqrt(lambdar.^4+k3n.^2.*v^2.*(lambdar.^2+lambdaex.^2));

%k3p=abs(k3p1); k3n=abs(k3n1); zetap=abs(zetap1); zetan=abs(zetan1);

%Intermediate terms
t1= (lambdar.^2*v.^2)./(lambdar.^2+lambdaex.^2);
t2= ((lambdar.^2.*lambdaex.^2)*(2*lambdar.^2+lambdaex.^2))./((lambdar.^2+lambdaex.^2).^2);
t3= sqrt((lambdar.^2.*lambdaex.^2)/(lambdar.^2+lambdaex.^2));

%Defining the regions

if (mu) > sqrt(4*lambdar.^2+lambdaex.^2)
    cshe1p= (e/(8*pi)).*t1.*((k3p.^2./zetap)-(k3n.^2./zetan))- ((e/(4*pi)).*t2.*((1./zetap)-(1./zetan)));
    cshe1n= -(e/(8*pi)).*t1.*((k3p.^2./zetap)-(k3n.^2./zetan))+ ((e/(4*pi)).*t2.*((1./zetap)-(1./zetan)));
    figure(1)
    surf( mu,lambdaex, cshe1p/(4*pi))
    hold on
   % surf( mu,lambdaex, cshe1n/(4*pi))
    hold on
elseif sqrt(4*lambdar.^2+lambdaex.^2) > (mu) > lambdaex
    cshe2p= (e/(8*pi)).*t1.*(k3p.^2./zetap.^2)- (e/(4*pi)).*t2.*((1./zetap)-(1./lambdar.^2));
    cshe2n= -(e/(8*pi)).*t1.*(k3p.^2./zetap.^2)+(e/(4*pi)).*t2.*((1./zetap)-(1./(lambdar.^2)));
    figure(1)
    surf( mu,lambdaex, cshe2p/(4*pi))
    hold on
  %  surf( mu,lambdaex, cshe2n/(4*pi))
    hold on
elseif lambdaex > (mu) > t3
    cshe3p= (e/(8*pi)).*t1.*((k3p.^2./zetap)-(k3n.^2./zetan))- (e/(4*pi)).*t2.*((1./zetap)-(1./zetan));
    cshe3n= -(e/(8*pi)).*t1.*((k3p.^2./zetap)-(k3n.^2./zetan))+(e/(4*pi)).*t2.*((1./zetap)-(1./zetan));
    figure(1)
    surf( mu,lambdaex, cshe3p/(4*pi))
    hold on
   % surf( mu,lambdaex, cshe3n/(4*pi))
    hold on
else -t3 < (mu) < t3
    cshe4=0;
    figure(1)
    surf( mu,lambdaex, cshe4/(4*pi))
end
hold off

colorbar %add colorbar in plot
shading interp
colormap('turbo')
%set(gca,'zlim',[-1.5 1.5])
box on

Answer 1

[Balaji]

Hi Swastik,

As per my understanding you have difficulties in implementing the results of the paper you mentioned.

There are a couple of errors in your code.

Firstly, in MATLAB, expressions like a<b<c are interpreted as (a<b)<c. you can use (a<b) & (b<c) instead.

Secondly, comparison of two matrices gives a Boolean matrix which compares corresponding elements of the matrices.

So, in order to get your desired results you can use the following instead of the if else statements

cshe1p= (e/(8*pi)).*t1.*((k3p.^2./zetap)-(k3n.^2./zetan))- ((e/(4*pi)).*t2.*((1./zetap)-(1./zetan)));
cshe1n= -(e/(8*pi)).*t1.*((k3p.^2./zetap)-(k3n.^2./zetan))+ ((e/(4*pi)).*t2.*((1./zetap)-(1./zetan)));
cshe2p= (e/(8*pi)).*t1.*(k3p.^2./zetap.^2)- (e/(4*pi)).*t2.*((1./zetap)-(1./lambdar.^2));
cshe2n= -(e/(8*pi)).*t1.*(k3p.^2./zetap.^2)+(e/(4*pi)).*t2.*((1./zetap)-(1./(lambdar.^2)));
cshe3p= (e/(8*pi)).*t1.*((k3p.^2./zetap)-(k3n.^2./zetan))- (e/(4*pi)).*t2.*((1./zetap)-(1./zetan));
cshe3n= -(e/(8*pi)).*t1.*((k3p.^2./zetap)-(k3n.^2./zetan))+(e/(4*pi)).*t2.*((1./zetap)-(1./zetan));

condition1 = (mu) > sqrt(4*lambdar.^2+lambdaex.^2);
condition2 = (sqrt(4*lambdar.^2+lambdaex.^2) > (mu)) & ((mu) > lambdaex);
condition3 = (lambdaex > (mu)) & ((mu) > t3);
condition4 = (-t3 < (mu)) & ((mu) < t3);

cshep = zeros(5);
cshep(condition1) = cshe1p(condition1);
cshep(condition2) = cshe2p(condition2);
cshep(condition3) = cshe3p(condition3);
cshep(condition4) = 0;

Hope this helps!

Thanks,

Balaji