A standard linear model (e.g., a simple regression model) can be thought of as having two 'parts'. These are called the *structural component* and the *random component*. For example:  
$$
Y=\beta_0+\beta_1X+\varepsilon  \\
\text{where }  \varepsilon\sim\mathcal{N}(0,\sigma^2)
$$
The first two terms (that is, $\beta_0+\beta_1X$) constitute the structural component, and the $\varepsilon$ (which indicates a normally distributed error term) is the random component. When the response variable is not normally distributed (for example, if your response variable is binary) this approach may no longer be valid. The [generalized linear model][1] (GLiM) was developed to address such cases, and logit and probit models are special cases of GLiMs that are appropriate for binary variables (or multi-category response variables with some adaptations to the process). A GLiM has three parts, a *structural component*, a *link function*, and a *response distribution*. For example:  
$$
g(\mu)=\beta_0+\beta_1X
$$
Here $\beta_0+\beta_1X$ is again the structural component, $g()$ is the link function, and $\mu$ is a mean of a conditional response distribution at a given point in the covariate space. The way we think about the structural component here doesn't really differ from how we think about it with standard linear models; in fact, that's one of the great advantages of GLiMs. Because for many distributions the variance is a function of the mean, having fit a conditional mean (and given that you stipulated a response distribution), you have automatically accounted for the analog of the random component in a linear model (N.B.: this can be more complicated in practice).  

The link function is the key to GLiMs: since the distribution of the response variable is non-normal, it's what lets us connect the structural component to the response--it 'links' them (hence the name). It's also the key to your question, since the logit and probit are links (as @vinux explained), and understanding link functions will allow us to intelligently choose when to use which one. Although there can be many link functions that can be acceptable, often there is one that is special. Without wanting to get too far into the weeds (this can get very technical) the predicted mean, $\mu$, will not necessarily be mathematically the same as the response distribution's *canonical location parameter*; the link function that does equate them is the *canonical link function*. The advantage of this "is that a minimal sufficient statistic for $\beta$ exists" ([German Rodriguez][2]). The canonical link for binary response data (more specifically, the binomial distribution) is the logit.  However, there are lots of functions that can map the structural component onto the interval $(0,1)$, and thus be acceptable; the probit is also popular, but there are yet other options that are sometimes used (such as the complementary log log, $\ln(-\ln(1-\mu))$, often called 'cloglog'). Thus, there are lots of possible link functions and the choice of link function can be very important.  The choice should be made based on some combination of:  

 1. Knowledge of the response distribution,  
 2. Theoretical considerations, and  
 3. Empirical fit to the data.  

Having covered a little of conceptual background needed to understand these ideas more clearly (forgive me), I will explain how these considerations can be used to guide your choice of link. (Let me note that I think @David's [comment][3] accurately captures why different links are chosen *in practice*.)  To start with, if your response variable is the outcome of a Bernoulli trial (that is, $0$ or $1$), your response distribution will be binomial, and what you are actually modeling is the probability of an observation being a $1$ (that is, $\pi(Y=1)$). As a result, any function that maps the real number line, $(-\infty,+\infty)$, to the interval $(0,1)$ will work.  

From the point of view of your substantive theory, if you are thinking of your covariates as *directly* connected to the probability of success, then you would typically choose logistic regression because it is the canonical link. However, consider the following example: You are asked to model `high_Blood_Pressure` as a function of some covariates. Blood pressure itself is normally distributed in the population (I don't actually know that, but it seems reasonable prima fascie), nonetheless, clinicians dichotomized it during the study (that is, they only recorded 'high-BP' or 'normal'). In this case, probit would be preferable a-priori for theoretical reasons. This is what @Elvis meant by "your binary outcome depends on a hidden Gaussian variable". Another consideration is that both logit and probit are *symmetrical*, if you believe that the probability of success rises slowly from zero, but then tapers off more quickly as it approaches one, the cloglog is called for, etc.  

Lastly, note that the empirical fit of the model to the data is unlikely to be of assistance in selecting a link, unless the shapes of the link functions in question differ substantially (of which, the logit and probit do not). For instance, consider the following simulation:  

    set.seed(1)
    probLower = vector(length=1000)
    
    for(i in 1:1000){      
        x = rnorm(1000)
        y = rbinom(n=1000, size=1, prob=pnorm(x))
      
        logitModel  = glm(y~x, family=binomial(link="logit"))
        probitModel = glm(y~x, family=binomial(link="probit"))
       
        probLower[i] = deviance(probitModel)<deviance(logitModel)
    }
    
    sum(probLower)/1000
    [1] 0.695

Even when we know the data were generated by a probit model, and we have 1000 data points, the probit model only yields a better fit 70% of the time, and even then, often by only a trivial amount. Consider the last iteration:  

    deviance(probitModel)
    [1] 1025.759
    deviance(logitModel)
    [1] 1026.366
    deviance(logitModel)-deviance(probitModel)
    [1] 0.6076806

The reason for this is simply that the logit and probit link functions yield very similar outputs when given the same inputs.  

![Enter image description here][4]  

The logit and probit functions are practically identical, except that the logit is slightly further from the bounds when they 'turn the corner', as @vinux stated. (Note that to get the logit and the probit to align optimally, the logit's $\beta_1$ must be $\approx 1.7$ times the corresponding slope value for the probit. In addition, I could have shifted the cloglog over slightly so that they would lay on top of each other more, but I left it to the side to keep the figure more readable.) Notice that the cloglog is asymmetrical whereas the others are not; it starts pulling away from 0 earlier, but more slowly, and approaches close to 1 and then turns sharply.  

A couple more things can be said about link functions. First, considering the *identity function* ($g(\eta)=\eta$) as a link function allows us to understand the standard linear model as a special case of the generalized linear model (that is, the response distribution is normal, and the link is the identity function). It's also important to recognize that whatever transformation the link instantiates is properly applied to the *parameter* governing the response distribution (that is, $\mu$), not the actual response *data*. Finally, because in practice we never have the underlying parameter to transform, in discussions of these models, often what is considered to be the actual link is left implicit and the model is represented by the *inverse* of the link function applied to the structural component instead. That is:  
$$
\mu=g^{-1}(\beta_0+\beta_1X)
$$
For instance, logistic regression is usually represented:
$$
\pi(Y)=\frac{\exp(\beta_0+\beta_1X)}{1+\exp(\beta_0+\beta_1X)}
$$
instead of:
$$
\ln\left(\frac{\pi(Y)}{1-\pi(Y)}\right)=\beta_0+\beta_1X
$$  

For a quick and clear, but solid, overview of the generalized linear model, see chapter 10 of [Fitzmaurice, Laird, & Ware (2004)][5], (on which I leaned for parts of this answer, although since this is my own adaptation of that--and other--material, any mistakes would be my own). For how to fit these models in R, check out the documentation for the function [?glm][6] in the base package.  

*(One final note added later:)* I occasionally hear people say that you shouldn't use the probit, because it can't be interpreted. This is not true, although the interpretation of the betas is less intuitive. With logistic regression, a one unit change in $X_1$ is associated with a $\beta_1$ change in the log odds of 'success' (alternatively, an $\exp(\beta_1)$-fold change in the odds), all else being equal. With a probit, this would be a change of $\beta_1\text{ }z$'s. (Think of two observations in a dataset with $z$-scores of 1 and 2, for example.) To convert these into predicted *probabilities*, you can pass them through the normal [CDF][7], or look them up on a $z$-table.  

(+1 to both @vinux and @Elvis. Here I have tried to provide a broader framework within which to think about these things and then using that to address the choice between logit and probit.)

  [1]: http://en.wikipedia.org/wiki/Generalized_linear_model
  [2]: http://data.princeton.edu/wws509/notes/a2s1.html
  [3]: http://stats.stackexchange.com/questions/20523/difference-between-logit-and-probit-models#comment37049_20529
  [4]: https://i.sstatic.net/glarC.png
  [5]: http://www.amazon.com/Applied-Longitudinal-Analysis-Probability-Statistics/dp/0471214876
  [6]: http://web.njit.edu/all_topics/Prog_Lang_Docs/html/library/base/html/glm.html
  [7]: https://en.wikipedia.org/wiki/Cumulative_distribution_function