As a former physicist I can see how it could have been derived. This is how physicists proceed:
when they encounter a finite integral of a positive function, such as beta function: $$B(x,y) = \int_0^1t^{x-1}(1-t)^{y-1}\,dt$$ they instinctively define a density: $$f(s|x,y)=\frac{s^{x-1}(1-s)^{y-1}}{\int_0^1t^{x-1}(1-t)^{y-1}\,dt}=\frac{s^{x-1}(1-s)^{y-1}}{B(x,y)},$$ where $0<s<1$
They do this to all kinds of integrals all the time so often that it happens reflexively without even thinking. They call this procedure "normalization" or similar names. Notice how by definition trivially the density has all the properties that you want it to have, such as always positive and adds up to one.
The density $f(t)$ that I gave above is of Beta distribution.
UPDATE
@whuber's asking what's so special about Beta distribution while the above logic could be applied to an infinite number of suitable integrals (as I noted in my answer above)?
The special part comes from the binomial distribution. I'll write its PDF using similar notation to my beta, not the usual notation for parameters and variables: $$ f'(x,y|s) = \binom {y+x} x s^x(1-s)^{y}$$
Here, $x,y$ - number of successes and failures, and $s$ - probability of success. You can see how this is very similar to the numerator in the Beta distribution. In fact, if you look for the prior for Binomial distribution, it'll be the Beta distribution. It's not surprising also because the domain of Beta is 0 to 1, and that's what you do in Bayes theorem: integrate over the parameter $s$, which is the probability of success in this case as shown below: $$\hat f(x|X)=\frac{f'(X|s)f(s)}{\int_0^1 f'(X|s)f(s)ds},$$ here $f(s)$ - probability (density) of probability of success given the prior settings of Beta distribution, and $f'(X|s)$ - density of this data set (i.e. observed success and failures) given a probability $s$.