I'm a physicist studying higher-spin theory. In my research, we work with fully symmetric tensors using a notation which is implicit both in the dimension of space and the order of the tensor, i.e. an order-$s$ fully symmetric tensor is simply written as
$$ \phi_{\mu_1 \cdots \mu_s} \equiv \phi. $$
The $n$-th gradient of $\phi$ is written as $\partial^n \phi$, the $n$-th divergence, as $\partial^n \cdot \phi$ and the $n$-th trace as $\phi^{[n]}$. Lower traces are simply written with a prime, e.g. $\phi''$ for the second trace.
All indices are implicitly symmetrized, without weight factors, using the minimal number of terms.
For example, if $s=2$,
$$ \partial^2 \phi \equiv \partial_\mu \partial_\nu \phi_{\sigma \rho} + \partial_\mu \partial_\sigma \phi_{\nu \rho} + \partial_\mu \partial_\rho \phi_{\sigma \nu} + \partial_\nu \partial_\sigma \phi_{\mu \rho} + \partial_\nu \partial_\rho \phi_{\sigma \mu} + \partial_\sigma \partial_\rho \phi_{\mu \nu} $$
$$ \partial (\partial \cdot \phi) \equiv \partial_\nu (\partial^\lambda \phi_{\mu \lambda}) + \partial_\mu (\partial^\lambda \phi_{\lambda \nu})$$
$$ \eta \partial^2 \cdot \phi \equiv \eta_{\mu\nu} \partial^\rho \partial^\sigma \phi_{\rho \sigma} $$
The formalism implies the following set of rules
$$ ( \partial^p \phi )' = \Box \partial^{p-2} \phi + 2 \partial^{p-1} \left( \partial \cdot \phi \right) + \partial^p \phi' $$
$$ \partial \cdot (\partial^p \phi) = \Box \partial^{p-1} \phi + \partial^p \left( \partial \cdot \phi \right) $$
$$ \partial^p \partial^q = {{p+q}\choose{q}} \partial^{p+q} $$
which you can simply take as axioms of the game. By the way, $\Box = \partial \cdot \partial$, which we are allowed to invert to $\frac{1}{\Box}$, but don't worry about what this means, you can just use the three rules written above.
PROBLEM
As you can see, these are very simple objects and I was able to write a C++ code implementing these objects. Basically I just created a data type specifying a symbol for the tensor, its rank, the number of traces, the number of divergences, the number of gradients, the number of $\Box$-operators and the numerical multiplicative factor, along with some obvious algebraic constraints (for example, 3-divergence of a 2-tensor is trivially zero).
I also implemented functions on these objects, which take gradients, divergences or traces, along with the multiplication by numbers and addition of such objects.
I would really like to implement this in Mathematica, since I can't get everything I need out of C++, but I'm not sure where to start with defining these objects and functions on them. I don't want to use a pre-existing package, I'd like to build it from scratch, with as little unnecessary baggage as possible for full functionality. Any ideas? All tips are greatly appreciated.
Please do tell me if you want me to clarify anything.
EDIT 1
By lower traces I simply mean a number of traces low enough to actually write out the primes explicitly, so, for example, $\phi''''$ could just be written as $\phi^{[4]}$.
EDIT 2
I provide here the explicit form of first two $\mathcal{F}_n$ after $n=0$.
$$ \mathcal{F}_1 = \Box \phi - \partial (\partial \cdot \phi) + \partial^2 \phi'$$
$$ \mathcal{F}_2 = \Box \phi - \partial (\partial \cdot \phi) + \frac{2}{3 \Box} \partial^2 (\partial^2 \cdot \phi) + \frac{1}{3} \partial^2 \phi' - \frac{1}{\Box} \partial^3 (\partial \cdot \phi') + \frac{1}{\Box} \partial^4 \phi''$$
EDIT 3
Note that explicit indices are not required in this formulation, I only wrote the explicit forms to show where the notation comes from. I do not need to implement indices in the code.
