Question: Do the divergences vanish or not?

Load the Tetrads package, and choose Schwarzshild spacetime:

with(Physics):
with(Tetrads):
Setup(metric = schwarzschild,quiet):

Consider then the following code/output:

expr := D_[~mu](e_[mu,a]);
seq(simplify(SumOverRepeatedIndices(expr             )),a = 1..4);
seq(simplify(SumOverRepeatedIndices(D_[~mu](e_[mu,a]))),a = 1..4);

In my mind, the last two code lines are entirely identical, but their corresponding outputs are certainly not. I don't get it. The above can be compared with the following code/output in which the vierbein indices mu and a have been interchanged:

expr := D_[~mu](e_[a,mu]);
seq(simplify(SumOverRepeatedIndices(expr             )),a = 1..4);
seq(simplify(SumOverRepeatedIndices(D_[~mu](e_[a,mu]))),a = 1..4);

Here, the outputs are identical. The mystery becomes total when realizing that according to Maple, e_[a,mu] and e_[mu,a] are considered completely identical [a notational choice I do not find reassuring].