I have a partial answer, sadly I can't get fast enough f-vector computations to confirm my suspicions.
Ziegler describes here some simple ways to get non-monotone f-vectors. Most of the constructions rely on projective transformation to get everything to fit. But in 0-1 polytopes, we can't do that.
Instead, we can distill the philosophy of this construction. If a polytope has:
- a region where the peak of that region's f-vector is low-dimensional
- a region where the peak of that region's f-vector is high-dimensional
- a small number of faces on the boundary of these regions
then the f-vector will have these two peaks.
This might be achievable in 0-1 polytopes. We can make a region that is simple, by keeping a part of the cube. Say, all vertices whose bitsum is less than or equal to $h$. We can make a region that is (close to) simplicial by choosing $n$ vertices at random. Then we just hope that the boundary is not too large.
These (random) polytopes have three parameters: dimension, height, and number of additional vertices. The higher the dimension, the further apart the peaks of the components.
For example: $d=13$, $h=3$, $n=3$
The f-vector of just the cube fragment (this is actually the matroid polytope for U(3,13)) is
- (1, 378, 5317, 20956, 48906, 80509, 97383, 87516, 58344, 28457, 9867, 2314, 338, 27, 1)
The f-vector when 3 vertices are added is
- (1, 381, 6177, 34618, 108515, 226559, 340412, 379265, 314077, 190136, 81180, 22918, 3791, 277, 1)
An f-vector with a similar number of faces where all vertices are chosen at random is
- (1, 23, 253, 1769, 8794, 32717, 92625, 198540, 316533, 365359, 293391, 153504, 46159, 5845, 1)
Notice that the peak in the f-vector of the matroid polytope is the 5-faces, whereas the peak of the f-vector for the random polytopes is the 8-faces. For the mixed polytope, the peak is the 5-faces.
I very strongly suspect with large $d$ and appropriate choices for $h$ and $n$, there are 0-1 polytopes with non-unimodular f-vectors.