Sampling from a $2$d normal with a given covariance matrix

How would one sample from the $2$-dimensional normal distribution with mean $0$ and covariance matrix $$\begin{bmatrix} a & b\\b & c \end{bmatrix}$$ given the ability to sample from the standard ($1$-dimensional) normal distribution?

This seems like it should be quite simple, but I can't actually find an answer anywhere.

Typically, given a covariance matrix $Q$, one perform the Cholesky decomposition on $Q$ i.e. get $Q = LL^T$ where $L$ is a lower triangular matrix.

Once we have $L$, generate a standard normal vector say $x$. Then $y=Lx$ is the desired random vector with covariance matrix $Q$.

This is because we have $$\mathbb{E}(y) = \mathbb{E}(Lx) = L\mathbb{E}(x) = L \times 0 = 0$$ and $$\mathbb{E}(yy^T)=\mathbb{E}(Lx(Lx)^T)=\mathbb{E}(Lxx^TL^T) = L\mathbb{E}(xx^T)L^T = L \times I \times L^T = LL^T$$

Note that any other symmetric square root decomposition of $Q$ also gives you the desired result i.e. if we can decompose $Q$ into $GG^T$, then $Gx$ will also gives us the random vector with zero mean and covariance $Q$.

In your case, if we perform Cholesky decomposition, we get $$L = \begin{bmatrix}\sqrt{a} & 0\\ \dfrac{b}{\sqrt{a}} & \dfrac{\sqrt{ac-b^2}}{\sqrt{a}} \end{bmatrix}$$ Now generate $x = \begin{bmatrix} x_1 \\ x_2\end{bmatrix}$ where $x_,x_2 \sim \mathcal{N}(0,1)$, then $$y = Lx = \begin{bmatrix}\sqrt{a} & 0\\ \dfrac{b}{\sqrt{a}} & \dfrac{\sqrt{ac-b^2}}{\sqrt{a}} \end{bmatrix} \begin{bmatrix}x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} \sqrt{a} x_1 \\ \dfrac{bx_1 + \sqrt{ac-b^2}x_2}{\sqrt{a}}\end{bmatrix}$$is the desired vector you are looking for.

Wikipedia has an outline. You need a square root (of sorts) $A$ of your covariance matrix, which you can get using the Cholesky decomposition. Then sample two iid $N(0,1)$ random variables, multiply them (as a vector) by $A$, and you get a bivariate draw from your bivariate normal distribution.

If you work in R, you can use the mvtnorm package.