What is the distribution of $e=Y-\mathbb{E}(Y)$ where $Y=\exp(u), \ \ \ u\sim\mathbb{N}\left(\mu,\sigma^2\right)$ As $Y$ is log-normal we've $Y\sim \mathbb{LN}\big(\exp(\mu+\sigma^2/2),\exp(2\mu+\sigma)(\exp(\mu^2)-1)\big)$.
Now I define $e = Y - \mathbb{E}(Y) = Y - \exp(\mu+\sigma^2/2)$. 
As $e$ is the centered version of $Y$ I would expect that $e$ has a first moment of zero. I'm not sure about the distribution because $e$ can obtain negative values so log-normal is out of the question. Looking at some empirical results it seems that normal is not a proper fit either (top left is the histogram of $\hat{u}=\log(Y)-\widehat{\mu}$, top right is $\exp(\hat{u})$, bottom left corresponds to $\hat{e} = Y - \exp(\hat{\mu}+s_{\hat{u}}^2/2)$. The red line gives the log-normal $\mathbb{LN}(\exp(\frac{0.5^2}{2}),\exp(0.5^2)(\exp(0.5^2)-1))$ and the blue line the normal $\mathbb{N}(\bar{\hat{e}}=0.36,s_{\hat{e}}^2=2.67)$. We see that neither does fit the observed distribution properly:

library(data.table)
set.seed(100)
n <- 100000
mu <- 1
sigma <- 0.5
e <- rnorm(n,mu,sigma)
y <- exp(e)
data <- data.table(y=y,e=e)
mod <- lm(log(y)~1)
summary(mod)
data[,u:=resid(mod)]
s <- data[,sqrt(mean(u^2))]
data[,y_hat_exp:=exp(fitted(mod)+(s^2/2))]
data[,e:=y-y_hat_exp]

par(mfrow=c(2,2))
u <- data[,u]
mean.u <- mean(u)
sd.u <- sd(u)
h1 <- hist(u,freq=FALSE,ylim=c(0,1),main="u = log(y)-y_hat")
xfit<-seq(min(u),max(u),length=40) 
yfit<-dnorm(xfit,mean=mean.u,sd=sd.u) 
lines(xfit, yfit, col="red", lwd=2)

u.exp <- data[,exp(u)]
h1 <- hist(u.exp,freq=FALSE,ylim=c(0,1),main="exp(u)")
xfit<-seq(min(u),max(u),length=40) 
yfit<-dlnorm(xfit,mean=mean.u,sd=sd.u) 
lines(xfit, yfit, col="red", lwd=2)

e <- data[,e]
h2 <- hist(e,freq=FALSE,ylim=c(0,0.5),main="e = y-exp(mu-s^2/2)")
xfit<-seq(min(e),max(e),length=40) 
mean.e <- exp(mean.u+sd.u^2/2)
var.e <- exp(2*mean.u+sd.u^2)*(exp(sd.u^2)-1)
sd.e <- sqrt(sd.e)
yfit<-dlnorm(xfit,mean=mean.e,sd=sd.e) 
mean.e <- mean(e)
sd.e <- sd(e)
yfit2<-dnorm(xfit,mean=mean.e,sd=sd.e) 
lines(xfit, yfit, col="red", lwd=2)
lines(xfit, yfit2, col="blue", lwd=2)

 A: One can derive that $e = Y - \mathbb{E}(Y)$ is log-normal. To see this, define $\mathbb{E}(Y):=\mu$, $Var(Y) := \sigma^2$ and simply write out the probability density functions (pdfs):
\begin{align}
f_Y(ln(x)) &= \frac{1}{\sigma \sqrt{2\pi}}e^
\frac{- \left( {x - \mu } \right)^2}{2 \sigma^2} \\
f_e(x) &= f_Y(x - \mu)
\end{align}
Now $f_e(ln(x)) = f_Y(ln(x - \mu))  = \frac{1}{\sigma \sqrt{2\pi}}e^
\frac{- \left( {(x - \mu) - \mu } \right)^2}{2 \sigma^2} = \frac{1}{\sigma \sqrt{2\pi}}e^
\frac{- \left( {x } \right)^2}{2 \sigma^2}$~$N(0,\sigma^2)$.
So the logarithm of $e = Y-\mu$ is standard normal, i.e. $ln(e)$~$N(0,\sigma^2)$.
A: The answer is simply
$$
f_E(e) = \frac{1}{ [e+\exp(\mu+\sigma^2/2)]\sigma \sqrt{2 \pi}}\exp\left[-\frac {(\mbox{ln}[e+\exp(\mu+\sigma^2/2)] - \mu)^{2}} {2\sigma^{2}}\right]
$$
This can be proven using the "change of variables" property
$$
f_E(e) = \left| \frac{d}{de} (g^{-1}(e)) \right| \cdot f_Y(g^{-1}(e))
$$   
Where $f_Y$ is the log normal density function
$$
f_Y(y) = \frac{1}{ y\sigma \sqrt{2 \pi}}\exp\left[-\frac {(\mbox{ln}y - \mu)^{2}} {2\sigma^{2}}\right]
$$
and $g$ is the monotonic function defined by
$$
g(y) = y-\exp(\mu+\sigma^2/2)
$$
Note $g^{-1}(e) = e+\exp(\mu+\sigma^2/2)$ and so $ \frac{d}{de} (g^{-1}(e))=1$.  Thus $f_E(e) = f_Y(g^{-1}(e))$.
