Compute vernal equinox in Astropy

Astropy
1

I need to compute the date and time of the vernal equinox for a given year. I have the following code:

from astropy.time import Time
from astropy.coordinates import get_sun, AltAz, EarthLocation
import numpy as np
import astropy.units as u

def find_vernal_equinox(year):
    # Make an array of times around the expected equinox (March 19-21)
    times = Time(f"{year}-03-19 00:00:00") + np.linspace(0, 3, 1000) * u.day
    
    # Get the declination of the Sun at those times
    sun_coords = get_sun(times)
    declinations = sun_coords.dec.deg

    # The equinox occurs when declination crosses 0
    zero_crossings = np.where(np.diff(np.sign(declinations)))[0]
    if not len(zero_crossings):
        raise ValueError("No zero crossing found!")

    idx = zero_crossings[0]
    # Interpolate to get a more precise time
    t1, t2 = times[idx], times[idx+1]
    d1, d2 = declinations[idx], declinations[idx+1]
    # Linear interpolation
    frac = -d1 / (d2 - d1)
    equinox_time = t1 + frac * (t2 - t1)
    return equinox_time.utc.iso

# Example usage:
print(find_vernal_equinox(2025))

Strangely, the results I get are about 8 or 9 hours off (all times UTC, actual data from https://data.giss.nasa.gov/modelE/ar5plots/srvernal.html):
Computed 2025-03-20 17:34:03.719 – Actual 2025 3/20 9:00
Computed 2026-03-20 23:39:44.858 – Actual 2026 3/20 14:49
Computed 2027-03-21 05:41:50.943 – Actual 2027 3/20 20:38

Moreover, the declination of the Sun at the time of the equinox is not zero:

print(get_sun(Time(f"2025-03-20 09:00:00")).dec.deg)
# returns -0.14103777488837554

What am I doing wrong?