A question on fractional derivatives

There are basically no interesting solutions to this equation beyond first and zeroth order operators, even if one only imposes the stated constraint for $n=2$.

First, we can depolarise the hypothesis $$ D^u(f^2) = \alpha_2 D^u(f) f \quad (1)$$ by replacing $f$ with $f+g, f-g$ for arbitrary functions $f,g$ and subtracting (and then dividing by $4$) to obtain the more flexible Leibniz type identity $$ D^u(fg) = \frac{\alpha_2}{2}( D^u(f) g + f D^u(g) ). \quad (2)$$

There are now three cases, depending on the value of $\alpha_2$:

1. $\alpha_2 \neq 1,2$. Applying (2) with $f=g=1$ we then conclude that $D^u(1)=0$, and then applying (2) again with just $g=1$ we get $D^u(f)=0$. So we have the trivial solution $D^u=0$ in this case.
2. $\alpha_2=2$. Then $D^u$ is a derivation and by induction we have $D^u(f^n) = n D^u(f) f^{n-1}$, just as with the ordinary derivative, so we just have $\alpha_n=n$ for all $n$ with no fractional behavior.
3. $\alpha_2=1$. Applying (2) with $g=1$ we obtain (after a little bit of algebra) $D^u(f) = mf$ where $m := D^u(1)$. Thus $D^u$ is just a multiplier operator, which obeys $D^u(f^n) = D^u(f) f^{n-1}$, thus $\alpha_n=1$ for all $n$.

Thus there are no linear solutions to your equation other than the usual derivations (e.g., $D^u(f) = a(x) \frac{d}{dx} f$ for any smooth symbol $a$) and multiplier operators $D^u(f) = mf$, i.e., first order and zeroth order operators.

On the other hand, fractional derivatives $D^u$ tend to obey a "fractional chain rule" $$ D^u( F(f) ) = D^u(f) F'(f) + E$$ for various smooth functions $F,f$, where the error $E$ obeys better estimates in various Sobolev spaces than the other two terms in this equation. In particular, for $F(t) = t^n$, we would have $$ D^u(f^n) = n D^u(f) f^{n-1} + E$$ for a "good" error term $E$. For instance, taking $u=n=2$ with $D$ the usual derivative, we have $$ D^2(f^2) = 2 D^2(f) f + E \quad (3)$$ with $E$ the "carré du champ" operator $$ E := 2 (Df)^2.$$ Note that the error $E$ is controlled uniformly by the $C^1$ norm of $f$ but the other two terms in (3) are not. See my previous MathOverflow answer at https://mathoverflow.net/a/94039/766 for some references and further discussion.

It appears you actually want $D^u(f^n)=\alpha f^{n-1} D^u f$, where $\alpha$ is a scalar.

There is no reason for this to be true, and this is indeed false in general. E.g., for $n=2$ and the Riemann--Liouville fractional derivative of $f:=\exp$ with $u=1/2$, $a=0$, and $x>0$ we have $$f(x)^{n-1}(D^uf)(x)=e^{2 x} \text{erf}\left(\sqrt{x}\right)+\frac{e^x}{\sqrt{\pi } \sqrt{x}},$$ whereas $$(D^u(f^n))(x)=\sqrt{2} e^{2 x} \text{erf}\left(\sqrt{2} \sqrt{x}\right)+\frac{1}{\sqrt{\pi } \sqrt{x}},$$ so that $$\frac{D^u(f^n)}{f^{n-1}\,D^uf}$$ is quite unlike any constant.

Moreover, the term $\text{erf}\left(\sqrt{2} \sqrt{x}\right)$ in the expression for $(D^u(f^n))(x)$ here versus the term $\text{erf}\left(\sqrt{x}\right)$ in the expression for $f(x)^{n-1}(D^uf)(x)$ seem to make it very unlikely that any other kind of fractional derivative will work as you want.

The generalized Leibniz formula applicable to the classic fractional integroderivative is

$$ D^{\omega}\; f(x)g(x) = \sum_{n \geq 0} \binom{\omega}{n} [D^{\omega-n}f(x)]D^ng(x)=(D_L+D_R)^{\omega} g(x)f(x),$$

where $D_L$ acts on the function on the left of the product and $D_R$ on the right function. See, e.g., Leibniz rules and integral analogues for fractional derivatives via a new transformation formula by Fugere, Gaboury, and Tremblay.

This generalized Leibniz rule applies to the fractional integroderivative satisfying the sensible axioms given by Pincherle described in "The Role of Salvatore Pincherle in the Development of Fractional Calculus" by Francesco Mainardi and Gianni Pagnini--those satisfied by the usual derivative raised to integral powers, negative or positive. Reps of this op are presented in this MSE-Q and can be used to define the confluent (see this MO-Q) and regular hypergeometric fcts.

These reps of $D^{\omega}$ are at the heart of the definitions of the Euler gamma and beta functions via integrals, generalizations of the integral factorials and integral binomial coefficients (see my reply to/refs in this MO-Q), which most researchers use frequently in their math endeavors--contrary to some opinions expressed on MO. See an example of the half-derivative in this MO-Q (which many users apparently confuse with some pseudo-differential operator defined by the Fourier transform).