Why is the sum of all the elements in a Gaussian-distributed list with zero mean not zero? If I generate a list of elements from a Gaussian distribution with zero mean using Python
List = np.random.normal(0, 1, 500)

my intuition (why is obviously wrong) tells me that when we sum all the elements up, the result would be close to zero
np.sum(List) ≃ 0

just like the integral of a sine function is zero when we integrate it across one period.
Where does my intuition go wrong?
 A: Suppose we are simulating $N$ values from a standard normal distribution, as you do. With some notation, we can write this as $X_1, \ldots , X_N \sim N(0, 1)$ (iid).
In this case, the expected value of the sum of these values is indeed zero:
$$E\left(\sum_{i=1}^{N} X_i\right) = \sum_{i=1}^{N} E(X_i) = \sum_{i=1}^{N} 0 = 0$$
However, when you add $N$ standard normal independent random variables, the sum of these terms has an additive variance:
$$\text{Var}\left(\sum_{i=1}^{N} X_i \right) = \sum_{i=1}^{N} \text{Var}(X_i) = \sum_{i=1}^{N} 1 = N$$
So when you add many terms, the expected value is zero, but the variance is extremely large (and gets larger as you add more and more terms). So in any given simulation, you will find that your sum of simulated data points is very different from zero.
One test you can do is to repeat your process of simulating 500 values, say 10,000 times, and plot a histogram of your simulated sums. You will probably see it centered around zero, with an estimated variance close to 500.
In Python, the test would look like this:
import numpy as np

sum_values = np.empty(10000)
for i in range(10000):
    sim_values = np.random.normal(0, 1, 500)
    sum_values[i] = np.sum(sim_values)

sum_values.mean()
sum_values.var()

A: Your intuition isn't too bad! In a certain sense, you should indeed expect np.sum(List)=0. More precisely, this sum's mean is zero since the expectation is linear. However, as you note, this is often not the case because this sum's variance is quite large.
Suppose you take $n$ samples of IID standard Gaussians (mean 0, variance 1 in your code). In that case, the standard deviation is equal to $\sqrt{n}$ (and in the more general case where the Gaussians have variance $\sigma^2$ the standard deviation goes as $\sqrt{n} \sigma$). Using the handy 68-95-99.7 Rule for Gaussian distributions (since the sum of Gaussians is still Gaussian), we can see that when $n=500$ about $\frac{1}{3}$ of the time the absolute value of np.sum(List) should be greater than $\sqrt{n} = \sqrt{500} = 22$.
To "correct" this behavior, i.e., to ensure that the sum stays "close" to its mean, we have to bring down its variance by dividing the sum by $n$ (or equivalently multiplying the sum by $\frac{1}{n}$). We are then guaranteed that the scaled sum will stay close to its mean value, in this case 0. This is related to your point that the integral of a sine function is zero when integrated across a period. Recall that when you integrate a function, you're not just summing its values over a specific region; first, you must "multiply" it by an infinitesimal area $dx$ (speaking roughly). Here you can think about the $\frac{1}{n}$ factor as the "discrete" version of the infinitesimal. This is made rigorous by the field of measure theory, where the factor of $\frac{1}{n}$ is called the "empirical measure"
