UnitQuaternions and ForwardDiff

When differentiating a function that takes a quaternion as an input, I get several different solutions for different setups.
I know it's not really correct to treat a quaternion as a regular 4D vector when differentiating, but I'm a bit confused about the dimensions and entries of some of the results.

Code:

using LinearAlgebra
using ForwardDiff
using Rotations

function testfoo1(q)
    q = UnitQuaternion(q) 
    Rotations.params(q)
end

function testfoo2(q)
    q
end

q1 = rand(4)
q1 = q1/norm(q1)
q2 = UnitQuaternion(q1)

ForwardDiff.jacobian(testfoo1,q1) # 4x4 "random" matrix
ForwardDiff.jacobian(testfoo1,q2) # 4x9 "random" matrix
ForwardDiff.jacobian(testfoo2,q1) # 4x4 identity matrix
ForwardDiff.jacobian(testfoo2,q2) # 9x9 identity matrix

What I'd expect (but that's definitely debatable):

ForwardDiff.jacobian(testfoo1,q1) # 4x4 identity matrix
ForwardDiff.jacobian(testfoo1,q2) # 4x4 identity matrix or 4x3 "proper" quaternion derivative
ForwardDiff.jacobian(testfoo2,q1) # 4x4 identity matrix
ForwardDiff.jacobian(testfoo2,q2) # 4x4 identity matrix or 4x3 "proper" quaternion derivative

Especially the first one ( ForwardDiff.jacobian(testfoo1,q1) ) where I pass in a regular vector and return a regular vector seems off to me. And dimension 9 for some of the matrices also seems off; does that come from writing the quaternions as rotation matrices?

Start by reproducing the four ForwardDiff.jacobian calls using testfoo1, testfoo2, q1, and q2 from the issue. Read the UnitQuaternion and Rotations.params entry points alongside ForwardDiff's handling of these values; done means explaining the differing dimensions and whether the observed Jacobians are expected or require a targeted change.