Intuition behind Correlated Terms in Mixed Effect Models In mixed effect models, one can specify a varying intercept term and a varying slope term that covary with each other. The second to last row in this vignette by Douglas Bates shows how to specify such an interaction.

I do not have an intuitive understanding for what allowing this interaction means. 
1) What does it intuitively mean to allow these particular coefficients to covary? 
2) If I introduced another varying slope $s_2$ term into the model, in addition to $s_1$, what would it intuitively mean for $s_1$ and $s_2$ to vary?
3) Finally, what resources can I use to better understand how to perform inference on such correlations (in a manner similar to say Inference in linear mixed model
)?

What do I know:
I have an imprecise mathematical understanding of what correlations between coefficients mean. In particular, I recognize that a varying slope term and varying intercept term are random variables, and so the usual definition of covariance $E[(X - \mu_{X})(Y - \mu_{Y})]$ applies to them. I also understand that the goal of fitting a mixed effect model can be articulated as trying to find a covariance matrix, among other parameters, which is the covariance of a vector of model coefficients. Off diagonal terms in the covariance matrix represent correlations between model coefficients; thus, allowing (or prohibiting) correlation is tantamount to letting the off-diagonal terms be non-zero (or $0$).
 A: Suppose, for concreteness, that we have a model 
y ~ 1 + x1 + x2 + (1|g)

where x1 and x2 are (for simplicity) continuous predictor variables (g is a categorical grouping variable). This model states that the expected value of y changes linearly with changes in x1 and x2 and that there are differences in the intercept between groups (i.e. $\hat y = (\beta_0+b_{0,i}) + \beta_1 x_1 + \beta_2 x_2$), with $b_{0,i} \sim \textrm{Normal}(0,\sigma^2_0)$); that's what the 1 in (1|g) means.
If we change the random effect to (1|g) + (x1|g) + (x2|g) (separate terms), or equivalently (1+x1+x2||g), that specifies that the intercept, slope with respect to x1, and slope with respect to x2 all vary across groups, but this variation is independent: we could write this model out as
$$
\begin{split}
\hat y_{ij} & = (\beta_0 + b_{0,i}) + (\beta_1 + b_{1,i}) x_1 + (\beta_2 + b_{2,i}) x_2 \\
b_{0,i} & \sim \textrm{Normal}(0,\sigma^2_0) \\
b_{1,i} & \sim \textrm{Normal}(0,\sigma^2_1) \\
b_{2,i} & \sim \textrm{Normal}(0,\sigma^2_2) \quad .
\end{split}
$$
So a particular group might have a higher-than-average intercept, a lower-than-average response to $x_1$, and an average response to $x_2$, but all of these term-specific effects are independent of each other.
If we instead write (1+x1+x2|g), we can no longer specify models for $b_{k,i}$ separately: instead we have to write 
$$
\boldsymbol b_i = \{b_{0,i}, b_{1,i}, b_{2,i} \} \sim \textrm{MVN}(\boldsymbol 0, \Sigma)
$$
(where MVN is "multivariate normal"). Now in addition to the separate variances for each varying term ($\sigma^2_0$, $\sigma^2_1$, $\sigma^2_2$) we also have to specify the covariances or correlations ($\rho_{01}$, $\rho_{02}$, $\rho_{12}$). For example, suppose that $\rho_{12}$, the correlation between the $x_1$ and the $x_2$ slopes, is negative. That means that groups that respond strongly (positively) to changes in $x_1$ are likely to respond weakly, or even in the opposite direction, to changes in $x_2$ -- similar logic applies to $\rho_{01}$ and $\rho_{02}$ (correlations between among-group variation in the intercept and the two slopes).
You can compare the likelihood of a model with the full ("unstructured" or "general positive-definite") variance-covariance matrix to one with the diagonal (independent) variance-covariance matrix using likelihood ratio tests (or AIC): the independent-terms model is properly nested within the full model (i.e. starting with the full model and constraining $\rho_{01}=\rho_{02}=\rho_{12}$ gets you the reduced/nested model). Alternatively, you can find confidence intervals for individual $\rho_{ij}$ parameters. (nlme::lme does this by computing Wald intervals on a constrained (hyperbolic-tangent) scale and back-transforming; lme4::lmer does it by computing likelihood profiles).
