Quaternion Converted to Rotation Matrix then Derived with Respect to this Quaternion

I was wondering how a derivative of a rotation matrix generated based on a quaternion and then differentiated with respect to this quaternion would be calculated.

$$q_{1\times 4} = [q_0 q_1 q_2 q_3]^T$$

$$d_{rot_{3\times3}} = \frac{\delta C(q)}{\delta q}$$

If I attempt to lay out the rotation matrix by its elements and derive it with respect to each quaternion element, I would need to differentiate a matrix w.r.t. a vector which results in a 3D tensor differentiation see here. But If it is possible to derive it from the symbolic expression it eventually should remain 2D.

As a second case: If I rotate a vector with this rotation and differentiate it afterward, then I get a 2D Jacobian.

$$d_{rot_{3\times3}} = \frac{\delta \left(C(q_1) \cdot v_{3x1}\right)}{\delta q}$$

I'm trying to remain in 2D if possible. Also, rotation quaternions do have norm = 1 if this helps to simplify the differentiation.

Furthermore, I'm also interested in the Rank and if this constellation always has full Rank, then I would need to prove it.

J. Kelly "Indirect Kalman filter for 3D Attitude Estimation"

Lemma 1 Shows that based on the unit quaternion fact, the result in his example should have full Rank. But this only accounts for a quaternion and not for a rotation matrix based on quaternions.

I am recently encountering a similar problem as yours. I managed to show (at least partly) that $$\dfrac{d}{dt}R(\mathfrak{q})=R(\mathfrak{q})S(\omega),$$ where $$R:\mathbb{S}^3\to\mathsf{SO}(3)$$ is defined by $$R(\mathfrak{q}):=I_3+2\mathfrak{q}_0 S(\mathfrak{e})+2S(\mathfrak{e})^2$$ for each $$\mathfrak{q}:=(\mathfrak{q}_0,\mathfrak{e}):=(\mathfrak{q}_0,\mathfrak{q}_1,\mathfrak{q}_2,\mathfrak{q}_3)\in\mathbb{S}^3$$, where $$S$$ is such that $$S(v)w=v\times w$$ for each $$v,w\in\mathbb{R}^3$$.
The way is to consider a partition of $$R$$ such that $$R(\mathfrak{q})=[R(\mathfrak{q})e_1\,R(\mathfrak{q})e_2\,R(\mathfrak{q})e_3]$$, then consider for each $$k\in\{1,2,3\}$$ that $$\dfrac{d}{dt}R(\mathfrak{q})e_k=\dfrac{\partial R(\mathfrak{q})e_k}{\partial \mathfrak{q}}\dot{\mathfrak{q}},$$ in which the matrix $$\dfrac{\partial R(\mathfrak{q})e_k}{\partial \mathfrak{q}}=\left[\dfrac{(\partial R(\mathfrak{q})e_k)_i}{\partial \mathfrak{q}_j}\right]_{ij}$$ and $$\dot{\mathfrak{q}}=\frac{1}{2}\begin{pmatrix}-\mathfrak{e}^\top\\\mathfrak{q}_0I_3+S(\mathfrak{e})\end{pmatrix}\omega.$$ I checked that for $$k=1$$, the above leads to $$R(\mathfrak{q})S(\omega)e_1$$.
Although this process does not answer your question, I believe it contains the elements you want, namely the matrix $$\dfrac{\partial R(\mathfrak{q})e_k}{\partial \mathfrak{q}}$$. Note that the quaternion here is in column vector form.