For what it’s worth, it seems that, after the two changes mentioned above (swapping the lines assigning to mu[j] and dL0_star[j] and adding lower bound to vector<lower=0>[T] dL0;), the problem with the initialization is with the gradient: Gradient evaluated at the initial value is not finite. I think this is caused by some of the Idt values being zeroes because some of the values in Y are zero, while the parameter of the Poisson distribution should be positive. Without actually thinking about what the model is doing, I don’t know what would be the correct solution to this problem (or if this even is the real issue).
(edit: ok, one final guess without understanding details of the model: it seems that Poisson with parameter -> 0 tends to point mass at 0, so based on this one could skip all likelihood terms where Y[i, j] and dN[i, j] are both zero.)