Do we also need to model a probability distribution for the decoder of a VAE?
Thanks @nbro for pointing this out.

The pictorial architecture in the slides uses the Gaussian loss, which when coupled with Maximum Likelihood Estimation gives the squared error loss (not to remove any tractability issues). The main reason we do the encoder Gaussian trick is to force the latent variable $z$ to be normal so that we can apply $KL$ $Divergence$ to optimise an otherwise intractable integral. You can get a better intuition and reasoning in this video. 

The pictorial architecture is basically taking the Gaussian loss so that the final loss becomes the squared error loss effectively. Also the loss term used in your blog link is exactly the same loss term used in the original paper, but the blog is using CE loss (it is the more common loss used for classification). I am not sure how they are using the CE loss though, since it is only valid for $0$ and $1$ values and AFAIK MNIST data-set has grayscale images. 

I am not exactly sure how they implement the randomness of the Gaussian loss in the decoder structure, but in the simplest of cases they just take the MSE

Check out this blog on VAE's (where they have taken the mean $\Sigma$ which they have abbreviated as mean, I have not checked their implementation detail to know what they exactly mean by that) and also this answer on Data Science on implementation of VAE's4 (both of which gives a more general form of loss). Also for the exact Mathematics check out Appendix C of the original paper.