Dear Edgardo,

great, many thanks. We've had a few emails about this back and forth in the meantime, also about generalisations where there are coefficients c[n] in the Hamiltonian, and I am very happy with the support you have included in the Physics package. As of v. 1200 it does exactly what I was hoping to do here. I consider this issue solved.

As I have mentioned above, this kind of detailed support is unheard of from other big commercial software, and really singles out Maple (and in particular the Physics package) from my perspective.

Cheers, Daniel

@ecterrab Thanks a lot Edgardo, I really appreciate the attention to detail and the time you've spent even updating the Physics package. That is a very rare treat in the software world and goes to show the strengths of Maple wrt other products.

A minor comment: I've upgraded the physics package and my Eq. (2) does not look like yours: mine still contains 4 KroneckerDeltas, even the two for the diagonal terms. Only after I say "Expand" and "Simplify" I get your Eq. (3).

To give some physics context for the whole question, I want to generate the right hand side of Eq. (3) from 10.1103/PhysRevLett.119.030402 using computer algebra.

I completely agree that the difficulty of simplifcation comes from the boundary terms here. While your suggested fix works, it is somewhat "manual". Ultimately, I need to look at higher powers of H, and having to stare at each term by hand and realising which boundary terms do and do not contribute is exactly what I want to avoid in the first place. So I wonder if there's a more elegant way to do this:

1) can we tell Maple to begin with that the dimension of Hilbert space is d+1 (I guess I should have posed the question originally to go only to d-1, so that the Hilbert space dimension would be d, but ok let's keep it d+1). Then Maple realises which boundary terms do and do not exist?

2) I was trying to add coefficients c[n] to the summation above, and then tell Maple that c[0]=0, c[d+1]=0. This way, it should be easier for Maple to see which terms do not contribute, even if it does not know the Hilbert space dimension. But I had (beginners) problems telling Maple that the c[n] are realnumbers. Obviously Setup(realobjects = {c[n]}) does not work, what would be the correct syntax here?