This paper discusses how to compute oscillon lifet...

I would not use a single “mean oscillon” as the main lifetime estimator. It is fine as a diagnostic plot or sanity check, but it is likely to give a biased lifetime.

The reason is that the map

$$(A, R, E, \text{profile shape}, M) \longrightarrow \tau_{\rm osc}$$

is highly nonlinear. In the paper, the authors explicitly find that lifetimes depend strongly on the initial oscillon shape and potential parameter, and they therefore scan the amplitude–radius plane rather than relying on one averaged object. They also show that oscillons with the same amplitude but slightly different radii can have lifetimes differing by about an order of magnitude; for example, at fixed $A_-^{(e)}=0.46$, changing $R_p^{(e)}$ from 3.25 to 3.65 changes the decay time from roughly $5\times10^4$ to $5.2\times10^5$ in units of $m_\phi t$. (1907.00611v2.pdf) (1907.00611v2.pdf)

Averaging $A$ and $R$ separately is especially dangerous because oscillon amplitudes and radii are correlated, often anti-correlated, and the breathing mode makes the instantaneous $A,R$ values phase-dependent. The paper notes that $A(t)$ and $R(t)$ oscillate oppositely during the breathing mode, and that the breathing amplitude itself depends on the particular oscillon shape. (1907.00611v2.pdf) A mean oscillon can therefore land in a part of parameter space that is not representative of the actual ensemble, or even artificially close to or far from the stability boundary.

A better manageable method is an ensemble-weighted representative sampling / surrogate approach.

For each value of your monodromy parameter $M$, I would do the following.

First, extract a catalogue of oscillons from the 3+1D simulation at a fixed post-stabilization time or scale factor. For each oscillon record not only radius and amplitude, but also energy, central frequency if available, breathing amplitude, asphericity, and nearest-neighbor distance. The paper’s extraction strategy is basically of this type: identify oscillons in 3D, fit them by Gaussian-like radial profiles, and use the fitted amplitude and radius as inputs to radial simulations. (1907.00611v2.pdf)

Second, do not simulate all oscillons individually unless the number is modest. Instead, cluster the oscillons in descriptor space,

$$\theta_i(M)= \left( A_i,\; m_\phi R_i,\; E_i,\; B_i,\; \epsilon_i,\; d_i \right),$$

where $B_i$ is a breathing-mode measure, $\epsilon_i$ an asphericity measure, and $d_i$ a nearest-neighbor or isolation measure. Then choose actual medoids or stratified representatives, not arithmetic averages. For example, choose representatives from the median, the high-energy tail, near the stability boundary, and the sparse tails of the $A$-$R$ distribution.

Third, run 1+1D radial simulations only for those representatives. Use their weights, meaning the number or energy fraction of 3D oscillons represented by each cluster. This gives you an ensemble lifetime distribution rather than one lifetime:

$$P_M(\tau) \approx \sum_{k=1}^{K} w_k(M)\, \delta\!\left[\tau-\tau_k(M)\right].$$

For cosmological impact, the more useful observable is often not the mean lifetime but the survival function,

$$S_N(t|M)= \frac{\sum_k w_k\,\mathbf{1}[\tau_k>t]} {\sum_k w_k},$$

and the energy-weighted survival fraction,

$$S_E(t|M)= \frac{\sum_k w_k E_k\,\mathbf{1}[\tau_k>t]} {\sum_k w_k E_k}.$$

That directly tells you how much of the oscillon population, or energy density, survives to a given time.

Fourth, build a surrogate lifetime map,

$$\tau_{\rm osc} = f(M, A, m_\phi R, E, B,\ldots),$$

using the representative radial simulations. A simple interpolator, Gaussian-process emulator, radial-basis interpolation, or even piecewise interpolation over the $A$-$R$-$M$ grid is enough. The literature already shows that the lifetime landscape has a stability/instability boundary in the $A$-$R$ plane, so you should sample more densely near that boundary. (1907.00611v2.pdf)

For your monodromy inflation study, I would separate the question into two parts:

$$M \longrightarrow \left[ N_{\rm osc}(M),\; P_M(A,R,E,B,\epsilon),\; f_{\rm osc}(M) \right]$$

from the 3+1D simulations, and

$$(M,A,R,E,B,\epsilon) \longrightarrow \tau_{\rm osc}$$

from the radial simulations. Then combine them. This will let you say not only “larger/smaller $M$ produces more oscillons,” but also whether those oscillons are typically longer-lived, shorter-lived, closer to the instability boundary, more energetic, or more strongly breathing.

A minimal version that is still scientifically defensible would be:

1. For each $M$, run several 3+1D seeds and extract all oscillons.
2. Plot the $A$-$R$-$E$ distribution.
3. Choose perhaps 10–30 representative oscillons per $M$: median, quartiles, high-energy tail, low-energy tail, near-boundary cases, and a few random checks.
4. Run radial simulations for those representatives.
5. Report median lifetime, energy-weighted mean lifetime, 10–90% lifetime range, and survival fraction $S_E(t|M)$.
6. Include the lifetime of the “mean oscillon” only as a comparison curve, not as your primary result.

This is better than a single mean oscillon because it preserves the finite-lattice ensemble information without requiring a radial simulation for every object.

One practical improvement over the paper’s approach would be to use the actual angle-averaged radial profiles of representative oscillons, not only Gaussian fits. The paper notes that the Gaussian approximation works well for their hilltop case but also acknowledges that real oscillons are not perfectly Gaussian or perfectly spherical. (1907.00611v2.pdf) (1907.00611v2.pdf) For monodromy potentials, the profile family may differ from the hilltop case, so keeping the actual radial profile for the cluster medoids is safer.

For comparing different $M$, normalize carefully. Use dimensionless variables such as

$$\tilde R = m_\phi(M) R_{\rm phys}, \qquad \tilde t = m_\phi(M) t,$$

but also convert lifetimes back to physical time or e-folds. The paper makes a similar point when comparing different hilltop parameters: lifetimes in $m_\phi t$ can change strongly, but $m_\phi$ itself also depends on the model parameters, so they also translate to e-folds. (1907.00611v2.pdf)

So my recommendation is: do not replace the ensemble by one averaged oscillon. Replace the ensemble by a weighted, representative lifetime distribution. That keeps the complexity manageable while still respecting the physics that lifetime is shape-sensitive, phase-sensitive, and potential-dependent.