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

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))

I have tried my best to avoid my feel of possibility of opinion based answers on this question and looked this up on internet. Found out it is quite hard to find any relevant info that I could trust, so I enumerate these posts (except wikipedia):

timeanddate is stating the exact time is 18:20 UTC on December 21, which is as you have calculated

winstars have stated the time when the planets will be at closest angle as 18:25 UTC and they mention that conjunction will occur at 13:30 UTC, I am not sure if that is the first time.

Not sure how relevant is this, but the conjunction here is stated to be 6.2 degreest at 17:32 GMT, thus 18:32 UTC

The most relevant source I was able to find was in the sky, where the time was estimated exactly to 13:24 UTC., based on calculations on data by Jet Propulsion Laboratory - source code can be checked here (c).

You can see that mostly not both types of calculation are used, and that the times vary. The reason of that is that in such calculations you need very long floats for best precision. As you are limited by the machine you use, the precision is not perfect. Such as @bad_coder has suggested, you might get better answer in Astronomy stack exchange.