Intuition for why sum of gaussian RVs is different from gaussian mixture I know that in the case of Gaussian mixture, the "intuition" is that you're drawing from a PDF which itself is just a sum of weighted Gaussian PDFs.
I don't understand the intuition behind how the sum of two iid Gaussian RVs leads to a Gaussian RV. I know the maths behind it, but I'm having a hard time visualizing how the sum of two Gaussian RVs leads to a Gaussian RV (and not a mixture).
 A: Forget the Gaussian part for a moment. Compare these two simple situations:
A) take a coin whose two sides are marked with 0 and 1 and a die with 20 sides numbered 1 to 20. Toss the coin and roll the die --- and add the results to get a total.
Consider these questions (and hints):
What's the chance you get a total of 0?  (you can't!)
What's the chance you get a total of 1?  (need 0 on the coin and a 1 on the die)
What's the chance you get a total of 11? (You need (0,11) or (1,10) here)
What's the chance you get a total of 21? (need (1,20) here)  
B) take the same coin and the same die and choose which one to use, by tossing a second coin labelled "coin" and "die". Now your total is the number showing on whichever object you tossed or rolled at the second step.
What's the chance you get a total of 0?  (get this by tossing ("coin",0) )
What's the chance you get a total of 1?  ( ("coin",1) or ("die",1) )
What's the chance you get a total of 11?  ( ("die",11) )
What's the chance you get a total of 21?  (you can't!)  
The first case (A) is a sum of random variables (convolution of the p.m.f.s). The second case (B) is a mixture (weighted average of p.m.f.s). 

They're entirely different kinds of things.
Now consider a Gaussian case. For example, if X ~ N(0,1) and Y ~ N(100,10) - where the two are independent - then X+Y has almost no chance of being between -1 and 1 (the contribution of Y makes the sum very far above values like those), but a 50-50 mixture of X and Y has quite a good chance of being between -1 and 1. 
In the case of iid Gaussians (as you suggest in comments), let's simplify further and take the example of the standard normal case (because it's the same up to scale and location shifts that don't alter the distribution shape of either the inputs or the output), though I'll talk about more general cases. First we have to see what calculating distributions for sums of random variables (particularly independent ones) involves; this is necessarily somewhat mathematical but I'll try to motivate what's going on rather than just do the algebra (which is easy to find done for the Gaussian case anyway so I'll just be explaining what happens rather than repeating what's easily found)
With discrete variables (as above) to compute the probability of a sum, you have to add the probabilities of all the distinct events that have that sum. So if $Z=X+Y$, $P(Z=z)=\sum_x P(X=x)P(Y=z-x)=\sum_x p_X(x)p_Y(z-x)$ (product because of independence of $X$ and $Y$, and $x$ ranges over the values possible for $x$ in that sum). 
This is discrete convolution.
Similarly for the continuous case -- we add over the cases that give the sum -- .
So $f_{Z}(z)=\int_x f_X(x)f_Y(z-x) dx$
... this is a convolution integral. In particular, see the animation in the second diagram here$^{[1]}$
For standard normals, that integral is just $\int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}}e^{-\frac12 x^2} \frac{1}{\sqrt{2\pi}}e^{-\frac12 (z-x)^2} dx$
But note that we have a quadratic-in-x in both exponents, so when we multiply them we're just adding the quadratics -- giving another quadratic (clearly this same idea applies for any independent Gaussians whatever their means and variances, and indeed, will generalize up to the jointly-normal case). 
So we end up with exp of (a quadratic in x with negative leading term). Because we can complete the square, we can split the quadratic into two terms -- the quadratic becomes  $k.(x-\text{something})^2$ + a term not in $x$ (which drops out the front of the integral, since as far as the integral in $x$ is concerned, it's a "constant"). 
In our simplified example that should look something like $2(x^2-z/2)^2+z^2/4$.
The first thing is (up to a scaling factor) a density in $x$, so if you pull out the required scaling factor from the scaling constant on the bivariate you started with you're left with an integral of a normal density. The remainder of the original scaling factor goes out the front, and the integral is $1$ (it's a density!). What's left out the front is $\exp$ of a quadratic in $z$ with a negative leading term (i.e. a Gaussian), times just the required scaling constant to make it a density. 
In short, the reason why convolutions of Gaussians are Gaussian is because products of $\exp$s are $\exp$s of sums and sums of quadratics are also quadratics. If your quadratics have negative leading terms (as they will for Gaussians), the sum is also going to have a negative leading term. So convolutions of Gaussians must be in the form of Gaussians.
Another simple/intuitive argument is as follows: Again consider the iid-bivariate Gaussian case $(X,Y)$, where the contours are circular. Without loss of generality, make it centered at the origin. Now if we rotate 45 degrees 
\begin{equation}
        \begin{pmatrix}
        U \\
        V \\
        \end{pmatrix} =
        \begin{pmatrix}
        \sqrt{\frac12} & \sqrt{\frac12} \\
        \sqrt{\frac12} &-\sqrt{\frac12} \\
        \end{pmatrix} 
        \begin{pmatrix}
        X \\
        Y \\
        \end{pmatrix} =
        \begin{pmatrix}
        \sqrt{\frac12}(X+Y) \\
        \sqrt{\frac12}(X-Y) \\
        \end{pmatrix} 
\end{equation}
we still have the same circular-bivariate Gaussian density centered at the origin. So $U$ and $V$ are also iid Gaussian with the same mean and variance as (X,Y). But $U$ is just a scaled $Z$, so $Z$ is also Gaussian. (We can generalize to the non-zero-mean case by simple shifts)
There are a number of proofs for the Gaussian case here; generally using characteristic functions (/Fourier transforms) or by using the convolution integral. I don't know that they're necessarily going to be intuitive for you. However, if you're used to characteristic functions/Fourier transforms that pretty much gives it to you immediately and might convey some of the intuition; in particular the Fourier transform of a Gaussian is Gaussian in form, and the product of two scaled Gaussian functions in the same variable is another scaled Gaussian, and convolution in the time domain equals multiplication in the frequency domain, so we see that the convolution of two Gaussians is Gaussian. 
In any case it's useful to have the proofs so you see in more detail where my prior handwaving about quadratics is coming from.
$[1]$  Weisstein, Eric W. "Convolution." From MathWorld--A Wolfram Web Resource. https://mathworld.wolfram.com/Convolution.html 
