Não expanda subexpressões até que seja necessário

Isto não é uma resposta, porque não sei o que você quer dizer. No entanto, entendo o seu problema. Você sabe que thetaé uma função de aX, aY, aZ, assim como os componentes de incQuat. E se trabalharmos com funções aplicadas indefinidas, como esta?

from sympy import *

aX, aY, aZ = symbols('a_X a_Y a_Z')
rotInc = Matrix(3,1,[aX, aY, aZ])
theta = Function("theta")(aX, aY, aZ)
incQuat = Quaternion.from_axis_angle(rotInc/theta, theta*2)

# substitution dictionary, will be used later
subs_d1 = {
    theta: sqrt((rotInc.T @ rotInc)[0,0])
}
symb_components = [Function(f"f_{i}")(aX, aY, aZ) for i in range(4)]
# substitution dictionary, will be used later
subs_d2 = {f: c for f, c in zip(symb_components, incQuat.to_Matrix())}
incQuat_symbolic = Quaternion(*symb_components)

qX, qY, qZ, qW = symbols('q_X q_Y q_Z q_W')
baseQuat = Quaternion(qW, qX, qY, qZ)
poseQuat = incQuat_symbolic * baseQuat

d4 = diff(poseQuat, aX)
d4

Esta é a expressão do seu primeiro passo. Então, quando você se sentir pronto para realizar o cálculo propriamente dito, substitua nos dicionários:

d4.subs(subs_d2).doit().subs(subs_d1).doit()

e você acabará com expressões muito longas...

Achei que isso estivesse em algum lugar do SymPy, mas talvez o uso de cse na preparação de funções compiladas (ou algo assim) seja onde ele é usado. Mas, seguindo sua ideia (e provavelmente redundante em relação ao que você pensa nesta questão ):

from sympy import symbols, Function, Derivative, cse

def dep(r, x):
    reps = {}
    for v, e in r:
        if e.xreplace(reps).has(x):
            reps[v] = Function(str(v))(x)
    return reps

def differentiate_with_cse(expr, x, backsub=False):
    r, e = cse(expr)
    reps = dep(r, x)
    dr = []
    for v, s in r:
        if v in reps:
            f = reps[v]
            df = Derivative(f, x).doit()
            ds = s.xreplace(reps).diff(x)
            dr.append((df, ds))
    dexpr = e[0].subs(reps).diff(x).subs(dr).expand()
    for k, v in reversed(reps.items()):
        dexpr = dexpr.subs(v, k)
    return dexpr.subs(list(reversed(r))) if backsub else dexpr,dict(r)

...your code up to derivative

>>> deriv, pats = differentiate_with_cse(poseQuat, aX)
>>> deriv.simplify()
(-qX*x7) + qW*x7*i + (-qZ*x7)*j + qY*x7*k
>>> pats
{x0: aX**2, x1: aY**2, x10: aZ*x7, x2: aZ**2, x3: x0 + x1 + x2, x4: sqrt(x3), x5: cos(x4), x6: 1/x3, x7: sin(x4)/(x4*sqrt(x0*x6 + x1*x6 + x2*x6)), x8: aX*x7, x9: aY*x7})

Então, com a dica das funções do @Davide_sd, criei um método genérico que me permite controlar facilmente como as subetapas são divididas. Basicamente, estou derivando manualmente as funções que separo, mas, assim como o CSE, mantenho os resultados em um dicionário para compartilhar entre todas as ocorrências.
As expressões base que compõem o cálculo são inseridas e nunca modificadas, a lista de derivação é semeada com o que você deseja derivar (múltiplas expressões são aceitáveis) e ele as derivará recursivamente, usando a lista de expressões conforme necessário.
No final, ainda posso usar o CSE para a) trazê-lo para esse formato, caso você precise, e b) fatorar ocorrências ainda mais comuns.
Funciona razoavelmente bem com meu pequeno exemplo; posso atualizá-lo conforme eu adicionar mais complexidade à função que preciso derivar.

from sympy import *

def find_derivatives(expression):
    derivatives = []
    if isinstance(expression, Derivative):
        #print(expression)
        derivatives.append(expression)
    elif isinstance(expression, Basic):
        for a in expression.args:
            derivatives += find_derivatives(a)
    elif isinstance(expression, MatrixBase):
        for i in range(rows):
            for j in range(cols):
                derivatives += find_derivatives(self[i, j])
    return derivatives

def derive_recursively(expression_list, derive_done, derive_todo):
    newly_derived = {}
    for s, e in derive_todo.items():
        print("Handling derivatives in " + str(e))
        derivatives = find_derivatives(e)
        for d in derivatives:
            if d in newly_derived:
                #print("Found derivative " + str(d) + " in done list, already handled!")
                continue
            if d in derive_todo:
                #print("Found derivative " + str(d) + " in todo list, already handling!")
                continue
            if d in expression_list:
                #print("Found derivative " + str(d) + " in past list, already handled!")
                continue
            if d.expr in expression_list:
                expression = expression_list[d.expr]
                print("  Deriving " + str(d.expr) + " w.r.t. " + str(d.variables))
                print("    Expression: " + str(expression))
                derivative = Derivative(expression, *d.variable_count).doit().simplify()
                print("    Derivative: " + str(derivative))
                if derivative == 0:
                    e = e.subs(d, 0)
                    derive_todo[s] = e
                    print("        Replacing main expression with: " + str(e))
                    continue
                newly_derived[d] = derivative
                continue
            print("Did NOT find base expression " + str(d.expr) + " in provided expression list!")
    derive_done |= derive_todo
    if len(newly_derived) == 0:
        return derive_done
    return derive_recursively(expression_list, derive_done, newly_derived)

incRot_c = symbols('aX aY aZ')
incRot_s = Matrix(3,1,incRot_c)
theta_s = Function("theta")(*incRot_c)
theta_e = sqrt((incRot_s.T @ incRot_s)[0,0])
incQuat_c = [ Function(f"i{i}")(*incRot_c) for i in "WXYZ" ]
incQuat_s = Quaternion(*incQuat_c)
incQuat_e = Quaternion.from_axis_angle(incRot_s/theta_s, theta_s*2)
baseQuat_c = symbols('qX qY qZ qW')
baseQuat_s = Quaternion(*baseQuat_c)
poseQuat_c = [ Function(f"p{i}")(*incRot_c, *baseQuat_c) for i in "WXYZ" ]
poseQuat_s = Quaternion(*poseQuat_c)
# Could also do it like this and in expressions just refer poseQuat_s to poseQuat_e, but output is less readable
#poseQuat_s = Function(f"pq")(*incRot_c, *baseQuat_c)
poseQuat_e = incQuat_s * baseQuat_s

expressions = { theta_s: theta_e } | \
    { incQuat_c[i]: incQuat_e.to_Matrix()[i] for i in range(4) } | \
    { poseQuat_c[i]: poseQuat_e.to_Matrix()[i] for i in range(4) }

derivatives = derive_recursively(expressions, {}, { symbols('res'): diff(poseQuat_s, incRot_c[0]) })
print(derivatives)

elements = cse(list(expressions.values()) + list(derivatives.values()))
pprint(elements)