I'm not aware of any Arduino library or ready-to-use code that does quite what you want. Maybe it's out there? But my guess is you're going to have to do some substantial programming work to build this.

At least you *can* get the micros() count pretty easily, and you can (with some not-very Arduino-like syntax) raise your interrupt to top priority (if that's even necessary) to really get a PPS capture timestamped to within 1 us.

Maybe you'll just subtract each one from the previous, to get the number of microseconds of "Teensy time" which actually elapsed using 1 GPS second? Or maybe you'll average the last several of those, or maybe run them through a low-pass filter or other algorithm to give you finer resolution and reject the likely +/- 1us measurement jitter?

How you're going to actually use that info to process your MPU9250 data is a good question. If you're already doing some sort of filter based on the assumed sample rate (especially of the gyro data), maybe the difference in actual versus assumed/nominal sample rate can be applied as merely a minor (ideally linear) correction factor if the difference between Teensy's crystal and the GPS true PPS time is very small?

You might also consider whether any of this is even meaningful. MEMS-based motion measurement typically has analog bandwidth in the dozens to (maybe) hundreds of Hz. Some chips allow you to select different filters. But internally, they all work on the same principle, where a tuning fork oscillators in the dozens to hundreds of kHz and motion modulates it. The changes are demodulated to a pulsing signal and then filtered. The final result necessarily has fairly low bandwidth to give you high resolution, or more noise if you want higher speed.

We've actually had several other people ask here about doing all sorts of elaborate things with MEMS motion sensors. I recall one recent thread involved trying to sample hundreds of (unspecified) on/off sensors (of unspecified bandwidth) at some extremely high speed and then trying to correlate that data with an IMU chip. It's easy to imagine building such incredibly high speed signal acquisition, without really considering whether the analog signal sources even has high enough speed/bandwidth and reproducible enough timing characteristics to make all that effort worthwhile.

Most of these sensors also have offsets & drift. That's why the filters are used, especially for accelerometers and magnetometers to correct for gyroscope drift. I'm having a hard time imagining how very minor timing corrections could end up being anywhere on the scale of the pretty substantial drift in these measurements.

Anyway, this question does come up from time to time. If you do find a want to use this data and somehow eek out a tiny bit extra performance that actually matters with real IMU data, I hope you'll consider sharing anything you learn along the way?