The Gaussian copula (with correlation matrix $R$):
$$
C_R(u)=\Phi_n(\Phi^{-1}(u_1),\dots,\Phi^{-1}(u_n);R).
$$
Then, letting $\zeta:=(\Phi^{-1}(u_1),\dots,\Phi^{-1}(u_n))^{\top}$,
$$
c_R(u)=\frac{\phi_n(\Phi^{-1}(u_1),\dots,\Phi^{-1}(u_n);R)}{\phi(\Phi^{-1}(u_1))\cdots \phi(\Phi^{-1}(u_1))} \\=(2\pi)^{-n/2}\lvert R \rvert^{-1/2}\frac{-\exp(\zeta^{\top}R^{-1}\zeta/2)}{\prod_{i=1}^n (2\pi)^{-1/2}\exp(-\zeta_i^2/2)}=\lvert R \rvert^{-1/2}\frac{\exp(-\zeta^{\top}R^{-1}\zeta/2)}{\exp(-\zeta^{\top}I_n\zeta/2)}.
$$
Here, $\Phi_n(\cdot;R),\phi_n(\cdot;R)$ are zero-mean multivariate normal CDF and PDF with covariance matrix $R:n\times n$.
In general, there is (as far as I know) only one universally applicable sampling method.
It goes by the name of "conditional distribution method", "conditional inverse method" or "conditional sampling".
Depending on the specific type of copula, other simulation methods can be available.
1) Conditional distribution method
I'll limit myself to the 2D case, and specifically I follow chapter 2.9 (page 40ff) in: Nelsen. An introduction to copulas. Springer 2006, 2nd edition.
For random variables $(U,V) \sim C$ (i.e. $U$ and $V$ have a the copula you want to simulate from as joint distribution) denote by $c_u$ the conditional distribution of $V$ given $U=u$.
It can be shown that $c_u$ has the following form
\begin{align}
c_u(v) = \Pr[V\leq v | U = u] = \frac{\partial C(u,v)}{\partial u}.
\end{align}
The function $c_u$ is furthermore non-decreasing, implying that a (generalized) inverse $c_u^{\leftarrow}$ exists.
You can now simulate from $C$ as follows:
- Generate two independent standard uniform random numbers $u$ and
$t$.
- Set $v = c_u^{\leftarrow}(t)$.
- The pair $(u,v)$ is now a sample from $C$.
In principle this method also works in the $d$-dimensional case, but it can in general be quite difficult to derive the necessary inverses of the conditional distribution functions.
The $d$-dimensional description is for example in
Cherubini, Luciano, Vecchiato. Copula Methods in Finance. John Wiley & Sons, Ltd, 2004, chapter 6.3 (page 183ff),
but is (in my opinion) often not the most practical way to simulate.
So instead of the general approach, specialized approaches depending on the specific copula in question are often the easier way to go.
2) Copula specific methods
Depending on the definition of the copula, more efficient simulation methods can be available.
For example, in case of the Gaussian copula $C_R$ you can simply simulate from a multivariate normal distribution with standard normal margins $\Phi$ and variance-covariance matrix $R$. For a realization $\mathbf{x} = (x_1,\ldots,x_d)$ you can now apply $\Phi$ to each component to get $\mathbf{u} = (\Phi(x_1),\ldots,\Phi(x_d))$. Then $\mathbf{u}$ is a sample from the Gaussian copula $C_R$ with parameter matrix $R$.
The same goes for the $t$-copula, where simulating from a multivariate $t$-distribution and then transforming with the margins is easier than simulating from the copula directly.
For (nested) Archimedean copulas simulation is tied to the Laplace transform of the copula generator.
You can find the details in:
Hofert, M. (2007). Sampling Archimedean copulas (https://www.uni-ulm.de/fileadmin/website_uni_ulm/mawi/forschung/PreprintServer/2007/preprintmariushofert.pdf)
Best Answer
Let $x=[\Phi^{-1}(u_1),\ldots,\Phi^{-1}(u_n)]^{\top}$. Then \begin{align} c(u_1,\ldots,u_n;\Sigma)&=\frac{\partial^nC(u_1,.\ldots,u_n;\Sigma)}{\partial u_1\cdots\partial u_n}=\frac{\Phi(x;0,\Sigma)}{\prod_{i=1}^n \phi(x_i)} \\ &=(2\pi)^{-\frac{n}{2}}|\Sigma|^{-\frac{1}{2}}\exp\!\left(-\frac{1}{2}x^{\top}\Sigma^{-1}x\right)\times\prod_{i=1}^n (2\pi)^{-\frac{1}{2}}\exp\!\left(\frac{1}{2}x_i^2\right) \\ &=|\Sigma|^{-\frac{1}{2}}\exp\!\left(-\frac{1}{2}x^{\top}\Sigma^{-1}x+\frac{1}{2}x^{\top}x\right) \\ &=|\Sigma|^{-\frac{1}{2}}\exp\!\left(-\frac{1}{2}x^{\top}(\Sigma^{-1}-I)x\right). \end{align}