Combining SHO and Rotation terms derivative issue

I've been having an issue combining SHO and Rotation terms for describing stellar light curves. I've tried on many light curves, in many environments, devices, and with the PyMC3 and PyMCv5 versions of exoplanet. Troubleshooting has led me to believe the issue is in the the derivative of the log probability when both terms are combined; the derivate is very "spiky", leading to the scipy optimization step (pmx.optimize) either failing to improve log p due to precision loss or jumping the GP parameters into unreasonable values given their priors. Any NUTS sampling also fails, as the step size is too small and max tree depth is reached with every step.

To further try to troubleshoot the issue, I just ran the tutorial workbook Gaussian Process Models for Stellar Variability (https://gallery.exoplanet.codes/tutorials/stellar-variability/) which incorporates a kernel that is the sum of an SHO and Rotation term. Running that notebook after in an environment with exoplanet==0.5.3 and celerite2==0.3.2 yields the following outputs.

From defining the model, and running pmx.optimize():

From attempting to sample:

The convergence checks and posterior:

I noticed that the published version of the notebook uses exoplanet version 0.5.4.dev27+g75d7fcc, so I created a cloned environment, upgraded exoplanet to that version, and found a commit of celerite2 pushed on the same day (celerite2-0.3.2.dev1+g50193ff)...and it worked!

Result from pmx.optimize():

Result of sampling:

Convergence check and posterior:

So it seems something since that version of celerite2 has broken the ability to add those two terms. Using either term alone seems to work fine, but not in combination.