My Jupiter/Saturn conjunction calculation with skyfield doesn't match wikipedia

Using Python-Skyfield to calculate the upcoming conjunction if Jupiter and Saturn.

Wikipedia Great conjunction times (1800 to 2100)

Using Right Ascension:

Dec 21,2020 13:22:00 UTC - Wikipedia.
Dec 21,2020 13:34:33 UTC - My Calculation.

Using Ecliptic Longitude:

Dec 21,2020 18:37:31 UTC - Wikipedia
Dec 21,2020 18:20:40 UTC - My Calculation.

from skyfield.api import load, tau, pi
from skyfield.almanac import find_discrete

planets = load('de421.bsp')
sun = planets['sun']
earth = planets['earth']
jupiter = planets['jupiter barycenter']
saturn = planets['saturn barycenter']

ts = load.timescale(builtin=True)

def longitude_difference(t):
    e = earth.at(t)
    j = e.observe(jupiter).apparent()
    s = e.observe(saturn).apparent()
    _, lon1, _ = s.ecliptic_latlon()
    _, lon2, _ = j.ecliptic_latlon()
    return (lon1.degrees - lon2.degrees) > 0

def longitude_difference1(t):
    e = earth.at(t)
    j = e.observe(jupiter).apparent()
    s = e.observe(saturn).apparent()

    jRa, _, _ = j.radec()
    sRa, _, _ = s.radec()
    return (sRa._degrees - jRa._degrees) > 0


longitude_difference.rough_period = 300.0
longitude_difference1.rough_period = 300.0

print()
print("Great conjunction in ecliptic longitude:")
t, b = find_discrete(ts.utc(2020), ts.utc(2021), longitude_difference)
for ti in t:
    print(t.utc_jpl())

print()
print("Great conjunction in right ascension:")
t, b = find_discrete(ts.utc(2020), ts.utc(2021), longitude_difference1)
for ti in t:
    print(t.utc_jpl())

I'm new to Skyfield so any help is appreciated.

Solution

Your conjection-finding code looks very good, and I suspect its results are far better than those of the Wikipedia — looking at the version history, it’s not clear where their numbers even came from, and unattributed astronomy calculations can’t easily be double-checked without knowing from which ephemeris and software they derived them.

I attach below a slightly improved version of your solver. Here are the tweaks I recommend:

Pass epoch='date' when computing coordinates, since conjunctions and oppositions are traditionally defined according to the celestial sphere of the year in which they happen, not the standard J2000 celestial sphere.
Before doing any math on the dates, force them into a range of zero to a full circle (360 degrees or 24 hours). Otherwise, you will see the result of the subtraction jump by ±360°/24h whenever one of the right ascensions or longitudes happens to cross 0°/0h and change signs. This will give you “phantom conjunctions” where the planets are not really swapping places but merely switching the sign of the angle they’re returning. (For example: jumping from -69° to 291° is really no motion in the sky at all.)
Note that both of our scripts should also find when Jupiter and Saturn are across the sky from each other, since the sign of their difference should also flip at that point.
In case any sources you track down cite the moment of closest approach or the angle between the planets at that moment, I've added it in.

Here’s the output, very close to yours, and again not agreeing well with those old unexplained uncredited Wikipedia numbers:

Great conjunction in ecliptic longitude:
['A.D. 2020-Dec-21 18:20:37.5144 UT']

Great conjunction in right ascension:
['A.D. 2020-Dec-21 13:32:04.9486 UT']

Great conjunction at closest approach:
A.D. 2020-Dec-21 18:21:00.2161 UT - 0.1018 degrees

And the script:

from skyfield.api import load
from skyfield.searchlib import find_discrete, find_minima

planets = load('de421.bsp')
sun = planets['sun']
earth = planets['earth']
jupiter = planets['jupiter barycenter']
saturn = planets['saturn barycenter']

ts = load.timescale(builtin=True)

def longitude_difference(t):
    e = earth.at(t)
    j = e.observe(jupiter).apparent()
    s = e.observe(saturn).apparent()
    _, lon1, _ = s.ecliptic_latlon(epoch='date')
    _, lon2, _ = j.ecliptic_latlon(epoch='date')
    return (lon1.degrees - lon2.degrees) % 360.0 > 180.0

def ra_difference(t):
    e = earth.at(t)
    j = e.observe(jupiter).apparent()
    s = e.observe(saturn).apparent()

    jRa, _, _ = j.radec(epoch='date')
    sRa, _, _ = s.radec(epoch='date')
    return (sRa.hours - jRa.hours) % 24.0 > 12.0

def separation(t):
    e = earth.at(t)
    j = e.observe(jupiter).apparent()
    s = e.observe(saturn).apparent()

    return j.separation_from(s).degrees

longitude_difference.step_days = 30.0
ra_difference.step_days = 30.0
separation.step_days = 30.0

print()
print("Great conjunction in ecliptic longitude:")
t, b = find_discrete(ts.utc(2020), ts.utc(2021), longitude_difference)
for ti in t:
    print(t.utc_jpl())

print()
print("Great conjunction in right ascension:")
t, b = find_discrete(ts.utc(2020), ts.utc(2021), ra_difference)
for ti in t:
    print(t.utc_jpl())

print()
print("Great conjunction at closest approach:")
t, b = find_minima(ts.utc(2020, 6), ts.utc(2021), separation)
for ti, bi in zip(t, b):
    print('{} - {:.4f} degrees'.format(ti.utc_jpl(), bi))