As explained in this answer, the probability of a certain outcome in a Logit model can be written as $$ P=\int_{\varepsilon=-\beta'x}^{\infty} f(\varepsilon)d\varepsilon\\ = 1- F(-\beta'x) = 1-\dfrac{1}{1+\exp(\beta'x)} = \dfrac{\exp(\beta'x)}{1+\exp(\beta'x)} $$ I think I understand the reasoning behind this; the logit model assumes that the error term is distributed logistically, so the CDF of the logistic distribution is used. The probit model assumes that the error term follows a normal distribution giving
$$ \Pr(Y=1 \mid X) = \Phi(X'\beta) =\int_{\varepsilon=-\beta'x}^{\infty} f(\varepsilon)d\varepsilon = \Phi(x_{i}'\beta) = \int_{-\infty}^{x_{i}'\beta} \dfrac{1}{\sqrt{2\pi}}e^{-\dfrac{1}{2}t^{2}}dt $$ where $\Phi(x)$ is the CDF of the standard normal distribution. Now the referred answer, along with other sources, say that because these probabilities aren't closed-form, unlike those of the logit model, probit models require simulation.
But why use something complicated as simulation when there are a variety of methods that very accurately approximate $\Phi(x)$? It would seem to me that this would both be less complicated and more accurate.