Does the Gibbs sampler converge to a global maximum in the presence of multiple modes? For example in case of a Gaussian mixture distribution?
First of all, I think you are conflating Gibbs-sampling-based simulated annealing (an optimization procedure that uses Gibbs sampling to draw the updates of the algorithm) and pure Gibbs sampling (a means to sample from a distribution).
Your question applies to simulated annealing that uses Gibbs sampling, but it does not make sense for general Gibbs sampling.
As Geman and Geman proved in their paper that introduced the widespread use of Gibbs sampling, a simulated annealing procedure based on a Metropolis/Gibbs type proposal distribution will converge to a stationary distribution that samples uniformly from all of the local modes, as long as the cooling schedule is slow enough.
So if you have an objective function that is the mixture of 3 Gaussians, say, with 3 distinct modes, then a simulated annealing process with a slow enough cooling schedule will converge to the uniform distribution on those three modes.
The cooling has to be inverse-logarithmic in the relative score of the objective function, which turns out to be way way too slow for practice. In practice, you give up the guarantee that you will converge to the distribution that samples all the modes (and hence would let you see the global minimum by simple inspection), and in return you get faster cooling schedules and approximate solutions with non-zero probability of getting stuck in a local mode that is not globally optimal.