Polynomial time mixing Markov chain for multimodal distribution 
Is there a discrete space Markov chain, starting from a fixed state, whose stationary distribution is a multimodal distribution and that mixes in polynomial time?

For example, Ising model on say a complete graph has a multimodal stationary distribution at low temperature. Critical $\beta$ (i.e. inverse temperature) is known to be $\beta_c=1$. So $\beta>1$ is low temperature regime. Single-site Glauber dynamics for Ising model at low temperature on a complete graph is known to mix exponentially slow. Here is a proof. I believe faster mixing Swendsen-Wang dynamics which update multiple spins in one step also mix exponentially slow at low temperature.
So I am wondering if there is a fundamental barrier and no Markov chain can mix in polynomial time, starting from a fixed state, at low temperature for Ising model on say a complete graph? Or is it that no one has found such a Markov chain yet?
Coupling is often used to prove upper bounds on mixing time. If there exists a contractive coupling, then a chain mixes in polynomial time. Here is a paper that proves that there exists a contractive coupling for a Markov chain on a continuous state space and a multimodal stationary distribution. For instance, see example 2.4 in this paper which talks about a Gaussian mixture.
 A: The OP seems specifically interested in ergodic Markov chains with a unique stationary distribution $\pi$.  The following example is admittedly a bit contrived, but is helpful to illustrate some general principles discussed more below.

Claim: There exists a Markov chain with a "bimodal" stationary distribution, such that the mixing time of the chain is logarithmic in $n$.

"Bimodal" just means there are two distinct states where $\pi$ is larger than in other states.
Proof. Let $\Delta_i(t)$ denote the total variation distance between a Markov chain starting at the $i$th state and the stationary distribution $\pi$.  For any accuracy parameter $\epsilon>0$, define the mixing time: $$
\tau(\epsilon) := \max_i \min\{ t > 0 : \Delta_i(t) < \epsilon \} \;.
$$ A chain that displays rapid mixing is given by the following lazy Markov chain with $n \times n$ transition matrix: $$
P_{ij} = \begin{cases} 1/2 & \text{if $i=j$} \\
1/4 & \text{if $i \notin \{1,n \}$ and $j \in \{1, n \}$} \\ 1/(2 (n-1)) & \text{if $i \in \{1,n\}$ and $j\ne i$} \\
0 & \text{else} 
\end{cases}
$$ The stationary distribution of this chain satisfies $$
\pi \propto \left( 1, \frac{2}{n-1}, \cdots, \frac{2}{n-1},1 \right) 
$$ which is bimodal (concentrated at the endpoints), and with spectral gap $\gamma=1/2$, and since the chain is reversible, the mixing time satisfies $$
\tau(\epsilon) \le 2 \log\left(\frac{2 n-3}{\epsilon}\right) 
$$ which is logarithmic in $n$. $\Box$
In general, energy and/or entropic barriers can cause chains that make only local moves to mix with exponential time.  This is, e.g., the case with Glauber dynamics applied to the mean-field Ising (or Curie-Weiss) model at low temperature; see, e.g. Theorem 3 of Levin/Luczak/Peres (2010) and references therein.  In contrast, the above chain is contrived to make global moves, and doesn't really have any type of barriers, which is why it mixes so fast.
A: As noted in the comment by James Martin, some assumption is needed on the Markov chain, e.g.,  that each step of the chain can be implemented on a Turing machine in polynomial time.
The Swendsen-Wang dynamics for the (Ferromagnetic) Ising model on a complete graph $K_n$ mixes in at most $O(n^{1/4})$ iterations, see [1]. More generally, on any graph, this dynamics has polynomial mixing time [2], although the sharp exponent is not known; it is possible that the bound  $O(n^{1/4})$ from [1] holds true for all graphs, but the known bounds are far greater.
The Swendsen-Wang dynamics for the Potts model on a complete graph $K_n$ has exponential mixing time at low temperatures [3].
perhaps most pertinent is the work of Alan Sly [4] that relates sampling to computational hardness.
[1] Long, Yun, Asaf Nachmias, Weiyang Ning, and Yuval Peres. A power law of order 1/4 for critical mean field Swendsen-Wang dynamics. American Mathematical Soc., 2014.   https://www.google.com/books/edition/A_Power_Law_of_Order_1_4_for_Critical_Me/YNlVBQAAQBAJ?hl=en&gbpv=1&pg=PP1&printsec=frontcover
[2] Guo, Heng, and Mark Jerrum. "Random cluster dynamics for the Ising model is rapidly mixing." In Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms, pp. 1818-1827. Society for Industrial and Applied Mathematics, 2017.
[3] Gheissari, Reza, Eyal Lubetzky, and Yuval Peres. "Exponentially slow mixing in the mean-field Swendsen-Wang dynamics." In Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms, pp. 1981-1988. Society for Industrial and Applied Mathematics, 2018.
[4] Sly, Allan. "Computational transition at the uniqueness threshold." In 2010 IEEE 51st Annual Symposium on Foundations of Computer Science, pp. 287-296. IEEE, 2010.
