Update: I expanded this answer into a paper Sorting with an average of $\mathrm{lg}(n!)+o(n)$ comparisons.
Yes, such an algorithm exists. I will only prove the $\mathrm{lg}(n!)+o(n)$ bound, but under a likely randomization assumption we also get $\mathrm{lg}(n!)+O(n^{1-ε})$. I will also describe an attempt for $n^{0.5+o(1)}$ and $O(n^{0.5-ε})$.
We can assume that all elements are distinct, by annotating them if necessary; the average case uses distinct elements in random order. We can compute the average number of comparisons by adding the entropy loss for each comparison relative to using a fair coin.
The starting point is insertion sort with a binary search to decide where to insert the next element into the sorted subset $S$. When $(1-ε)2^m ≤ |S| ≤ 2^m-1$, an insertion uses at most $m$ comparisons, which (in terms of entropy) is optimal up to an $O(ε)$ additive factor (and for average-case complexity, $2^m ≤ |S| ≤ (1+ε) 2^m$ also works). Now, when $|S|$ is not close to a power of 2, insertion of an element $A$ is suboptimal (even in the average case and regardless of how we balance each query), but if wasting $o(1)$ comparisons, we could steer $A$ to an approximately uniform distribution over an interval of $S$ of length close to a power of 2, we get the desired optimality.
We achieve this by adding elements in batches, and sometimes efficiently comparing elements of the batch with each other, such that the interval of $S$ corresponding to an element $A$ decreases in a quasi-random manner (and with the probability distribution of $A$ inside the interval nearly uniform), and when the interval length is close enough to a power of 2, doing the binary search to insert $A$.
Common constructs
We will keep a subset $S$ of sorted elements, and for each unsorted element $A$, we will keep track of the minimal interval $I_A$ of $S$ where $A$ is known to be located. $|I_A|$ is the length of $I_A$; $I_A=I_B$ is by the identity of the intervals.
Let $\mathrm{Compare}(A,B)$ be: Compare $A$ with $B$, and then (in random order) compare $A$ and $B$ against corresponding elements of $S$ until their intervals are disjoint (or have length 1). The element of $S$ is chosen (in a consistent manner) to make the probabilities for the comparison as close to 1/2 as possible, assuming that when $\mathrm{Compare}$ is called, $(A,B)$ is uniformly distributed on $I_A⨯I_B$. Because of the disjointness in the end, $\mathrm{Compare}$ preserves the uniformity assumption.
The following sections can be read independently of each other.
A $\mathrm{lg}(n!)+o(n)$ algorithm
Given: A sorted list $S$, and a batch of $m$ unsorted elements; $m∈ω(1)∩o(|S|)$; the unsorted elements are random relative to $S$.
Repeat (1)-(3) while possible:
1. Pick two elements $A$ and $B$ from the batch with $I_A=I_B$ (any choice will work).
2. Run $\mathrm{Compare}(A,B)$.
3. If $|I_A|$ is close enough to a power of 2,(note 1) remove $A$ from the batch (without forgetting $I_A$); and do similarly with $B$.
Finally: Insert all elements into $S$ and complete the sort.
Note 1: For "close enough", any $o(1)$ relative error (as a function of $m$) works as long as $m-o(m)$ elements will be removed in step (4) (possible by note 2). Under a conjectured randomization assumption, using $c \log \log m / \log m$ relative error captures $m(1-\log^{-Θ(c)}m)$ elements, allowing a $\mathrm{lg}(n!)+O(n \log \log n / \log n)$ average comparison sorting algorithm.
Note 2: Because the same sequence of comparisons leads to the same bounding interval, almost all elements will go through step (1) $Ω(\log m)$ times (unless removed in step 4). In the beginning, if $A < B$ and we pick $A$, we compare $A$ against element $S[≈(1-1/\sqrt{2})|S|]$, and each application of step (3) to $A$ has $O(1)$ probability of reducing $|I_A|$ in $≈1/(1-1/\sqrt{2})$ times. Now for every ratio $a>1$ that is not a rational power of 2, we have $∀ε>0 ∀d>0 ∃m,n∈\mathbb{N} \,\, 1-ε < \frac{a^m}{d2^n} < 1+ε$, and so we get the $o(n)$ bound.
A likely $\mathrm{lg}(n!)+O(n^{1-ε})$ algorithm
Modulo a randomization assumption, we can achieve $\mathrm{lg}(n!)+O(n^{1-ε})$ average comparisons as follows.
Randomly shuffle the items, and sort the first half into a list $S$, while keeping the second half as an unsorted batch.
Repeat until the batch is empty:
Randomly choose $A∈\text{batch}$. Let $G = \{ B∈\text{batch}: |P(A < B) - 0.5| < n^{-0.51ε} \}$. If $G$ is empty, remove $A$ from the batch and insert into $S$. Otherwise:
- If there is $B∈G$ such that with probability $Θ(1)$ (say ≥0.05), $\mathrm{Compare}(A,B)$ makes $|I_A|$ within $n^{-ε}$ relative error of a power of 2, run $\mathrm{Compare}(A,B)$ and if successful (i.e. $|I_A|$ is within $n^{-ε}$ relative error of a power of 2), remove $A$ from the batch and insert into $S$.
- If there is no such $B∈G$, run $\mathrm{Compare}(A,B)$ for a random $B∈G$.
If our randomization assumption works out (i.e. the distribution of interval lengths and positions is random enough), then throughout much of the process, a typical $A$ can be efficiently compared with a choice of $n^{Θ(1)}$ elements (with $n^{Θ(1)}$ different interval lengths). Thus, we can typically choose a comparison for (1) above, and if we are unlucky with the comparison result, we still get $Θ(\log n)$ chances, thus achieving (if $ε$ is small enough, say 0.01) a $\mathrm{lg}(n!)+O(n^{1-ε})$-comparison algorithm. With some changes and approximations, the total compute can be made quasilinear: Given an element $A$, compute promising interval lengths, and then look up $B$s with the right approximate center and interval lengths.
There are a number of ways to optimize comparisons, but the obstacle is that each comparison may end up being unlucky and we have a limited number of comparisons. If after optimization, $\mathrm{Compare}(A,B)$ does an average of 4 comparisons and 'succeeds' with 1/4 probability, we get $ε≈(1-ε)/4/\log_{4/3} 2 ≈ 0.09$.
A perhaps much better approach is to wait until an interval is close to a power of 2, controlling not the individual interval lengths but the distributions of lengths.
An attempt at a $\mathrm{lg}(n!)+n^{0.5+o(1)}$ algorithm
Suppose that $|S|=n$ and we are given an unsorted batch of $n$ elements with the intervals $I_A$ also given, with $|I_A|$ typically $n^{1-o(1)}$ and with $\frac{|I_A|}{2^{\lfloor \mathrm{lg} |I_A| \rfloor}}$ distributed uniformly (up to a random error, and holding with sufficient precision even if conditioned on $A < S[i]$). Then, we can sort the items wasting an average of $n^{0.5+o(1)}$ comparisons as follows:
(*) Insert all elements in the order of their initial $\frac{|I_A|}{2^{\lfloor \mathrm{lg} |I_A| \rfloor}}$. This way all elements are inserted when their interval length is close to a power of 2.
The sorting algorithm will be: Randomly shuffle the list and sort the first half $S$. To insert the second half, make the distribution right, and do the (*) above.
To make the $\frac{|I_A|}{2^{\lfloor \mathrm{lg} |I_A| \rfloor}}$ distribution right, we can make a 'random' distribution, and then withhold the right fraction of the elements for each $|I_A|/2^{\lfloor \mathrm{lg} |I_A| \rfloor}$ while randomizing the rest (repeating if necessary). However, while this should correct $\frac{|I_A|}{2^{\lfloor \mathrm{lg} |I_A| \rfloor}}$ globally, we do not know whether it can be controlled locally with the required precision (hence the word "attempt" above).
To make a 'random' distribution, we can randomly use $\mathrm{Compare}(A,B)$ with $P(A < B)≈0.5$, except that with the initial $I_A$ all identical, we do not expect randomization at a sublogarithmic depth (i.e. with $I_A$ long enough). However, I conjecture that we get randomization at a sublogarithmic depth using generalizations (probably any reasonable choice will work) of $\mathrm{Compare}$ to $k=ω(1)$ elements: If we keep $k=ω(1)$ elements entangled (i.e. connected using comparison results), we should have about $k$ noncommuting choices for each comparison with $S$. This should allow $O(\log_k n + \log k)$ randomization depth, as desired (assuming that $k$ is not too large as we need depth $Θ(\log k)$ to distentangle the elements). I expect that the compute can be made quasilinear if using a small enough $k$.
Since a comparison with $1/2+n^{-0.5}$ yes probability only wastes $O(1/n)$ entropy, the initial randomization and the slight nonuniformity of the elements in their bounding intervals should only need $n^{o(1)}$ entropy waste. If the distribution shaping succeeds well enough, the entropy waste stems primarily from interval length mismatches during (*) (hence the $n^{0.5+o(1)}$).
A possible $\mathrm{lg}(n!)+O(n^{0.5-ε})$ combination: If the distribution shaping works well enough and we make the batch size equal $|S|+n^{0.5+ε}$ and selectively reject $≈n^{0.5+ε}$ elements in (*) (above), we can insert all but these $≈n^{0.5+ε}$ elements with entropy waste $n^{0.5-ε/2+o(1)}$ as follows. Split $S$ into $n^ε$ nearly equal intervals, and when during insertion, $I_A$ settles on an interval, reject (i.e. cancel the insertion) if the interval is too long, thus reducing the variation in the lengths of these intervals $Θ(n^{ε/2})$ times, which in turn reduces the length variations of random length $n^{1-o(1)}$ intervals in $n^{ε/2-o(1)}$ times, as needed. Now, we can use the above $\mathrm{lg}(n!)+O(n^{1-ε})$ algorithm to insert the remaining elements with $O(n^{0.5-ε'})$ waste if $ε$ is small enough.
Worst-case complexity of sorting: Most likely, there is a sorting algorithm with $\mathrm{lg}(n!)+o(n)$ worst-case comparisons. For finding median, there is a linear gap between the average case ($1.5n+o(n)$ comparisons) and the worst case (at least $(2+ε)n-O(1)$ comparisons). However, for sorting, there is plenty of freedom for arranging comparisons and for finding new sorting algorithms.
