Astropy `SkyCoord` gives incorrect results when transforming points near the Declination pole

Question

The Astropy SkyCoord class provides incorrect values of a coordinate that is located near, or at, a coordinate pole when using its internal functions to transform between coordinate systems.

For example: start with a point that is near the Declination pole in ICRS coordinates, transform it to galactic coordinates, then transform it back. When the point approaches (or is set to be at) the Declination pole, the transformations yield a Right Ascension value that is wrong by several degrees.

Minimal example: start with a point offset from the north Declination pole by dec_delta (originally at Dec=80 degrees), slowly decrease dec_delta while converting to galactic coordinates and back, and see the resulting Right Ascension value deviate more from the correct value (300 degrees).

>>> from astropy.coordinates import SkyCoord
>>> ra_test = 300; dec_test = 90

>>> dec_delta = 1
>>> icrs_test = SkyCoord(ra_test, dec_test - dec_delta, unit='deg', frame='icrs'); icrs_test
<SkyCoord (ICRS): (ra, dec) in deg
    (300., 89.)>
>>> icrs_test.galactic.icrs    # convert to galactic coordinates and back to ICRS
<SkyCoord (ICRS): (ra, dec) in deg
    (300., 89.)>
# correct

>>> dec_delta = 1E-10
>>> icrs_test = SkyCoord(ra_test, dec_test - dec_delta, unit='deg', frame='icrs'); icrs_test
<SkyCoord (ICRS): (ra, dec) in deg
    (300., 90.)>
>>> icrs_test.galactic.icrs
<SkyCoord (ICRS): (ra, dec) in deg
    (299.99725044, 90.)>
# slightly wrong

>>> dec_delta = 1E-15
>>> icrs_test = SkyCoord(ra_test, dec_test - dec_delta, unit='deg', frame='icrs'); icrs_test
<SkyCoord (ICRS): (ra, dec) in deg
    (300., 90.)>
>>> icrs_test.galactic.icrs
<SkyCoord (ICRS): (ra, dec) in deg
    (291.41843232, 90.)>
# very wrong

>>> dec_delta = 0
>>> icrs_test = SkyCoord(ra_test, dec_test - dec_delta, unit='deg', frame='icrs'); icrs_test
<SkyCoord (ICRS): (ra, dec) in deg
    (300., 90.)>
>>> icrs_test.galactic.icrs
<SkyCoord (ICRS): (ra, dec) in deg
    (291.41843232, 90.)>
# problem remains at the pole

This doesn't seem to be a problem when near/at the Right Ascension pole:

>>> ra_test = 360
>>> ra_delta = 1E-15
>>> icrs_test = SkyCoord(ra_test - ra_delta, dec_test, unit='deg', frame='icrs'); icrs_test
<SkyCoord (ICRS): (ra, dec) in deg
    (0., 80.)>
>>> icrs_test.galactic.icrs
<SkyCoord (ICRS): (ra, dec) in deg
    (360., 80.)>
# correct since 0=360 for periodic RA

>>> ra_delta = 0
>>> icrs_test = SkyCoord(ra_test, dec_test, unit='deg', frame='icrs'); icrs_test
<SkyCoord (ICRS): (ra, dec) in deg
    (0., 80.)>
>>> icrs_test.galactic.icrs
<SkyCoord (ICRS): (ra, dec) in deg
    (360., 80.)>

The issue persists when using the underlying function that SkyCoord.galactic (or SkyCoord.icrs) use (SkyCoord.transform_to(...)). It is present for astropy versions 4, 5, and 6. I haven't tested other coordinate systems yet since I don't use them in my code.

Is this intended? Is there a way to properly transform coordinates near/at the Declination pole?