Ehrhart period collapse for $123\ldots k$-avoiding Birkhoff polytope?

I strongly suspect that the polytope in question is type of a (skew) Gelfand-Tsetlin polytope. The corresponding Ehrhart functions are polynomials, even though the polytopes are non-integral in general.

The classical Birkhoff polytope is an example of this; you can get it from skew Kostka polynomials associated with the skew shape, which consists of $n$ disconnected boxes. Dilation of the polytope with factor $k$ corresponds to having $n$ disjoint rows of length $k$.

Moreover, it is conjectured that the Ehrhart polynomials here always have non-negative coefficients (in fact, a more general conjecture states this for the Berenstein-Zelevinsky polytopes).

The Ehrhart function $L(\mathcal B^n_r;t)$ is an honest polynomial. We will show this by following Per's suggestion and proving that it coincides with the Ehrhart function of a certain Gelfand-Tsetlin polytope. The following two steps can probably be combined into one, but I thought it was natural to think of it this way.

Step 1: For every $k\in \mathbb N$, the lattice points of $k\mathcal B_r^n$ are in bijection with the set of arrays $\tau\in \mathbb Z_{\geq 0} ^{n\times n}$ that are plane partitions, meaning that $i\geq i', j\geq j'$ implies $\tau_{i,j}\le \tau_{i',j'}$, that also satisfy $kr\geq \tau_{1,1}$ and for each $-n\le s\le n$ we have the traces $$tr_s(\tau)=\sum_{i}\tau_{i,i+s}=k(n-|s|).$$ Proof: This bijection is simply RSK applied to the array $\pi$. The entry $\tau_{1,1}$ records the length of longest increasing subsequence in the biword recording $\pi$, which is exactly the largest possible sum of $\pi_{i,j}$'s along an upper left-lower right chain. This is where the restriction $kr\geq \tau_{1,1}$ comes from.

Let's introduce the partitions $\lambda^n_r=\{r^n0^n\}$ and $\mu^n_r=\{(r-1)^n1^n\}$.

Step 2: The plane partitions from the previous lemma are in bijection with the lattice points in the Gelfand-Tsetlin polytope $GT(k\lambda^n_r, k\mu^n_r)=kGT(\lambda^n_r, \mu^n_r)$, therefore they are counted by the Kostka number $K(k\lambda^n_r, k\mu^n_r)$.

Proof: If we write a Gelfand-Tseltin pattern with top row given by $k\lambda^n_r=\{kr,kr,\dots, kr, 0, 0,\dots, 0\}$ we can erase the leftmost entries that are forced to be $=kr$ and the rightmost entries that are forced to be $=0$ to be left with a ($45^{\circ}$ rotated) $n\times n$ array. This array is a plane partition with largest part $\le kr$ and the appropriate traces.

For $k\mathcal B^{4}_2$ this bijection looks like simply selecting the array of $*$'s. $$\begin{matrix} 2k && 2k && 2k && 2k && 0 && 0 && 0 && 0\\ & 2k && 2k && 2k && * && 0 && 0 && 0 \\ && 2k && 2k && * && * && 0 && 0\\ &&& 2k && * && * && * && 0 \\ &&&& * && * && * && * \\ &&&&& * && * && * \\ &&&&&& * && * \\ &&&&&&& * \end{matrix}$$

Combining both steps gives $L(\mathcal B^n_r ;t)=L(GT(\lambda^n_r, \mu^n_r);t)$ which implies the desired polynomiality. It is also expected that positivity of coefficients holds, but I believe it is not known yet.