First, there is a fundamental error in your code. earth_sun_gcrs = get_sun(...) is the location of Sun center in the GCRS frame, so earth_sun_icrs = earth_sun_gcrs.transform_to(ICRS) is the location of Sun center in the ICRS frame. That is, the vector you return is merely the vector from the origin of ICRS (the solar-system barycenter) to Sun center, with no connection to the Earth at all.

Second, if what you want is the apparent Earth-Sun vector relative to ICRF axes, including the effects of stellar and planetary aberration, then all you need to do is to skip the transformation to ICRS. The origin of GCRS is Earth center, so earth_sun_gcrs.cartesian is the Earth-Sun vector relative to ICRF axes, including the above observer effects.

Alternatively, if what you want is instead the true/geometric Earth-Sun vector relative to ICRF axes, you shouldnâ€™t use get_sun() at all. Instead, you can call get_body_barycentric() to get the true location of a solar-system body in ICRS, e.g.:

You can use ICRF axes with any origin. Setting aside some details:

ICRS = ICRF axes with the origin at the solar-system barycenter

GCRS = ICRF axes with the origin at Earth center

HCRS = ICRF axes with the origin at Sun center

etc.

My understanding of MJ2000 is that itâ€™d be very similar to GCRS.

If you have the Sunâ€™s location as a SkyCoord, you can transform to Astropyâ€™s TEME, which has the origin at Earth center, so the vector will be the Earth-Sun. Iâ€™ll caution that there are multiple definitions of TEME, so itâ€™d be prudent for you to confirm whether Astropyâ€™s definition of TEME is what you need.

The code you have for getting the true/geometric Earth-Sun vector relative to ICRF axes - how might I modify it to use GCRS, to be Earth-centered? Like this?:

Again, if what you want is the true/geometric Earth-Sun vector relative to ICRF axes, you shouldnâ€™t use get_sun() at all. Note that as the difference of two positions, this vector is the same irrespective of the notion of the coordinate-frame origin, e.g., SSB for ICRS or Earth center for GCRS.

However, if what you want is the Earth-Sun vector relative to ICRF axes, but including the effects of stellar and planetary aberration for an observer at Earth center, then yes, your new code block gives you that.