Stanley's non-Ehrhart-positive order polytopes

We say a polytope is Ehrhart positive if all the coefficients in its Ehrhart polynomial are positive. Answering an Ehrhart positivity question posed on Mathoverflow, Stanley provided an example of a non-Ehrhart-positive order polytope of dimension $21$. Stanley's example comes from a certain family of order polytopes. In this paper, we study the Ehrhart positivity question on this family of polytopes. By giving explicit formulas for the coefficients of the Ehrhart polynomials of these polytopes in terms of Bernolli numbers, we determine the sign of each Ehrhart coefficient of each polytope in the family. As a consequence of our result, we conclude that for any positive integer $d \ge 21,$ there exists an order polytope of dimension $d$ that is not Ehrhart positive, and for any positive integer $\ell$, there exists an order polytope whose Ehrhart polynomial has precisely $\ell$ negative coefficients, which answers a question posed by Hibi. We finish this article by discussing the existence of lower-dimensional order polytopes whose Ehrhart polynomials have a negative coefficient.