That's an excellent but highly unresolved question. The problem is that the *physicists* tend to be not so interested in mathematical foundations once a theory (statistical mechanics in this case) is successful and there is some plausible heuristic justification for it, whereas *ergodic theorists* quite often have little clue about the physical relevance and have an altogether different reason for being interested in Gibbs measures (starting with a "non-physical" analogy with statistical mechanics made by Sinai, Ruelle, ...). Meanwhile *mathematical physicists* have been busy with already highly challenging mathematical problems that arise if we take the physical meaning of Gibbs measures for granted.

I will tell you my take (à la Boltzmann) on how to approach the problem of the dynamical justification of Gibbs measures as "equilibrium states". Others will correct me if they don't agree.

Statistical mechanics is all about the relationship between the microscopic and macroscopic properties of systems. Let us say we have a set $X$ containing all the possible instantaneous *microscopic states* of a physical system, and a transformation $T:X\to X$ describing the time dynamics of the system, which we take to be discrete for simplicity. We also have a collection $\mathscr{O}$ of macroscopic *observables*, which are simply functions $f:X\to\mathbb{R}$. Knowing the *macroscopic state* of the system amounts to knowing the value of all the macroscopic observables, that is, a map $\pi:\mathscr{O}\to\mathbb{R}$. Each macroscopic state $\pi$ corresponds to a set $X_\pi$ of microscopic states that realize $\pi$, that is, $X_\pi=\lbrace x\in X: \text{$f(x)=\pi(f)$ for each $f\in\mathscr{O}$}\rbrace$.

A typical setting in statistical mechanics is a lattice model, in which the microscopic states of a physical system are represented by symbolic configurations on the lattice $\mathbb{Z}^d$, that is, by functions $x:\mathbb{Z}^d\to\Sigma$ for some finite set of symbols $\Sigma$. A natural candidate for the time dynamics $T$ is a cellular automaton, that is, a continuous map on $X=\Sigma^{\mathbb{Z}^d}$ that has all the translations as its symmetries (i.e., it commutes with the shifts). The choice of the macroscopic observables in this setting is somewhat debatable, but let's say the frequency of a finite word in the configuration is a macroscopic observable, and everything that can be written as a linear combination of such frequencies is also a macroscopic observable. Then, the macroscopic states can be represented by shift-invariant probability measures on $X$ (or better, shift-*ergodic* probability measures).

Suppose now that the system has a collection of conserved quantities like *energy*, that is, macroscopic observables $e\in\mathscr{O}$ such that $e(Tx)=e(x)$ "for all" $x\in X$. What can we say about the macroscopic state of the system? The second law of thermodynamics *suggests* that if the system is "in equilibrium", then the macroscopic state of the system has to maximize the "entropy" subject to the constraints imposed by the conservation laws, and if the system is not "in equilibrium", it gradually approaches the "equilibrium", provided the system is sufficiently "chaotic".

The physicists' explanation of the second law of thermodynamics in a finite-state system is that if the dynamics goes through all the configuration space in one cycle (it is "ergodic"), then in equilibrium (long long time after it has started its evolution) it is "most likely" to be found in a configuration whose macroscopic state encompasses the largest portion of the state space (i.e., $|X_\pi|$ is the largest), and if it is not in equilibrium, it is "most likely" to evolve in time towards macroscopic states with larger volume. In this case, the entropy of a macroscopic state $\pi$ is $\log|X_\pi|$.

Translating this heuristic to infinite systems (say in our lattice setting with cellular automaton dynamics) is not obvious. One is inclined to replace $|X_\pi|$ with the configuration space volume of $X_\pi$ according to some notion of volume, say the uniform Bernoulli measure in $\Sigma^{\mathbb{Z}^d}$ in case of a lattice model. Although this is how the mathematical notion of *ergodicity* came to being, it doesn't help with the problem of identifying macroscopic equilibrium states, because for all but a single macroscopic state $\pi$, the measure of $X_\pi$ with respect to the uniform Bernoulli measure is zero.

The correct way is to measure the the size of $X_\pi$ is by the amount of information *per site* required to describe a microscopic configuration in $X_\pi$, that is, by the entropy per site of $\pi$. Note that in the language of dynamical systems, this is the Kolmogorov-Sinai entropy of $\pi$ under the dynamics of the shift action (i.e., the *space dynamics*), and a priori has nothing to do with the time dynamics. Equilibrium statistical mechanics (e.g., the book of Robert Israel) tells us that the macroscopic states maximizing entropy under the constraints of the conservation laws are Gibbs measures.

A reasonable ergodic-theoretic justification of considering Gibbs measures as equilibrium states involves (1) showing that they are invariant under the time dynamics (this is the easy part) and (2) showing that starting from other macroscopic states the system evolves towards states with larger entropy at least under reasonable assumptions on the starting state, and provided the dynamics has a sufficiently strong chaotic behaviour.

Like I said at the beginning, this is largely open, but there is a beautiful example (and its extensions) for which some mathematical results have been found. The XOR cellular automaton is defined as $T:\lbrace 0,1\rbrace^{\mathbb{Z}}\to\lbrace 0,1\rbrace^{\mathbb{Z}}$ with $(TX)_i:=x_i+x_{i+1}\pmod{2}$. This system has no conserved quantity, so the second law suggests maximum entropy measure (the uniform Bernnoulli measure) as the only macroscopic equilibrium state. It is easy to see that the uniform Bernoulli measure is indeed invariant under $T$. Miyamoto and Lind (and later others) have shown that starting from any macroscopic state that has sufficient mixing property under the space dynamics (e.g., any Bernoulli or Markov measure), the system gradually approaches under $T$ towards the uniform Bernoulli measure (in a suitable sense).