3
$\begingroup$

I am trying to plot the density of states $$ N(E)= \frac{1}{N}\Sigma_{k} (\delta{(E - \epsilon_{k}))} $$ where $\epsilon_{k} = -2t*[cos(k_x)+cos(k_y)]$ and $k$ goes from $-\pi$ to $\pi$. For this plot, we will get a singular point when E = 0 and it will peak at this point. This is done for a 2D Hubbard Model without considering the disorder. This is the plot I am expected to get -

enter image description here

My code:

energy[kx_, ky_] := -2  t  (Cos[kx] + Cos[ky]);

nPoints = 100;
kRange = Range[-Pi, Pi, 2  Pi/(nPoints - 1)];

energyValues = 
  Flatten[Table[energy[kx, ky], {kx, kRange}, {ky, kRange}]];

densityOfStates[
   Ee_] := (1/Length[energyValues])  Sum[
    DiracDelta[Ee - en], {en, energyValues}];

Plot[densityOfStates[Ee], {Ee, -4, 4}, PlotRange -> All, 
 AxesLabel -> {"E", "N(E)"}, PlotRange -> {{-4, 4}, {0, 0.5}}]

This is the plot I am getting -

enter image description here

I am implementing this paper - https://www.cond-mat.de/events/correl16/manuscripts/scalettar.pdf

Improve this question
$\endgroup$
5
  • $\begingroup$ DiracDelta is the implementation of the $\delta$-distribution in Mathematica, not a usual function. Its plots make no sense. $\endgroup$
    – user64494
    Commented Apr 27 at 12:59
  • $\begingroup$ @user64494 Oh okay, could you please suggest any other way I can implement this? Thanks $\endgroup$
    – CuriousMind
    Commented Apr 27 at 13:15
  • $\begingroup$ What is t? How can you plot anything when you miss the definition of t? $\endgroup$
    – azerbajdzan
    Commented Apr 27 at 13:52
  • $\begingroup$ @azerbajdzan I set t to 1, its the hopping parameter of the particle between two any adjacent sites. $\endgroup$
    – CuriousMind
    Commented Apr 27 at 13:54
  • $\begingroup$ In the paper they say " In the continuum limit (large number of sites), the sum over discrete momenta is replaced by an integral". So if you know how to do it maybe you can get a better approximation than I got using sum. $\endgroup$
    – azerbajdzan
    Commented Apr 27 at 15:04

2 Answers 2

Reset to default
3
$\begingroup$

I used function dirac as approximation to Dirac delta. Dirac delta is limit as a->0. I used value a->1/3 and nPoints = 30. For larger number of points the code is slow and for smaller value of a there are overflows errors.

The plot is quite similar but still not same as in OP. So I guess there might be some mistakes in OP code or they used in the paper a different approximation for Dirac delta.

t = 1.;
dirac[x_] := 1/(a Sqrt[Pi]) E^-(x/a)^2 /. a -> 1/3

energy[kx_, ky_] := -2 t (Cos[kx] + Cos[ky]);

nPoints = 30;
kRange = Range[-Pi, Pi, 2 Pi/(nPoints - 1)];

energyValues = 
  Flatten[Table[energy[kx, ky], {kx, kRange}, {ky, kRange}]];

densityOfStates[
   Ee_] := (1/Length[energyValues]) Sum[
    dirac[Ee - en], {en, energyValues}];

Plot[densityOfStates[Ee], {Ee, -4, 4}, PlotRange -> All, 
 AxesLabel -> {"E", "N(E)"}, PlotRange -> {{-4, 4}, {0, 0.5}}]

enter image description here

With a->1/10 in dirac[x_] := 1/(a Sqrt[Pi]) E^-(x/a)^2 /. a -> 1/10 and nPoints = 200 the plot resembles the OP image quite well. For smaller a and larger nPoints I am sure it can be even better but the code is slow.

enter image description here

Improve this answer
$\endgroup$
5
  • $\begingroup$ This is wishful thinking. Replacing a->1/3 by a->1/10, one obtains something strange. For the user's convenience I add the result to your answer. $\endgroup$
    – user64494
    Commented Apr 27 at 17:47
  • 5
    $\begingroup$ Can moderators block this user to stop modifying my answer with useless edits? $\endgroup$
    – azerbajdzan
    Commented Apr 27 at 17:51
  • $\begingroup$ See the screen with my addition to the azerbajdzan's answer. $\endgroup$
    – user64494
    Commented Apr 27 at 17:58
  • $\begingroup$ To those who can not read with comprehension. The code depends also on number of nPoints not only on approximation factor a of Dirac delta. When a is decreased the number of points nPoints must be increased appropriately. $\endgroup$
    – azerbajdzan
    Commented Apr 27 at 18:11
  • $\begingroup$ Indeed, with a->1/10 and nPoints=100; the result is better than yours. I repeat that a good code is a commented code. $\endgroup$
    – user64494
    Commented Apr 27 at 18:21
2
$\begingroup$

This is an extended comment to support the answer by @azerbajdzan (and his comments). I have voted for @azerbajdzan answer. Using the dirac delta distribution approximation. As @azerbajdzan comments it depends on choices. It is a different implementation (probably slower) but the motivation is illustrative. The point is better approximations of the density of states for finer meshes require better delta distribution approximations.

d[x_, e_] := Exp[-(x/e)^2]/(e   Sqrt[Pi])
mesh[n_] := -2 (Cos[#1] + Cos[#2]) & @@@ 
  Tuples[Subdivide[-Pi, Pi, n - 1], 2]
func[en_, e_, n_] := Chop[Total[d[en - #, e] & /@ mesh[n]]/n^2]
gr = Show[
   Quiet@Table[
     ListPlot[Table[{j, Evaluate@func[j, k, 5/k]}, {j, -5, 5, 0.1}], 
      Joined -> True, GridLines -> {{-4, 4}, None}, 
      PlotRange -> {0, 0.4}, 
      PlotStyle -> RGBColor[10 k, 0, 0]], {k, {0.1, 0.05, 0.01}}]];


Legended[gr,
 LineLegend @@ 
  Transpose[
   Table[{RGBColor[10 k, 0, 0], 
     Row[{"e = ", k, ", grid size = ", 25/k^2}]}, {k, {0.1, 0.05, 
      0.01}}]]]

enter image description here

Improve this answer
$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.