Discrepancy between my results of Great Conjunction 2020 and Wikipedia Using python-skyfield to calculate the upcoming conjunction of Jupiter and Saturn. I'm off a few minutes with the times on Wikipedia.
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) &gt; 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) &gt; 0


longitude_difference.rough_period = 300.0
longitude_difference1.rough_period = 300.0

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

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

Can anyone verify the correct UTC time?  Is there a bug in my code?
There are no bugs in your codes. But you can make it more accurate with the following corrections (additions):

Using new ephemeris series (de430 or de440) most recent one is de440.
Use True equinox of True epoch of the date. Adding epoch='date' as below:

.radec(epoch='date')
.ecliptic_latlon(epoch='date')


use the following codes for loading timescale instead of ts = load.timescale(builtin=True).

url = load.build_url('finals2000A.all')
with load.open(url) as f:
    finals_data = iers.parse_x_y_dut1_from_finals_all(f)

ts = load.timescale()
iers.install_polar_motion_table(ts, finals_data)

These codes load the updated values of ΔT and also apply polar motion. Of course, because your calculations are geocentric, the polar motion does not affect the result.
You can reach the following result by the corrections and some extra calculations for other attributes.
JUPITER -- SATURN Conjunction 2020:

A.D. 2020-Dec-21 18:20:37.3129 UT
Angular Separation: +  0° 06ʹ 06.38154ʺ
JUPITER is -  0° 06ʹ 06.98608ʺ of ('Below', 'Right') of SATURN in direction of (+  0° 00ʹ 02.56845ʺ, N)