  We investigate the effect of quenched bond disorder on the two-dimensional
three-color Ashkin-Teller model, which undergoes a first-order phase transition
in the absence of impurities. This is one of the simplest and striking models
in which quantitative numerical simulations can be carried out to investigate
emergent criticality due to disorder rounding of first-order transition.
Utilizing extensive cluster Monte Carlo simulations on large lattice sizes of
up to $128 \times 128$ spins, each of which is represented by three colors
taking values $\pm 1$, we show that the rounding of the first-order phase
transition is an emergent criticality. We further calculate the correlation
length critical exponent, $\nu$, and the magnetization critical exponent,
$\beta$, from finite size scaling analysis. We find that the critical
exponents, $\nu$ and $\beta$, change as the strength of disorder or the
four-spin coupling varies, and we show that the critical exponents appear not
to be in the Ising universality class. We know of no analytical approaches that
can explain our non-perturbative results. However our results should inspire
further work on this important problem, either numerical or analytical.
