What is the correct bayesian formulation for the zero-truncated Poisson lognormal model? In ecology we use compound distributions to describe species-abundance data. One example is the Poisson Lognormal (PLN) distribution which is a Poisson distribution with rate parameter $\lambda$ that follows a lognormal PDF $f(\lambda;\mu,\sigma)$. The total probability of an observation (count) $x$ is then given by:
$$ p(x;\mu,\sigma) = \int_{\lambda = 0}^{+\infty}{\frac{\lambda^x e^{-\lambda}}{x!} f(\lambda;\mu,\sigma) d\lambda} $$
We can't observe zeroes so, we work with a truncated distribution to then construct a likelihood function and use the MLE frequentist estimator to obtain point estimates of $\mu$ and $\sigma$. The truncated PLN distribution is:
$$ p_{T}(x;\mu,\sigma) = \frac{p(x;\mu,\sigma)}{p(x>0;\mu,\sigma)} = \frac{p(x;\mu,\sigma)}{1-p(0;\mu,\sigma)} $$
The conditioning is done based on the total probability (the total probability of $x>0$) not on the conditional probability of $x>0$ (conditioned on a known value of $\lambda$). How can I formulate this model in Bayesian? I would like to use Stan to sample from the posterior $p(\mu,\sigma|\mathbf{x})$ without defining a new distribution.

Stan natively supports truncated distributions such that you can express the zero truncated likelihood $P(x|\lambda)$ with a sampling statement like this: x ~ poisson(lambda) T[1,];.
The Poisson Lognormal (PLN) model that you described can be expressed such that there is a different Poisson rate parameter $\lambda_i$ corresponding to each observation $x_i$.
$$
\begin{aligned}
log(\lambda_i) \sim & \ \mathrm{Normal}(\mu, \sigma) \\
x_i \sim  & \ \mathrm{ZeroTruncatedPoisson}(\lambda_i) \\
\end{aligned}
$$
We also need to define priors on the model parameters $\mu, \sigma$. In practice, we usually have some information specifying a reasonable range for model parameters to specify an informative prior. However, in the absence of any prior information, we can specify uninformative priors on these parameters. For normal distribution parameters, we can use the improper Jeffrey's prior:
$$
\begin{aligned}
\mu \sim & \ \mathrm{Uniform}(-\infty, +\infty) \\
log(\sigma) \sim & \ \mathrm{Uniform}(-\infty, +\infty) \\
\end{aligned}
$$
We can express this model in a Stan program, without defining any custom distributions, as follows:
data {
  int<lower=0> N;    // number of observations
  int<lower=1> x[N]; // observations
}
parameters {
  real mu;
  real log_sigma;
  vector[N] log_lambda;
}
transformed parameters {
  vector[N] lambda = exp(log_lambda);
  real sigma = exp(log_sigma);
}
model {
  // implicit uninformative improper uniform priors on `mu` and `log_sigma`

  // lognormal prior on lambda
  log_lambda ~ normal(mu, sigma);

  for (i in 1:N) {
    // zero-truncated poisson likelihood for x
    x[i] ~ poisson(lambda[i]) T[1,];
  }
}

Then, using this Stan model, we can sample from the posterior distribution $P(\mu, \sigma|x)$. Here's an example where the "true" values are $\mu=1.5, \sigma=0.75$. The range samples from the posterior distribution appear to overlap well with the "true" values in this case.
import pystan
import scipy.stats
import seaborn as sns

N = 200
mu = 1.5
sigma = 0.75
lambda_ = np.exp(np.random.normal(size=N) * sigma + mu)
x = scipy.stats.poisson.rvs(lambda_)
x = x[x > 0]

data = dict(N=len(x), x=x)

model = pystan.StanModel(model_code=MODEL_CODE)
samples = model.sampling(data=data)

df = samples.to_dataframe()
sns.pairplot(df[["mu", "sigma"]], plot_kws=dict(alpha=0.05))


