Let me first re-write your notations as little bit to make it easier for me.
Let $\xi\in \Xi \subset\mathbb{R}^n$. Define the functions $\pi_k:\Xi \to\mathbb{R}$ by
$$ \pi_k(\xi) = \begin{cases} \xi_k & k \in \{1, \dots, n\} \\ 1 - \sum_{1}^n \xi_i & k = 0 \end{cases} $$
which to me is a more natural way to write your function $p_\xi: \{0, \ldots, n\} \to \mathbb{R}$.
Your metric is defined to be
$$ \mathbf{g} = \sum_{k = 0}^n \frac{1}{\pi_k} ~\mathrm{d}\pi_k \otimes \mathrm{d}\pi_k = \sum_{k = 0}^n \mathrm{d}\pi_k \otimes \mathrm{d} \log \pi_k $$
What you write as the Christoffel symbols is not really the Christoffel symbol, it is the Christoffel symbol with one index lowered. Note that the Christoffel symbol for the Levi-Civita connection is given by
$$ g_{k\ell} \Gamma^{\ell}_{ij} = \frac12 \left( \partial_i g_{kj} + \partial_j g_{ki} - \partial_k g_{ij}\right) $$
one check indeed
$$ g_{k\ell} \Gamma^{\ell}_{ij} = \sum_{m = 0}^n \partial_k \pi_m \partial^2_{ij} \log\pi_m = - \sum_{m = 0}^n \frac{1}{\pi_m^2} \partial_k \pi_m \partial_i \pi_m \partial_j \pi_m $$
where in the second equality we used that $\pi_m$ are linear in the coordinates and hence the Hessian vanishes.
Let $\gamma:\mathbb{R}\to \Xi$ be a geodesic. The geodesic equation in coordinates can be written as
$$ g_{k\ell}\circ\gamma~ \ddot{\gamma}^\ell + g_{k\ell}\circ\gamma~ \Gamma^\ell_{ij}\circ\gamma ~\dot{\gamma}^i \dot{\gamma}^j = 0$$
plugging in the definitions we have
$$ \sum_{m = 0}^n (\partial_k \pi_m)\circ\gamma \left[ (\partial_\ell \log \pi_m) \circ\gamma ~\ddot{\gamma}^{\ell} + (\partial^2_{ij} \log\pi_m)\circ\gamma~ \dot{\gamma}^i \dot{\gamma}^j \right] = 0 $$
which we simplify as
$$ \sum_{m = 0}^n (\partial_k \pi_m) \circ \gamma \frac{d^2}{dt^2}( \log\circ\pi_m\circ \gamma) = 0 $$
Now, $\partial_k \pi_m = \delta_{km}$ if $m > 0$ and $-1$ otherwise. So we have
$$ \frac{d^2}{dt^2}( \log\circ\pi_k\circ \gamma) - \frac{d^2}{dt^2}(\log\circ\pi_0 \circ \gamma) = 0 \tag{1}$$
This is almost what you want. To finish, we observe the following:
$$ \sum_{m = 0}^n \pi_m\circ \gamma = 1 \implies \sum_{m = 0}^n \frac{d^k}{dt^k} \pi_m\circ\gamma = 0, k \geq 1 \tag{2}$$
Observe that
$$ \frac{d^2}{dt^2} \log \circ\pi_m \circ \gamma = \frac{1}{\pi_m\circ\gamma} \frac{d^2}{dt^2} \pi_m\circ\gamma - \frac{1}{(\pi_m\circ\gamma)^2} (\frac{d}{dt} \pi_m\circ\gamma)^2 $$
This says that
$$ \sum_{m = 0}^n \pi_m\circ\gamma \frac{d^2}{dt^2} \log\circ\pi_m\circ\gamma = - \sum_{m = 0}^n \frac{1}{(\pi_m\circ\gamma)} (\frac{d}{dt} \pi_m\circ\gamma)^2 \tag{3}$$
Applying (3) to (1) gives
$$ \sum_{m = 0}^n \pi_m\circ\gamma \frac{d^2}{dt^2} \log\circ\pi_0 \circ\gamma = - \sum_{m = 0}^n \frac{1}{(\pi_m\circ\gamma)} (\frac{d}{dt} \pi_m\circ\gamma)^2 $$
so from (2), and (1) again, we get, for every $k \in \{0,\dots, n\}$,
$$ \underbrace{\frac{d^2}{dt^2} \log\circ\pi_k \circ\gamma}_{\ddot{l}_t(x_k)} + \sum_{m = 0}^n \underbrace{\frac{d}{dt}(\pi_m\circ\gamma)}_{\dot{\gamma}_t(x_m)} \underbrace{\frac{d}{dt} \log\circ\pi_m\circ\gamma}_{\dot{l}_t(x_m)} = 0 $$
which is exactly the equation claimed in the question.
I should remark that part of the above answer is reconstructed in differential geometry language from the linked paper. In some ways the original derivation is a bit slicker, especially when it comes to the derivative of what I wrote as equation (3).
Addendum to address this final paragraph:
What you are treating is really the Fisher information metric. Consider the set $\Xi$ which is the interior of the simplex. If you look at the functions $\eta_i$ such that $(\eta_i)^2 = \xi_i$, the simplex is defined as $\sum (\eta_i)^2 = 1$, with $\eta_i \in (0,1)$. Notice that this can be interpreted as a quadrant of the unit sphere with $n$ dimensions (so $i$ runs from $0$ to $n$).
The metric you wrote down in $\xi$ coordinates turns to to be exactly the pull back metric of the standard unit sphere (up to some constant multiples) to $\Xi$ using the mapping $(\eta_i)^2 = \xi_i$.
The geodesics on the unit sphere are the great circles and are thus the intersections of the sphere with hyperplanes. Let $\zeta \in \mathbb{R}^{n+1}$ be a fixed vector, there is a geodesic in $\eta$ coordinates satisfying $\sum \zeta_i \eta_i = 0$. We can choose the affine parametrisation $\eta_i(t)$ to be such that $\sum (\frac{d}{dt}\eta_i(t))^2 = E$ is constant. It is well-known (since Newton) that in this parametrisation the acceleration of $\eta$ is in the direction of $\eta$ itself, with the proportionality factor $E$: that is
$$ \ddot{\eta} + E \eta = 0. $$
Doing the reverse change of variables you arrive back to the formula for the expression of the geodesics in the $\xi$ coordinates.
The fact that $\dot{\ell}$ in the language of the paper is a tangent vector, is exactly the statement that in $\eta$ coordinates, any curve along the sphere must satisfy $\sum_i \eta_i \dot{\eta}_i = 0$. By itself it has nothing to do with the fact that $\ell$ is geodesic. The "slick" derivation that I referred to is essentially what I wrote here, but cast in the language of information geometry where for some reason one prefers the coordinates $\xi$ instead of $\eta$.
