  We use Monte Carlo techniques and analytical methods to study the phase
diagram of multicomponent Widom-Rowlinson models on a square lattice: there are
M species all with the same fugacity z and a nearest neighbor hard core
exclusion between unlike particles. Simulations show that for M between two and
six there is a direct transition from the gas phase at z < z_d (M) to a demixed
phase consisting mostly of one species at z > z_d (M) while for M \geq 7 there
is an intermediate ``crystal phase'' for z lying between z_c(M) and z_d(M). In
this phase, which is driven by entropy, particles, independent of species,
preferentially occupy one of the sublattices, i.e. spatial symmetry but not
particle symmetry is broken. The transition at z_d(M) appears to be first order
for M \geq 5 putting it in the Potts model universality class. For large M the
transition between the crystalline and demixed phase at z_d(M) can be proven to
be first order with z_d(M) \sim M-2 + 1/M + ..., while z_c(M) is argued to
behave as \mu_{cr}/M, with \mu_{cr} the value of the fugacity at which the one
component hard square lattice gas has a transition, and to be always of the
Ising type. Explicit calculations for the Bethe lattice with the coordination
number q=4 give results similar to those for the square lattice except that the
transition at z_d(M) becomes first order at M>2. This happens for all q,
consistent with the model being in the Potts universality class.
