I'll describe below a controversy in statistical mechanics in the 1980's: the case of the lower critical dimension of the Ising model with an applied random magnetic field.
Background
Let me give a little background, though you might want to read Terry Tao's discussion of basic statistical mechanics instead. The Ising model is a statistical mechanical model of "spins" on a hypercubic lattice. The energy functional is: $E=\sum_{\langle ij\rangle}\frac{1}{2}(1-S_iS_j)-\sum h_iS_i$ where the first sum is taken over nearest neighbor pairs on the lattice and the second is taken over all sites, and $S_i$ is a $\pm1$ valued variable on each site called the spin and $h_i$ is the real-valued "externally applied magnetic field" applied to each site. Each possible configuration of spins on the lattice is assigned a probability proportional to its Boltzmann weight $e^{-\beta E}$ where $\beta>0$ is a parameter that is interpreted physically as the inverse temperature $T$.
Given such a model, one question is to determine the "phase behavior", or what are the typical properties of the ensemble of configuration at a given $\beta$, and how does this change with $\beta$.
Considering at the moment just the Ising model with $h_i=0$, one might expect that for large $\beta$, the typical configuration will tend to have lower energy, and hence have all its spins aligned to either all $+1$ or all $-1$. At small $\beta$, all the Boltzmann factors tend to 1 and the typical configuration will have random spins. This rough argument is just meant to guide the intuition that there might be a phase transition between "mostly aligned" configurations to "mostly random" configurations at some special value of $\beta$.
As it turns out, what happens is highly dependent on the dimensionality of the lattice.
The lower critical dimension $d_L$ of a model is the dimension below which no phase transitions can occur because even as $\beta\rightarrow\infty$, there is not enough of an energy gain from ordering to create a phase with long-range correlations. In the ordinary Ising model (with all $h_i=0$), the lower critical dimension is 1, and hence at any finite $\beta$, the average $\langle S_iS_{j}\rangle$ over configurations weighted with the Boltzmann distribution will approach zero (exponentially fast, even) as the distance between sites $i$ and $j$ approaches $\infty$. For two dimensions and above, it can be shown that above a certain $\beta_c$ (depending on dimension) this average will be finite in that long-distance limit.
Controversy
In the 1980's there was a controversy in the physics literature over the value of $d_L$ for the Random Field Ising model, a model where the $h_i$ are independent Gaussian random variables with zero mean and constant variance $\epsilon^2$.
I'm not in a position to describe the history accurately, but I believe that there were physical arguments by Imry and Ma originally that $d_L\leq 2$, which were disputed when an amazing connection between random systems in $d$ dimensions and their pure counterparts in $d-2$ dimensions was found, known as the "Parisi-Sourlas correspondence". My understanding of Parisi-Sourlas is that it is based on a hidden supersymmetry in some series representation of the model which yields order-by-order agreement in the "epsilon expansions" of the two systems. Their argument was also made rigorous by Klein, Landau and Perez (MR). Based on this, since the Ising model has $d_L=1$, the RFIM was argued to have $d_L=3$ by various authors, though this was never a consensus view.
This controversy was settled by work of John Imbrie (MR) and later work of Bricmont and Kupianen (MR) building off his results that proved rigorously that $d_L\leq2$ in this system. Apparently terms like $e^{-1/\epsilon}$ become important and the epsilon expansion breaks down in low dimensions, though I'm not sure if this has been made precise, and even today the RFIM is far from being completely understood.