Firstly excuse any sloppiness here -- I'm not a mathematician by training so I've had a difficult time formalizing my question and tracking down relevant material.
Consider a point in a smooth manifold, $p \in M$. Linear approximations to functions passing through that point span the first-order jet space, $J_p$ which can be interpreted as linear maps over the tangent space at that point, $T_p M$. Compositions of these functions manifests in a straightforward compositional algebra of jets where in local coordinates the corresponding Jacobians are multiplied together.
Multiplying a compositional Jacobian by an element of the tangent space gives directional derivatives. This action is exactly the action of first order forward mode automatic differentiation. We can also take the adjoint of this compositional Jacobian and multiply it by an element of the cotangent space to give gradients. This is similarly the action of first order reverse mode automatic differentiation. Everything at first order is wonderful.
My question concerns what happens for higher-order approximations to functions and the compositional structure of the corresponding higher-order jets.
At higher-orders the partial derivatives of each order will in general mix whenever two jets are composed together, and then this composite object acts on a direct product of higher-order tangent bundles. For example, at second-order the composite map eats two elements of the first-order tangent space, $u_1, u_2 \in T_p M$ and one element of $v \in T^{2}_p M$ and returns two first-order directional derivatives, $J \cdot u_1$ and $J \cdot u_2$, and one second-order directional derivative, $u_1^{T} \cdot H \cdot u_2 + J v$, where $H$ is the Hessian. There is clearly some interesting structure here where the image of the second-order jets seems to decompose into $T_p M \otimes T_p M \otimes T^{2}_p M$ but my math isn't good enough to work out the general theory or find the right references. Does anyone know what kind of algebraic structure is at work here, either for the jets or the action of the jets, or can suggest appropriate references?
The immediate follow up is then how does this structure admit adjoints? If the above gives higher-order forward mode automatic differentiation, how would we be able to take the adjoint of various subspaces of jets (or jet tangent products) to enable higher-order reverse mode automatic differentiation?
Using higher-order dual numbers I worked out the second and third-order behavior heuristically in Chapter 1 of https://github.com/stan-dev/nomad/tree/master/manual (again please excuse the mathematical sloppiness or poor notation) but I would love to have a better geometric/algebraic structure of what's going on.
Thanks!