ode DormandPrince solver second derivative support

When using the DormandPrince solver for ode, the second derivative is always nan. Is it a bug or is it not supported?
The BDF solver supports the second derivative.

tensorflow version: 2.8.0
tensorflow-probability: 0.16.0

Example:

import tensorflow as tf
import tensorflow_probability as tfp

def f(t, x):
    return 2 * tf.expand_dims(t, axis=0)

x_init = tf.ones((1,))

t_list = [0.0, 1.0]
with tf.GradientTape(persistent=True) as tape:
    tape.watch(x_init)
    #solver = tfp.math.ode.BDF()
    solver = tfp.math.ode.DormandPrince(atol=1E-5)
    times, states, _, _ = solver.solve(f, initial_time=t_list[0], initial_state=x_init, solution_times=t_list)
    diff = tape.gradient(states[1], x_init)
    print(diff) # is 1.0
diff = tape.gradient(diff, x_init)
print(diff) # is Nan, it shall be 0.0
# tf.Tensor([1.], shape=(1,), dtype=float32)
# tf.Tensor([nan], shape=(1,), dtype=float32)

Start by running the supplied TensorFlow 2.8.0 and TensorFlow Probability 0.16.0 example, comparing DormandPrince with BDF and checking the reported first- and second-derivative values. Trace the DormandPrince solver path to determine whether higher-order gradients are supported; done means the second derivative behaves as expected or the limitation is clearly documented.