I'm trying to get the Riemann tensor with four covariant indices from the standard form given by,
$$R^l_{ikj}=\Gamma^l_{ij,k}-\Gamma^l_{ik,j}+\Gamma^m_{ij}\Gamma^l_{mk}-\Gamma^m_{ik}\Gamma^l_{mj}$$
Contracting with the metric, we have
$$R_{likj}=g_{ls}R^s_{ijk}=g_{ls}(\Gamma^s_{ij,k}-\Gamma^s_{ik,j}+\Gamma^m_{ij}\Gamma^s_{mk}-\Gamma^m_{ik}\Gamma^s_{mj})$$
Using the identity
$$g_{ls}\Gamma^s_{ij,k}=\Gamma_{ijl,k}-\Gamma^s_{ij}(\Gamma_{lks}+\Gamma_{skl})$$
it becomes
$$R_{likj}=\Gamma_{ijl,k}-\Gamma_{ikl,j}-\Gamma^s_{ij}(\Gamma_{lks}+\Gamma_{skl})+\Gamma^s_{ik}(\Gamma_{ljs}+\Gamma_{sjl})+\Gamma^m_{ij}\Gamma_{lmk}-\Gamma^m_{ik}\Gamma_{lmj}$$
I should arrive to,
$$R_{likj}=\Gamma_{ijl,k}-\Gamma_{ikl,j}+\Gamma^m_{ik}\Gamma_{ljm}-\Gamma^m_{ij}\Gamma_{lkm}$$
as given here: https://www.cmu.edu/biolphys/deserno/pdf/diff_geom.pdf
However, I have no idea on how to combine the last four terms to give this, as I already have the first two terms. I've checked several pages here and on the web but so far I haven't found any proper proof using only indices. I've checked the following pages:
Deriving Riemann tensor with four lower indices
https://www.zweigmedia.com/diff_geom/Sec10.html
https://www.youtube.com/watch?v=UVBxE4IJKHs
https://physics.stackexchange.com/questions/378009/lower-raising-riemann-curvature-tensor
https://physics.stackexchange.com/questions/485583/lowering-index-of-riemann- tensor
What is the definition of $R_{ijkl}$ in terms of metrics on a manifold?
Riemann tensor with all indices lowered on a 2D manifold https://physics.stackexchange.com/questions/484815/indices-of-the-riemann-tensor-of-the-first-kind
I also found someone try a proof here: https://physics.stackexchange.com/questions/309022/covariant-riemann-tensor-indices
However, he uses some weird identities such as
$$g_{ij,k}=\Gamma_{kij}+\Gamma_{kji}$$
which is false, since we can directly proof that,
$$\Gamma_{kij}+\Gamma_{kji}=g_{ki,j}+g_{kj,i}-g_{ij,k}$$
I really appreciate any help, I've been trying to prove this for over a month now, so this is my last resource.
Edit: one can also find this formula in Wikipedia (https://en.wikipedia.org/wiki/List_of_formulas_in_Riemannian_geometry) as:
$$R_{ilkm}=\frac{1}{2}(g_{im,kl}+g_{kl,im}-g_{il,km}-g_{km,il})-g_{np}(\Gamma^n_{kl}\Gamma^p_{im}-\Gamma^n_{km}\Gamma^p_{il})$$
which is analogous to the one shown above, just with an added metric tensor to have mixed indices in the Christoffel symbols (the expansion in terms of the metric tensor follows trivially from expanding the first two terms of Christoffel symbols with derivatives).