In this answer, Ed Gorcenski says:
the Gaussian distribution (according to Wikipedia) was first studied by de Moivre.
and
Gauss used the distribution in astronomy
And Gorcenski uses this as an example of what he says often happens:
It seems that in many cases, naming the body of work was given to the person who first applied its study to some other field
The history is somewhat interesting, so I'm posting this present answer.
Much of this is based on things that I (pardon the expression) “read somewhere”, and a bit of it on an assertion in a course I once took.
Abraham de Moivre (1667–1754) wrote about the probability distribution of the number of “heads” that appear when a coin is tossed $1800$ times. As secondary school pupils are taught today, this is a binomial distribution that can be approximated by a normal, or “Gaussian” distribution with expectation $900$ and variance $450$ (thus with standard deviation $\sqrt{450} \approx 21$. Obviously computing precisely things like $ \sum_{k\,=\,880}^{905} \binom{1800}{k} \left(\frac12\right)^{1800}$ was prohibitively expensive, so de Moivre's discovery that areas under the bell-shaped curve $y=\text{constant}\times\exp(-z^2/2)$ approximates this was a major advance. I think that may be the first time anyone realized this function was so important. Computing areas under that curve is itself somewhat expensive in time and effort, but nowhere near computing sums like the one above.
According to Wikipedia, de Moivre discovered that $$n! \sim \text{constant} \cdot n^{n+1/2} e^{-n},$$ but did not discover that the “$\text{constant}$” was $\sqrt{2\pi\,}$, nor did he find that
$$
\int_{-\infty}^{+\infty} e^{-z^2/2} \, dz = \sqrt{2\pi\,}.
$$
It is somewhat routine to show that the “$\text{constant}$” referred to above and the value of this integral are the same; it is substantially simpler than actually finding that it's $\sqrt{2\pi\,}$. If I'm right about this so far, then de Moivre did find accurate numerical values for this without realizing that it's $\sqrt{2\pi\,}$.
Then de Moivre fled France to escape the persecution of Protestants in France, and went to England, where he met James Stirling. Stirling was the one who found that it's $\sqrt{2\pi\,}$. De Moivre then wrote a book about probability theory in English, titled The Doctrine of Chances, and included his results including, this time, $\sqrt{2\pi\,}$ rather than just something computed numerically.
“Optional” paragraph: That that phrase, the doctrine of chances, may have become 18th-century English for what today we call the theory of probability is suggested by the title of the famous posthumous paper of Thomas Bayes, “An Essay toward Solving a Problem in the Doctrine of Chances.” (The result in that paper, paraphrased in today's language, is that if $P$ is uniformly distributed between $0$ and $1$, and if $P$ then becomes the probability of success in a sequence of independent Bernoulli trials (so they are conditionally independent given $P$) then the conditional distribution of $P$ given the numbers of successes and failures is a certain beta distribution.)
Why do we use the standard deviation, which is the root-mean-square deviation from the mean, rather than the mean absolute deviation from the mean, which seems like the most obvious measure of dispersion? I think de Moivre was probably the first to understand the answer to that question. And he needed the answer for the coin-tossing problem above. It is this: The variance of the sum of independent random variables is the sum of their variances. Without that, de Moivre would not have know that the aforementioned standard deviation is $\sqrt{450}$.
So: What did Carl Gauss do?
The normal, or “Gaussian”, distribution with expectation $\mu$ and standard deviation $\sigma$ is
$$
\text{constant} \times \exp\left(-\frac12 \left( \frac{x-\mu} \sigma \right)^2 \right) \frac{dx}\sigma. \tag0\label{473075_0}
$$
The joint distribution for $n$ independent observations is proportional to
\begin{align}
& \prod_{i\,=\,1}^n \exp\left(-\frac12 \left( \frac{x_i-\mu} \sigma \right)^2 \right) \frac{dx_i}\sigma \tag1\label{473075_1} \\[8pt]
\propto {} & \underbrace{ \frac1{\sigma^n} \exp\left( -\frac1{2\sigma^2} \sum_{i\,=\,1}^n \left( x_i-\overline x\right)^2 \right) }_\text{No “$\mu$” appears here.} \, \exp\left( -\frac n{2\sigma^2} \left( \overline x-\mu \right)^2 \right)
\end{align}
where “$\propto$” means proportional as a function, not of the $x$s, but of $(\mu,\sigma)$, and $\overline x = (x_1+\dotsb+x_n)/n$.
Since “$\mu$” appears in this expression only in $(\overline x-\mu)^2$, the value of $\mu$ that maximizes \eqref{473075_1} is $\overline x$.
So here's the part I heard in an applied statistics course I once took: Gauss derived the expression \eqref{473075_0} from the fact that he wanted $\overline x$ to be the value of $\mu$ that maximizes \eqref{473075_1}. I.e. he wanted something that looked like $f((x-\mu)/\sigma)$, but what function should $f$ be? That was his question and line \eqref{473075_0} above was his answer.
Again, given the fact that I don't remember what the various sources that I read were, and that the instructor in the aformentioned course cited none, it's possible I'm substantially wrong about some of this, so tell me if I am.
