What are the Time Zones on Mercury?
Partial answer about the retrograde motion only:
Yes, the Sun has yearly retrograde motion as seen from Mercury
I've written a simple script to generate the orbital motion of Mercury based on its eccentricity. I've used reduced units with $GM = a = 1$ and period $T = 2 \pi$.
Mercury's rotational period is 2/3 of its orbital period.
I've plotted the rotational angle of mercury minus the angle of the Sun as seen from Mercury's position. To address retrograde motion, it's not necessary to consider from where on Mercury you are looking,
There are tiny dips below zero every year.


import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

def deriv(X, t):
    x, v = X.reshape(2, -1)
    acc = - x * ((x**2).sum())**-1.5
    return np.hstack((v, acc))

halfpi, pi, twopi, fourpi, sixpi = [f*np.pi for f in (0.5, 1, 2, 4, 6)]
to_degs, to_rads = 180/pi, pi/180

# Mercury
# https://space.stackexchange.com/questions/43213/what-is-the-analytical-closed-form-solution-of-the-two-body-problem-to-verify-it/43237#43237
# https://solarsystem.nasa.gov/planets/mercury/by-the-numbers/


e = 0.20563593
a = 1
peri = (1 - e) * a
v_peri = np.sqrt(2/peri - 1/a)

X0 = np.array([peri, 0, 0, v_peri])
times = np.linspace(0, sixpi, 1001)
dt_years = (times[1] - times[0]) / twopi

answer, info = ODEint(deriv, X0, times, full_output=True)
xy = answer.T[:2]
x, y = xy
r = np.sqrt((xy**2).sum(axis=0))

theta_rot = np.mod((3/2) * times, twopi)
sun = np.zeros(2)[:, None] - xy
theta_sun = np.arctan2(sun[1], sun[0])

if True:
    plt.figure()
    plt.subplot(2, 1, 1)
    angle = theta_rot - theta_sun
    dangle = np.mod(angle[1:] - angle[:-1]+halfpi, twopi)-halfpi
    dangle_dt = dangle / dt_years
    angle = np.mod(angle+halfpi, twopi)-halfpi
    plt.plot(times/twopi, angle)
    plt.ylabel('angle (rad)')
    plt.subplot(2, 1, 2)
    plt.plot(times[1:]/twopi, dangle_dt)
    plt.plot(times[1:]/twopi, np.zeros_like(dangle_dt), '-k', lw=0.5)
    plt.ylabel('dangle/dt (rad/year)')
    plt.xlabel('years')
    plt.show()


if True:
    plt.figure()
    plt.subplot(3, 1, 1)
    plt.plot(x, y)
    plt.plot(x[:1], y[:1], 'ok', ms=2)
    plt.plot([0], [0], 'oy', ms=8)
    plt.gca().set_aspect('equal')
    plt.subplot(3, 1, 2)
    plt.plot(times/twopi, x)
    plt.plot(times/twopi, y)
    plt.plot(times/twopi, r, '-k')
    plt.subplot(3, 1, 3)
    plt.plot(times/twopi, theta_rot)
    plt.plot(times/twopi, theta_sun)
    plt.show()