Why are extreme correlations so common in this Monte Carlo simulation? Consider the following simple Monte Carlo:
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import pearsonr
r = np.arange(3)+1
R = [pearsonr(np.random.normal(0,1,3), r)[0] for i in range(5009)]
plt.hist(R, bins=100)
plt.xlabel('Pearson\'s Correlation Coefficient')
plt.ylabel('Frequency')
plt.show()

Which gives a diagram similar to the following:

In words, 5009 IID standard normal variables are sampled in triplicate (n=3) and are paired to the values (1,2,3). The code above computes the Pearson's product-moment correlation coefficient between the triplicate standard normal variable measures and the sequence (1,2,3). I am curious about the distribution that results, which appears to my surprise to be non-uniform. I suspect there is a mathematical reason behind this reproducible effect.
Is there a mathematically formal reason why the extreme correlations of $\pm 1$ are more common than the less extreme correlation scores?

Vector length seems to be too small, $n$, trying $30$. Minor changes to your sample code:
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import pearsonr
r = np.arange(30)+1
R = [pearsonr(np.random.normal(0,1,30), r)[0] for i in range(5009)]
plt.hist(R,bins=100)
plt.xlabel('Pearson\'s Correlation Coefficient')
plt.ylabel('Frequency')
plt.show()

Also, if we think that a noise should be uncorrelated, so correlations should look like normal centred at zero, not uniform:. This is the consequence of law of large numbers, intuitively.
Edit: Good point from @user20160. The length argument can be thought in terms of Autocorrelation Functions in time-series. Because observation is valid for two random variates too. So increasing length will decrease ACF to 0, mean value in the simulation. For small n=2, it is all fully correlated then it decays with increasing length.
