Hopefully this is a quick fix, but I've been working on this for 3 days to no avail, so perhaps some external perspective will help.

Basically, I'm building a MatLab module to calculate mean particle $\frac{dE}{dx}$ from the Bethe-Bloch equation. To verify that I'd implemented it correctly, I used the particle Data Group chapter on "Particle Interaction with matter" (see link at end for a pdf thereof). Unfortunately, as you can see in the following plot, my implementation does not match the data very well at all.

enter image description here

The PDG data was digitized from the widely available plot seen below:

enter image description here

To generate the blue data series in the first plot, I merely calculated the following quantities for a range of (total) proton energies $E_n$ in MeV:

$$\gamma_n = \frac{E_n}{E_0}$$

where $E_0$ is the rest energy of the incident particle (proton) in MeV.

$$\beta_n = \sqrt{1-\bigg(\frac{1}{\gamma_n}\bigg)^2}$$

$$p = \sqrt{E_n^2 - \bigg(\frac{E_0}{c^2}\bigg)^2} = \frac{E_0}{c^2} \beta \gamma$$

$$T_{max}=\cfrac{2 m_e c^2 \beta^2 \gamma^2}{1+\cfrac{2 \gamma m_e}{M}+\bigg(\cfrac{m_e}{M}\bigg)^2}$$

Where $p$ is the particle momentum, $M$ is the rest mass of the incident proton in g, $m_e$ is the rest mass of an electron in g,and $c$ is, of course, the speed of light in m/s. I calculated the density effect, $\delta$ with the Sternheimer parameterization where:

$$x = \log{\beta \gamma}$$ $$ \delta(\beta \gamma) = \begin{cases} 2(\ln 10)x -\bar{C} & \text{if }x\geq x_1; \\ 2(\ln 10)x -\bar{C}+a(x_1-x)^k & \text{if }x_0\leq x < x_1; \\ 0 & \text{if }x<x_0\ \text{(nonconductors)}; \\ \delta_010^{2(x-x_0)} & \text{if }x<x_0\ \text{(conductors)}; \\ \end{cases} $$

Finally, the Bethe-Bloch equation I used (straight from the PDG pdf) was as follows:

$$-\frac{dE}{dx} = K z^2 \frac{Z_{target}}{A_{target}} \frac{1}{\beta^2} \bigg[\frac{1}{2} \ln\frac{2 m_e c^2 \beta^2 \gamma^2 Tmax}{I^2}-\beta^2 - \delta/2\bigg]$$

For lead, the following material parameters were used: $Z_{target} = 82$ amu, $A_{target}=207.2$ g/mol, $I= 823$ eV. For the calculation of $\delta$, the following parameters were used (see link below for source): $C=6.2018$, $X_0=0.3776$, $X_1=3.8073$, $a=0.09359$, $k=3.1608$,and $\delta_0=0.14$. $K$ was taken to be constant at a value of 0.307075 $\frac{MeV mol}{g}$.

Unfortunately I still get the error evident in the first plot where the Fermi Plateau doesn't seem to be in evidence. Does anybody see anything wrong in the steps I took to calculate $\frac{dE}{dx}$?

For the record, this is NOT a problem with programming/MatLab, my trouble lies with the computation of $\frac{dE}{dx}$. The issue is simply that the Bethe Bloch equation, when computed, does not appear to reflect the behavior seen in the PDG plot from literature.

PDG Link: http://pdg.lbl.gov/2009/reviews/rpp2009-rev-passage-particles-matter.pdf

Sternheimer Parameterization Link: https://www.sciencedirect.com/science/article/pii/0092640X84900020

  • $\begingroup$ Have you open sourced your MATLAB module ? I'd to take a look at the source code so I can gain some insight into calculating $dE/dx$ $\endgroup$
    – Zophikel
    May 30, 2019 at 14:56

1 Answer 1


Murhphy's law. Found out what was wrong as soon as I posted to the exchange. my $I$ was in the wrong units (eV instead of MeV).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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