optimization fails for a constrained parameter

Consider the following simple model, which should return the average of two positive numbers:

data {
  vector<lower=0>[2] X;
  real<lower=0> sigma;
}
parameters {
  real<lower=0> mu;
}
model {
  X ~ normal(mu, sigma);
}

R code:

avg_model <- stan_model("avg.stan")
optimizing(avg_model, list(X=c(0.3, 0.5), sigma=1e-4), init=0)$par

If sigma is not too low (1 or 0.1), it works fine. It also works fine if I remove the lower bound on mu. But if the lower bound is present and sigma is low (1e-3, 1e-4), optimizing() almost always returns mu = 0.

The manually transformed model has the same behavior:

data {
  vector<lower=0>[2] X;
  real<lower=0> sigma;
}
parameters {
  real log_mu;
}
transformed parameters {
  real mu = exp(log_mu);
}
model {
  X ~ normal(mu, sigma);
  target += log_mu;
}

Yet the likelihood looks fine, and R's optim() has no problem maximizing it despite using the approximate gradient.

liky <- Vectorize(function(log_mu) {
  X <- c(0.3,0.5)
  sum(dnorm(X, mean=exp(log_mu) , sd=1e-4, log=T)) + log_mu
})
curve(liky(x), from=-10, to=-0.2)
optim(0, liky, control = list(fnscale=-1))$par

rstan version 2.16.2

Start by reproducing the constrained model through R's optimizing() entry point with sigma values from 1 down to 1e-4, then compare it with the manually transformed model and R's optim() result. Done means identifying why low-sigma constrained optimization returns mu = 0 and making the constrained case reliably recover the average of the two observations.