The answer I'm posting (after two years) is based on an idea by @dominicprice (former student of Kasper, and developer of many aspects of `cadabra`

).

# The code

```
def expand_nabla(ex, *, index_name=None, index_set=None):
if index_name is None:
index_name = "expandnablaidx"
Indices($expandnablaidx, expandnablaidx#$, $expandnablaindices$)
index_counter = 1
for nabla in ex[r'\nabla']:
nabla.name=r'\partial'
dindex = nabla.indices().__next__()
for arg in nabla.args():
ret:=0;
for index in arg.free_indices():
t2:= @(arg);
newidx = Ex(index_name + str(index_counter))
index_counter += 1
if index.parent_rel==sub:
t1:= -\Gamma^{@(newidx)}_{@(dindex) @(index)};
t2[index]:= _{@(newidx)};
else:
t1:= \Gamma^{@(index)}_{@(dindex) @(newidx)};
t2[index]:= ^{@(newidx)};
ret += Ex(str(nabla.multiplier)) * t1 * t2
nabla += ret
distribute(ex, repeat=True)
if index_set is not None:
rename_dummies(ex, "expandnablaindices", index_set)
return ex
```

# How does it work?

The problem I found with the algorithm `expand_nabla`

reported in the cadabra book https://cadabra.science/notebooks/ref_programming.html, was that when considering higher order derivatives, the name of the dummy index wouldn't change, resulting in a "repeated indices" error.

It is desirable that the algorithm could use one of the declared indices!

Therefore, we have to give a name to the set of indices:

```
\nabla{#}::Derivative.
\partial{#}::PartialDerivative.
{a, b, c, d, e, f, g, h, i, j, k ,l}::Indices(space, position=fixed).
{\mu,\nu,\lambda,\rho}::Indices(spacetime, position=fixed).
```

**NOTE:** there are two types of indices: `space`

and `spacetime`

.

Now, we declare an expression and expand the nabla operator (covariant derivative).

```
foo := \nabla_{c}{ h^{a}_{b} };
expand_nabla(_, index_set="space");
```

Since we have asked to expand with dummy indices of type `space`

, the result is

$$\partial*{c}{h^{a},*{b}}+\Gamma^{a},*{c d} h^{d},*{b}-\Gamma^{d},*{c b} h^{a},*{d}$$

However, if we expand indices of type `spacetime`

, we get

$$\partial*{c}{h^{a}\,*{b}}+\Gamma^{a}\,*{c \mu} h^{\mu}\,*{b}-\Gamma^{\mu}\,*{c b} h^{a}\,*{\mu}$$

## Higher orders

Higher orders (I was able to calculate up to third order) can be computed with ease, for example:

```
foo := \nabla_{e}{\nabla_{d}{\nabla_{c}{ h^{a}_{b} }}};
#expand_nabla(_, index_set="vector");
expand_nabla(_, index_name=r'\lambda');
```

The `index_name=r'\lambda'`

ask the algorithm to use the name `\lambda`

to expand the dummy indices.