Is the shortest day duration constant? Is the exact duration of the shortest day of a year (at the winter solstice in the northern hemisphere) always the same, for given latitude and longitude, or does it fluctuate? If it does, what is the range?
By looking at the link provided by James K, I was able to write a Python 3 script that calculates the duration of the longest night for a predefined location (Boston, MA) for 160 years. The script uses the PyEphem module. The longest night becomes shorter by ~72 msec per year, and its length also oscillates with the period of ~19 years.



Here's the script that I used:

import ephem
import matplotlib.pyplot as plt
import matplotlib
import numpy as np

location = ephem.Observer()
location.lat, location.lon = '42.3601', '-71.0589'

years = range(1900, 2060)
longest_nights = []
for year in years:
    location.date = ephem.next_winter_solstice('{}-12-12'.format(year))
    sunset = location.next_setting(ephem.Sun())
    sunrise = location.next_rising(ephem.Sun())
    if sunrise &lt; sunset: # Look at the night before!
        location.date -= 1       
        sunset = location.next_setting(ephem.Sun())
    longest_nights.append(sunrise - sunset)

longest_nights = np.array(longest_nights)
longest_nights = (longest_nights - longest_nights.mean()) * 24 * 3600 # In seconds

# Find the slope
xp = np.poly1d(np.polyfit(years, longest_nights, 1))
np_years = np.array(years).reshape(len(years), 1)

# Plot
matplotlib.style.use("ggplot")
plt.plot(years, longest_nights, "-o")
plt.ylabel("Deviation from the mean, sec")
plt.xlabel("Year")
plt.title("Duration of the longest night in a year in Boston")
plt.plot(years, xp(years), "-")
plt.arrow(2017, plt.ylim()[0], 0, plt.ylim()[1] - plt.ylim()[0], color='k')
#plt.savefig("solstice.png")
plt.show()