In[1]:= Quiet @ Needs["HighEnergyPhysics`FeynCalc`"]

Loading FeynCalc from C:\Users\Frank\AppData\Roaming\Mathematica\Applications\HighEnergyPhysics

FeynCalc 8.2.0 For help, type ?FeynCalc, open FeynCalcRef8.nb or visit www.feyncalc.org

Loading FeynArts, see www.feynarts.de for documentation

FeynArts 3.7 patched for use with FeynCalc

In[2]:= DeclareNonCommutative[a, ad]

In[3]:= Commutator[a, ad] = 1

Out[3]= 1

In[4]:= DotSimplify[ Commutator[a, ad.ad]]

Out[4]= 2 ad

In[5]:= DotSimplify[ Commutator[a, ad.ad.ad]]

Out[5]= 3 ad.ad

In[6]:= DotSimplify[ Commutator[a, ad.ad.ad.ad]]

Out[6]= 4 ad.ad.ad

Looks like FeynCalc is the way to go.

Here's an example problem from an online course I'm taking:

Recalling that [a^,a^†]=1, calculate the commutator:

[a^,(a^†)n]

a is an annihilation operator and a dagger is a creation operator.