Ah, sorry… I wrote my own and only skimmed yours.
(∂1[1] * da + ∂2[1] * dx,∂1[2] * da + ∂2[2] * dx)
is not correct. ∂1 already multiplies with da internally: ∂Γ_∂a(_a,x,ind)*partials(a).
The ChainRules rules should go in SpecialFunctions I suppose, which then will need to depend on HypergeometricFunctions… Don’t know if the maintainers will go for that.
Same with ForwardDiff. They define dgamma_inc/dx themselves, and that could be extended, but again would necessitate another dependence.
There is also ForwardDiffChainRules to translate ChainRules frules to ForwardDiff.
I’m not a contributer to either ecosystem though, so I’d say talk to the maintainers before authoring any PRs 