Covariance between a normal variable (x1) and a sum including it (x1 + x2) I'm interested in the joint distribution of two variables, 
$x_1$ and $x_2$.
$$
x_1 \sim Normal(0, \sigma^2_1);\\
\epsilon \sim Normal(0, \sigma^2_{\epsilon});\\
x_2 = x_1 + \epsilon;
$$
as a bivariate normal distribution
$$
\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}
 \sim Normal( 0, \Sigma); \\
$$
$$
\Sigma =  \begin{bmatrix} 
  \sigma^2_1 &  \rho\sigma_1\sigma_2 \\ 
  \rho\sigma_1\sigma_2 & \sigma^2_2
\end{bmatrix}
$$
$$
\sigma^2_2 = \sigma^2_1 + \sigma^2_{\epsilon}
$$
I can simulate from this distribution quite easily, 
but would like the closed-form solution for the correlation coefficient $\rho$.
I can see that the $\rho$ should be proportional to
 $\sigma^2_1$: 
with higher $\sigma^2_1$, the signal present in $x_1$ is stronger and less influenced by the noise $\epsilon$.
Simulations show that the covariance is actually very close to the variance $\sigma^2_1$...
library(tidyverse)
vals = 2^seq(1, 8)
vals = seq(1, 10, 1)
d = expand.grid(var1=vals, var2=vals)

d$cov = pmap_dbl(d, function(var1, var2){
  # Generate random variables and check covariance between x1 and x3 (= x1 + x2)
  n = 1000
  x1 = rnorm(n, 0, sqrt(var1))
  x2 = rnorm(n, 0, sqrt(var2))
  x3 = x1 + x2
  cov(x1, x3)
})
ggplot(d, aes(var1, cov, color=factor(var2))) +
  geom_path() +
  geom_abline(intercept=0, slope=1, size=1, linetype='dashed') +
  labs(x='Var(x1)', y='Cov(x1, x2)', color='Var(x2)') +
  coord_equal()

jet.colors <- c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000")
ggplot(d, aes(var1, var2, fill=cov)) +
  geom_tile() +
  scale_fill_gradientn(colors=jet.colors) +
  labs(x='Var(x1)', fill='Cov(x1, x2)', y='Var(x2)') +
  coord_equal()



...but not quite right either, presumably as I need to incorporate $\sigma^2_{\epsilon}$ into the equation somewhere.
Is there a known analytic solution to this?
 A: Suppose $x_1$ and $\epsilon$ are independent, then 
$$\mathrm{Cov}\left(\matrix{x_1\\ \epsilon}\right) = \left (\matrix{\sigma_1^2 & 0 \\0&\sigma_\epsilon^2}\right)$$
$$\left(\matrix{x_1\\ x_2}\right) = \left (\matrix{1 & 0 \\1 &1}\right)\left(\matrix{x_1\\ \epsilon}\right) $$
So 
$$\mathrm{Cov}\left(\matrix{x_1\\ x_2}\right) = \left (\matrix{1 & 0 \\1 &1}\right)\left (\matrix{\sigma_1^2 & 0 \\0&\sigma_\epsilon^2}\right)\left (\matrix{1 & 1 \\0 &1}\right) = \left (\matrix{\sigma_1^2 & \sigma_1^2 \\\sigma_1^2 &\sigma_1^2+\sigma_\epsilon^2}\right)$$
So the correlation coefficient $ρ = \frac {\sigma_1^2}{\sqrt{\sigma_1^2 (\sigma_1^2+\sigma_\epsilon^2)}} = \frac {\sigma_1}{\sqrt{\sigma_1^2+\sigma_\epsilon^2}}$
A: Thank you to @a_statistician for the correct analytical solution.
I've since realised why my simulation results haven't agreed with this analytic result (cov(x1, x2) == var(x1)): the effect of sampling error is surprisingly strong here, even with large N, and so the correlation between $x_1$ and $\epsilon$ weren't necessarily 0, and the sample variances of $x_1$ and $\epsilon$ weren't necessarily their population values of $\sigma^2_1$ and $\sigma^2_{\epsilon}$.
Here's the corrected simulation code.
library(tidyverse)
orthogonalise = function(x1, x2){
  X = data.frame(x1, x2) %>% as.matrix
  pca = prcomp(X)
  pX = X %*% pca$rotation
  return(list(x1=pX[,1], x2=pX[,2]))
}
exact.orthogonal.samples = function(var1, var2){
  n = 1000
  x1 = rnorm(n, 0, 1)
  x2 = rnorm(n, 0, 1)
  # Remove correlation
  ox = orthogonalise(x1, x2)
  x1 = ox$x1
  x2 = ox$x2
  # Set variance
  x1 = x1 * (sqrt(var1)/sd(x1))
  x2 = x2 * (sqrt(var2)/sd(x2))
  return(list(x1=x1, x2=x2))
}

vals = 2^seq(1, 8)
vals = seq(1, 10, 1)
d = expand.grid(var1=vals, var2=vals)
d$cov = res = pmap_dbl(d, function(var1, var2){
  # Generate random variables and check covariance between x1 and x2 (= x1 + e)
  X = exact.orthogonal.samples(var1, var2)
  x1 = X$x1
  e = X$x2
  x2 = x1 + e
  cov(x1, x2)
})


jet.colors <- c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000")
ggplot(d, aes(var1, var2, fill=cov)) +
  geom_tile() +
  scale_fill_gradientn(colors=jet.colors) +
  labs(x='Var(x1)', fill='Cov(x1, x2)', y='Var(x2)') +
  coord_equal()


