Error on mean from measurements made from a distribution with a possible long tail I am using Langevin dynamics simulations to measure a first hitting time of particles leaving a region in parameter space. I would like to measure the mean first hitting time (and not the median for example) because after making some simplifying assumptions, I can calculate the mean of the first hitting time explicitly (but not the whole distribution), and want to compare this with data to test whether those simplifying assumptions put me in the right ball park. 
My problem is that the mean first hitting time measured from my simulation seems to increase with the time that I run for (more time means more trajectories/samples used in calculating this mean). Between a run of length $T$ and a run of length $\sim 1.2T$, the mean increased by more than 3 standard errors. A similar increase exists between a run of length $\sim 0.8T$ and $T$. So the tail of the empirical distribution seems to affect the mean, and there seems to be a systematic underestimation of both the mean and the standard error, which made me wonder about a long tailed distribution. However theoretically it is impossible for a set of non-zero measure of trajectories to remain inside the region forever, therefore the mean time of exit must be finite.
The histogram of my current samples looks like this, the average is ~1.4:

I am limited in computational power as to how long I can run my simulation for. I expect that if I were to continue to run my simulations, the mean would eventually asymptotically reach some value, however this is not yet visible from my current samples. Given the data that I have, is there a way to generate the confidence intervals on my current mean to ensure that this asymptotic value is included in the interval?
A few more details about the simulation/system: The region the particles are escaping from is inhomogenous and has other occupants, so theoretically it is a black box where not much can be done theoretically other than the few simplifying assumptions I am trying to test out. Multiple particles are escaping at the same time, however they interact minimally with each other and their density inside the region is low, so low chances of encountering another particle. Therefore the exit times of two particles is not exactly independent, but any correlations between the times of two particles will be extremely weak. The system is also at steady state, so on average there is the same number of particles inside the region for all times. 
 A: The difference of the empirical means of two sets of samples might well be larger than the standard error of the mean. Asymptotically, the standard error defines a confidence interval to the level 0.67 or so.
If the distribution has long tails and the sample mean is unstable, the sample standard deviation should be large. As long as the underlying distribution has a well-defined mean and variance, things should work out more or less.
Finally, if the sample mean is unstable, maybe it makes more sense to estimate and report a more robust measure of location, such as the median.
A: I think the answer to your problem may be more related to your field of study (say, molecular dynamics) than statistics. Statistics deals with probability distribution of random variables and, it does not matter what your distribution is, for applying central limit theorem, for making inferences about population (with large samples). However, since you have mentioned that such simulation does not work in your case with mean changing randomly every time, it is quite possible that the 'time spent in the region' may be a random variable with no finite mean at all (for example, a random normal distribution assumes that, even though the variable is random, the variable will tend towards a central value i.e. mean, so is the case with other known standard distributions which mimic some natural phenomenon). But if the 'time spent' variable is completely random, whose value is affected by other variables in your model, you cannot expect 'any mean value' whatever time you spend in simulation. For example, movement of stock prices over a time interval are completely random in nature, and it cannot be modeled under any known probability distribution and it cannot have any mean related to it. 
If your theory, however, indicates that there should be finite mean for time spent (which I intuit very unlikely for molecular movements) and if you are sure about it, then you should work on finding out the distribution function and its moment generating functions etc., to arrive at the mean function. 
Note: If finite mean is a sure theoretical possibility, another statistically sound method to estimate the mean is MLE estimation methods (especially in data science paradigm). But, MLE methods, that estimate population parameters, mandatorily require an assumption about the nature of population distribution and its density function (which is not required for CLT based inferential statistics that makes inferences on sample). But there could be other techniques available in data science paradigm, which could estimate population parameters, but they may rely on quantitative and computationally intensive algorithms.
A: The behavior that the mean keeps increasing when more samples are added, typically means that the mean time is infinite. You could for example be facing a pareto distribution with $\alpha\le 1$. 
You may not be facing this exact distribution, but it is the most standard example of a fairly simple distribution with infinite mean, so it is worth taking a look at; and it will have the same behavior as you are describing: it only takes finite values (mostly 1 or slightly above) so the mean of any finite sample will be finite, but it will grow as the number of samples grows.
To fortify my argument, I will explain why your "I expect that if I were to continue to run my simulations, the mean would eventually asymptotically reach some value" must be wrong. Let's say your mean after the first $n$ simulations is $\mu_1$, and after the first $2n$ simulations it is $\mu_2<\mu_1$. Then, the mean among simulations $n+1,\ldots,2n$ was $\mu_3=2\mu_2-\mu_1<\mu_2$, meaning we have $\mu_3<\mu_2<\mu_1$. However, since your simulations are independent, $E[\mu_3]=E[\mu_1]$, hence you cannot consistently find this behavior if the mean was finite and well-defined.
A: EDIT:
If the mean is changing significantly when you add more data then you may not have a finite mean or not unimodal. hist(c(rnorm(100, 0, 1) + rnorm(100, 3, 1)))`. But I'm guessing you've already made a histogram of your data or have been able to visualize it. So Here's what I would do. 1.) Calculate the mean with an increasing number of data samples. The more you can do the better but only do what is feasible. 2.) Look to see if the means follow some type of infinite series that converge. Maybe a geometric series, p-series or something like that. If you can find that the mean follows some kind of pattern like $const.\frac{n}{n+1}$ would be a simple one. Even if you're mean is moving up and down and it doesn't seem to follow a pattern you may have $-1^{i}$ in the sequence or something similar. 
I hope this is a little closer to what you were looking for.  
When you say simulations I'm assuming that you are either running Monte Carlo simulations or constructing bootstrap samples. Either way here is what you want to do. 
1. run b simulations
2. during each bth simulation draw a random sample that is the size of the data set with replacement. 
3. Calculate the mean of the sample drawn. 
4. The simulation is complete so now you have an array with b sample means. Take the mean of that array and this is your mean. You can also create credible regions (when you run a simulation confidence intervals no longer apply because your regions are empirical). 
B <- 1000
n <- dim(data)[1]
ind <- seq_len(n)
means_vector <- double(B)
for ( b in seq_len(B)) {
    rand_samp <- sample(ind, n, replace = TRUE)
means_vector[B] <- mean(data[rand_samp])
}
mean(means_vector)
quantile(means_vector, c(0.25, 0.975))

This is the bootstrap method. It's pretty simple, and there's an R function "boot()" that will make things crazy fast.
