# Alternative Representation for Energy Levels and Energy Level Density (II)

Following my previous question, I want a graph that must look like this (the lower image is taken from the Wikipedia page for the Microcanonocal Ensemble-->> Quantum Mechanical),

The code for solving Schrodinger's equation:

a0 = 0.02; a1 = 0.64; L = 10;
U[x_] := a0*x^4 - a1*x^2 + a1^2/(4 a0)
{vals, funs} =
NDEigensystem[{-1/2 Laplacian[u[x], {x}] + U[x]*u[x],
DirichletCondition[u[x] == 0, True]}, u[x], {x, -L, L}, 25,
Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions" \
-> {"MaxCellMeasure" -> 0.001}}}}];


How do we plot the graphs showing level density as shown in the above figure?

N.B.

The author used Python to plot the graph shared above. Please see the Source.

## 1 Answer

You can use ColorFunction to plot the wavefunction densities:

numEig = 25;
Show[Plot[U[x], {x, -2 L/3, 2 L/3}, PlotStyle -> Red,
AspectRatio -> 1.2, PlotHighlighting -> None],
Table[Plot[{vals[[i]]}, {x, -10, 10}, PlotStyle -> Thick,
ColorFunction -> Function[{x, y}, Opacity[2 funs[[i]][x]^2, Blue]],
ColorFunctionScaling -> False, PlotHighlighting -> None], {i, 1,
numEig}]]


As presented in the previous answer, use SmoothHistogram for the density of states. To rotate the plot, see #145789.

SmoothHistogram[vals, 0.1]


• Thank you @Domen for this quick response. the first graph is absolutely beautiful. However, I believe, the second graph is not correct. The adjacent graph in the Wikipedia page is a probability graph with inverted axes, i.e. the x-axis shows the state probability while y shows its position on the double well as a function of energy. Feb 29 at 10:53
• Please see this Python code for the above Wikipedia graph : commons.wikimedia.org/wiki/… Feb 29 at 11:02
• @user444, the plot in your question shows just the density of states, $g(E)$, and that is shown on my plot . The plot in the link you've pasted now shows a canonical ensemble, where you weight the states by the Boltzmann factor to get $\rho(E) = g(E) e^{-\beta E}$. So the two plots are different. Feb 29 at 11:05
• For the canonical ensemble, you could probably use something like β = 1/10; SmoothHistogram[WeightedData[vals, Exp[-β vals]], 0.1]. Feb 29 at 11:33
• Thank you @Domen, for you kind help. Feb 29 at 12:12