since this thread got some recent recognition, I'll post the answer that I found out some time ago through some additional literature research and the one I personally was most satisfied with.
As mze's answer proposed, the PREM model is used as the density profile for the matter potential and an analysis in Joachim Kopp's diploma thesis in chapter 4.2.4 shows that a simple description with an average density does not provide the same results as a better Mantle-Core-Mantle model with three different average regions or even the full PREM model (see e.g. fig. 4.6).
Kopp also develops an analytic expression (see eq. 4.13), which is quite lengthy.
In the end, i chose to calculate the appearance probability of neutrinos going completely through the earth numerically in the three-layer Mantle-Core-Mantle model with $\rho_\text{mantle}=4.66 \frac{\text{g}}{\text{cm}^3}, \rho_\text{core}=11.8 \frac{\text{g}}{\text{cm}^3}$, via
\begin{align}
P_{\nu_\mu \to \nu_e}(E)=|\langle \nu_e | \exp{(-\mathrm{i} H_\text{MCM}(E, \bar{x}_1) \cdot x_1)} \exp{(-\mathrm{i} H_\text{MCM}(E, \bar{x}_2) \cdot x_2)} \exp{(-\mathrm{i} H_\text{MCM}(E, \bar{x}_3)\cdot x_3)}| \nu_\mu \rangle|^2
\end{align}
, where in flavor space we have
\begin{align}
\rho_\text{MCM}(x) &= \left\lbrace \array{ \rho_\text{mantle}, & 0 \leq x \leq 0.5 \\
\rho_\text{core}, & 0.5 < x \leq 1.5 \\
\rho_\text{mantle}, & 1.5 < x \leq 2 \\
0 & \text{else} } \right. \, , \\
V(E,x) &= \frac{1}{2E} \begin{pmatrix} 2 \sqrt{2} \frac{Y G_F }{m_N} c_\text{conv} \rho(\frac{x-R_\text{Earth}}{R_\text{Earth}}) E & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix} \, ,\\
H_\text{MCM}(E, x) &= \frac{1}{2E} U_\text{PMNS} \begin{pmatrix} 0 & 0 & 0 \\ 0 & \Delta m^2_{21} & 0 \\ 0 & 0 & \Delta m^2_{31}\end{pmatrix} U_\text{PMNS}^\dagger + V(E,x) \, .
\end{align}
In this case for the first layer we have the layer thickness $x_1=0.5R_\text{Earth}$ and $\bar{x}_1 \in [0,0.5R_\text{Earth}]$ to set the Hamiltonian to pick the correct average density and analogously $x_{2}= R $, $\bar{x}_2 \in (0.5R_\text{Earth}, 1.5R_\text{Earth}]$ and $x_3=0.5R_\text{Earth}$ as well as $\bar{x}_3=(1.5R_\text{Earth}, 2R_\text{Earth}]$. In the case of the PREM model, one would have to discretize to even more layers to get more and more accuracy.
I also assumed the average electron number per nucleon to be $Y \approx 0.5$ and the mass of a nucleon at about $m_N \approx 1\,$GeV, so you can see that I use natural units. In the spirit of this unit system, don't forget the conversion factor $ c_\text{conv}$ from $\frac{\text{g}}{\text{cm}^3}$ to $\text{eV}^4$.
Of course, when you want to incorporate sterile neutrinos, you'll have to add the neutral current contributions from the aforementioned paper by Linder to the sterile slots in the potential $V(E, x)$.
This reproduces the results found by Kopp to a good degree.
Cheers!
