The theoretical asymptotic bounds provided for many Monte Carlo methods cannot be calculated using known quantities during search.
Often search results provide no principled measure of confidence and may be sub-optimal due to premature termination.

In this work, we prove two sets of bounds on the error rate of Monte Carlo search over non-stationary bandits and Markov decision processes.
The presented bounds hold for general Monte Carlo solvers meeting mild convergence conditions.
These bounds can be directly computed at the conclusion of the search and do not require knowledge of the true action-value, allowing them to be used as search stopping criteria.
We also provide a simple sub-optimality probability estimation method based on the presented bounds.
We empirically test the tightness of the bounds and accuracy of the estimator through experiments on a multi-armed bandit and a discrete Markov decision process.

Monte Carlo, Tree Search, MCTS, Bandit, MDP, POMDP

§ INTRODUCTION

Monte Carlo (MC) estimation can be used to solve decision making problems when the dynamics of the environment are not fully known.
Solvers draw samples from a generative model of the environment to estimate the values of different actions.
Typically, solvers that use Monte Carlo estimation sample actions during simulation according to a non-stationary search policy.
In these cases, empirical error bounds cannot be computed using standard methods that assume independent identically distributed ($iid$) samples.
Because of this, many solvers return action value estimates without a confidence measure.

Solvers based on Monte Carlo sampling have been proposed for single-step and sequential decision problems [Sutton and Barto, 2018].
Online search methods such as Monte Carlo tree search (MCTS) have solved a variety of complex sequential decision making tasks [Silver et al., 2016, Silver et al., 2017], though they are computationally expensive to execute.
The expected error of a Monte Carlo estimate

tends to decrease with the number of samples drawn [Audibert et al., 2009].

When expected error cannot be computed during search, arbitrary, heuristic stopping conditions are often used.
This can lead to high computation cost from excessive sampling or premature termination and sub-optimal results.
Measurable bounds or estimates of error probability would provide a more useful stopping criteria.

Several Monte Carlo solvers have theoretical asymptotic bounds on measures such as regret <cit.> PAC-correctness [Kearns et al., 1999, Kaufmann and Koolen, 2017], and error probability [Shah et al., 2020].
These theoretical bounds are defined in terms of parameters not generally known during search.
For example, regret bounds for the upper confidence tree (UCT) MCTS algorithm are defined in terms of the difference between the converged Monte Carlo estimate and the true optimal discounted return [Kocsis and Szepesvári, 2006].
These theoretical expressions help to understand algorithm behavior, but they cannot be numerically calculated.

Empirical bounds may be calculated using only values observed from simulation and can be used to evaluate solver performance during search.
Empirical bounds have been proven for some MC solvers under special circumstances [Grill et al., 2016, Shah et al., 2020], but are not generally available.

§.§ Our Contributions

In this work, we present two sets of finite-sample, empirical bounds on error probability for general Monte Carlo solvers.
One set of bounds can be applied to any Monte Carlo estimate of problem with bounded value.
The second set of bounds provides tighter limits for problems in which the central limit theorem (CLT) can be applied.

In both cases, we present bounds on the probability that the action value estimated by the search overestimates the true value by more than $\epsilon$ as

\begin{equation}
P(\bar{Q}(s, a) \geq Q(s, a) + \epsilon) \leq \alpha
\end{equation}

where $Q(s, a)$ is the expected sum of discounted returns from taking action $a$ from state $s$ and $\bar{Q}(s, a)$ is the action value estimate from the MC search.
We also bound the probability that the action returned from search $a$ is worse than any other considered action $a'$ by more than $\epsilon$ as

\begin{equation}
P(Q(s, a) \leq Q(s, a') - \epsilon) \leq \alpha
\end{equation}

when following the search policy for the remaining trajectory.

Both sets of bounds can be numerically calculated during search for use as a confidence measure or as a stopping criteria.
We also present terms to estimate both of the presented error probabilities without upper-bound guarantees.

We evaluated the tightness of the bounds and the accuracy of the estimation through numerical experiments.
We tested performance on a multi-armed bandit and a discrete action space navigation task.
A simple, fixed-depth Monte Carlo solver was used for evaluation on both tasks.
Additionally, UCT-MCTS was tested on the sequential navigation task.
The general Monte Carlo bounds hold for all tested configurations, and the CLT bounds hold when sufficient samples are drawn.
Bound gaps tended to decrease with increased sample set size.
The error estimate tends to produce a closer approximation of the actual error rates than the bounds.

§.§ Related Work

Multi-armed bandit problems have been extensively studied by the statistics and learning communities [Lai and Robbins, 1985].
Prior work in statistical analysis of random processes provide foundational measures toward generally applicable non-asymptotic bounds.
Prior works develop empirical Bernstein-type bounds based on sample variance [Maurer and Pontil, 2009].
Early work proved bounds for sums of independent, identically distributed random random variables based on sample mean and variance [Beygelzimer et al., 2011].
Further work extended these ideas to functions over martingales [Peel et al., 2013, Zhang et al., 2021] with some restrictions.

Many algorithms have been proposed to approximately solve the bandit problem.
Foundational algorithms, such as UCB1, have been extensively studied [Auer et al., 2002].
Similarly, many MCTS methods have been proposed to solve MDPs and partially observable MDPs (POMDPs).
UCT is a well-studied algorithm [Kocsis and Szepesvári, 2006] that has been used as the basis for several variants [Silver and Veness, 2010].
Convergence guarantees and non-asymptotic performance bounds have been proved for UCT and other Monte Carlo methods.
These bounds are often theoretical in that they are defined in terms of the unknown parameters such as the true expected return value.

Few algorithms provide performance bounds that can be numerically calculated using only empirically observed values.

Some solvers have been proposed with calculable non-asymptotic bounds.
Shah et al., 2020 present non-asymptotic analysis of UCT with improved regret-bounds, and additionally propose a modified, fixed-depth tree search with policy-improvement guarantees.
Early work on approximate sampling by Kearns et al., 1999 proved finite-sample convergence bounds for fixed-depth search, with limited scalability to large problem spaces.
An MCTS approach to simple regret minimization with finite-sample guarantees was proposed by Feldman and Domshlak, 2014.
Several other MCTS approaches have been proposed that provide finite-sample error guarantees [Grill et al., 2016, Kaufmann and Koolen, 2017, Grill et al., 2019].
While the bounds in these works hold for their respective algorithms, they are not generally applicable to other search methods.

§ PRELIMINARIES

Monte Carlo search estimates the values of taking actions from a given state.
The value of taking action $a$ in state $s$ and then following policy $\pi$ is denoted $Q^\pi(s,a)$.
Values associated with the optimal policy are denoted $Q^*(s,a)$.
Monte Carlo search is used as a policy by repeatedly searching from each encountered state and taking the action with the highest estimated value.
We will denote values associated with a Monte Carlo policy with superscript $\pi$.
For time-discounted problems, the time-discount rate is denoted $\gamma$.

The set of all available actions is denoted $\mathcal{A}$ and the set of all actions included in a Monte Carlo search is $\mathcal{A}_T \subseteq \mathcal{A}$.
For ab action $a_i$ the value is estimated using the sample set $\mathcal{C}_i \gets (q_i^{(1)}, \dots, q_i^{(n)})$.
The value estimate of action $a_i$ is then $\bar{Q}^{(n)}_i = \frac{1}{n_i}\sum_{q \in \mathcal{C}_i} q$, where $n_i = |\mathcal{C}_i|$.
The notation $\bar{X}^{(n)}$ denotes the Monte Carlo estimate of true quantity $X$ with $n$ samples.
The sample variance of the estimator is denoted $V_i$.

The bounds presented in this work may be applied to Monte Carlo search over sequential decision problems with non-stationary dynamics.
This includes multi-armed bandits (MABs) with non-stationary payout distributions and Markov decision processes (MDPs).
The bounds may be applied to Monte Carlo estimates over problems satisfying the following conditions almost surely.

* Bounded Returns: The sum of discounted rewards is bounded by some non-zero constant $a$ almost surely as $|\bar{Q}^{(n)}_i| \leq b \ \forall \ i, n$.

* Convergent Search: The estimated action values converge to the value of the search policy with infinite samples as $\mathbb{E}[\bar{Q}(s, a)] = Q^\pi(s,a)$.

* Accurate Simulation: The generative model used for simulation exactly matches the true environment dynamics.

For MDPs, we additionally assume that the search simulates trajectories from a given starting state $s$ to some potentially variable depth $h$ and uses a given estimator for the value the remaining trajectory.
We refer to the bounds that only require satisfying these conditions as general bounds.

We propose an additional set of bounds for problems to which the central limit theorem applies.
We refer to these bounds as CLT bounds.
In addition to the above conditions, application of these bounds assumes that distributions of MC estimates approaches a normal distribution with infinite samples as

\begin{equation}
\sqrt{n}\big(\bar{Q}(s,a) - Q(s,a)\big) \xrightarrow{d} \mathcal{N}(0, \sigma^2)
\end{equation}

where $\sigma^2$ is the variance of returns under the optimal policy.
For samples generated under a stationary policy, the CLT should hold.
Early works by Markov, 1951, Bernstein, 1927, and Dobrushin, 1956 proved that the CLT holds for payoffs from non-stationary Markov chains meeting certain ergodic conditions.
Proving the CLT for additional cases is an active area of research [Sinn and Chen, 2013, Arlotto and Steele, 2016].
These may be extended to return estimates under non-stationary policies, however, this is left for future work.

§ EMPIRICAL BOUNDS AND ESTIMATOR

We present two sets of error bounds for Monte Carlo search.
The first set of general bounds apply to search over any problem meeting the conditions in <ref>.
The second set of bounds are only valid for problems in which the central limit theorem holds.
We also propose terms to estimate the errors presented without enforcing an upper bound.

§.§ General Monte Carlo Bound

We present bounds on two types of errors that may result from a Monte Carlo search returning an action $a$ with estimated optimal value $\bar{Q}^*(s,a)$.
The first is a bound on the probability that a value is over-estimated by the Monte Carlo value.
We refer to this as an value error bound.
The second is a bound on the probability that an action $a' \in \mathcal{A}_T/a$ is better than the action recommended by search.
We refer to this as a bound on the action error probability.

The value error bound gives the probability that an estimated action value is over-estimates the optimal action value.
It is stated formally in <ref>.

For a Monte Carlo search returning action value estimate $\bar{Q}_i = q_i$, the probability that the true optimal action value is overestimated $\bar{Q}_i - Q^*_i \geq \epsilon$ is no greater than

\begin{equation}~\label{eq: sub bound}
\exp\Big(\frac{-n_i(|\epsilon - \zeta_i|^+)^2}{2\hat{\sigma}^{2, \alpha}_i}\Big) + \alpha
\end{equation}

where $n_i$ is the number of sampled trajectories in $\mathcal{C}_i$.
An upper bound on the return variance is given by $\hat{\sigma}^{2, \alpha}$ with significance $\alpha \in (0, 1)$.
The $\zeta$ term is an upper bound on the expected bias of the baseline value estimate.
Here, $|x|^+$ represents the $\max(0, x)$ function.

The variance upper bound in <ref> is defined for bounded random variables as

\begin{equation}~\label{eq: var bound}
\hat{\sigma}^{2, \alpha} \gets b^2\Bigg(\frac{1}{b}V_n^\frac{1}{2} + \sqrt{\frac{-\ln{\alpha}}{n - 1}}  \Bigg)^2
\end{equation}

where $V_n$ is the sample variance over $n$ sampled trajectory returns and $b$ is the range of the sample values.
Taking the minimum of <ref> over $\alpha$ gives the tightest bound.
In the bandit setting, the bias term $\zeta = 0$, since the values are sampled directly from the true distributions.
In the MDP setting, the $\zeta$ term is an upper bound on the bias, defined as

\begin{equation}
\zeta = \frac{1}{|\mathcal{C}|}\sum_{k \in \mathcal{C}}\gamma^{H_k} \epsilon^U_k
\end{equation}

where $H_k$ is the depth of the Monte Carlo rollout, $\gamma \in [0, 1]$ is a discount factor, and $\epsilon_k^U$ is an upper bound on the error of the baseline value estimate at the end of the trajectory.
At worst, $\epsilon_k^U \gets \max(Q_{max} - \hat{Q}_k, \hat{Q}_k - Q_{min})$, where $\hat{Q}_k$ is the value of the baseline estimate at the end of the trajectory $k$.
The error is zero for trajectories reaching termination.

The presented bound is entirely empirical and requires no unknown values to compute.
It can be used as a stopping criteria to ensure a policy meets a return target with high probability.
As can be seen, the probability decays exponentially with the number of trajectories sampled for the given action.
For MDPs, the search-bias decreases with the depth $H$ of each trajectory sample as $\gamma^H$.

The action error bound limits the probability that an action from the search set is better than the recommended action when the Monte Carlo search policy is followed for subsequent steps.
It is stated formally in <ref>.

Let a Monte Carlo search return an action $a_i$ with estimated value $\bar{Q}_i$.
For any other action $a_{j} \in \mathcal{A}_T/a_i$ with value estimate $\bar{Q}_j = \bar{Q}_i - \delta$, the probability that the true value of $a_j$ is $Q^\pi_j \geq Q^\pi_i + \epsilon$ is no greater than

\begin{equation}~\label{eq: gen ae bound}
\exp\Big(\frac{-n_i n_j(|\delta + \epsilon - \zeta_i - \zeta_j|^+)^2}{2(n_i\hat{\sigma}^{2, \alpha}_i + n_j\hat{\sigma}^{2, \alpha}_j)}\Big) + \alpha
\end{equation}

where $n_i$ and $n_j$ are the number of sampled trajectories in $\mathcal{C}_i$ and $\mathcal{C}_j$, respectively.
An upper bound on the return variance is given by $\hat{\sigma}^{2, \alpha}$ with significance $\alpha \in (0, 1)$.
The $\zeta$ terms are upper bounds on the expected bias of the baseline value estimate.

Intuitively, the action error provides a limit on how likely it is that an action from the search set is better than the action with the maximum estimated value.

Though these bounds decay exponentially with the number of samples, empirical experiments show they are fairly loose.
We therefore present an additional set of bounds and error approximations that hold more tightly for appropriate problems.

§.§ Central Limit Theorem Bound

The previously presented inequalities only assumed the values estimates were bounded.
Though widely applicable, the presented bounds were found to be fairly loose.

We develop tighter bounds by assuming the distribution of return estimates follows a normal distribution for cases in which the central limit theorem (CLT) applies.

The central limit theorem is valid for estimates over stationary distributions, but not necessarily for non-stationary process.
The CLT has been shown to hold for payouts from non-stationary Markov chains for a wide range of cases.
By applying the CLT with a finite sample penalty, we can derive tighter versions of the bounds in <ref> and <ref>.

For a Monte Carlo search satisfying the central limit theorem returning action value estimate $\bar{Q}_i$, the probability that the true optimal action value is overestimated $\bar{Q}_i - Q^*_i \geq \epsilon$ is no greater than

\begin{equation}~\label{eq: clt se bound}
1 - \Phi(\epsilon - \zeta_i; 0, \hat{\sigma}_i^{2, \alpha}/n_i) + \alpha +  \frac{0.4748}{n_i}
\end{equation}

where $\Phi(x; 0, \sigma^2)$ gives the cumulative distribution function for a zero-mean Gaussian with variance $\sigma^2$.

This bound is built on the Gaussian assumption of the CLT with two penalties for finite samples.
The first is the variance upper bound previously presented in <ref>.
The second penalty is introduce in the far right term.
This term is a correction for CLT convergence based on the Berry-Esseen theorem [Berry, 1941, Esseen, 1942] with the constant term proved by Shevtsova, 2011.

As before, we present the corollary bound on error probability.

Let a Monte Carlo solver generate estimates such that the central limit theorem holds.
For a search that returns action value estimates over two actions $\bar{Q}_i - \bar{Q}_j = \delta$, the probability that the true action value gap $Q_i - Q_j \leq -\epsilon$ is no greater than

\begin{equation}~\label{eq: clt ae bound}
1 - \Phi(\delta + \epsilon - \zeta_i - \zeta_j; 0, \bar{\sigma}^{2, \alpha}) + \alpha +  \frac{0.4748}{n_i + n_j}
\end{equation}

where the variance upper bound is given by $\bar{\sigma}^{2, \alpha} = \hat{\sigma}^{2,\alpha}_i/n_i + \hat{\sigma}^{2,\alpha}_j/n_j$.

§.§ Error Estimation

Instead of computing an upper bound, which may be loose, it may be beneficial to instead compute an estimate of the error probability.
An error estimate that holds to the true value more tightly than a bound may be more useful for certain applications, such as for a stopping criteria during anytime search.
We can model the distribution of mean estimates of samples drawn from an unknown distribution using Student's t-distribution [Gosset, 1908].

To approximate the value error, a similar estimate is defined by

\begin{equation}~\label{eq: t subopt est}
1 - T_{n-1}\Big(\frac{\epsilon - \zeta_i}{(\sqrt{V_i}/n_i)}\Big)
\end{equation}

which follows from the single-variable t-distribution with $n-1$ degrees of freedom.

We can similarly estimate the action error as

\begin{equation}~\label{eq: t ae est}
1 - T_{\nu}\Big(\frac{\delta + \epsilon - \zeta_i - \zeta_j}{(\sqrt{V_i}/n_i + \sqrt{V_j}/n_j)}\Big)
\end{equation}

where $T_{\nu}(x)$ is the cumulative distribution function for the t-distribution with $\nu$ degrees of freedom.
For estimates over the difference of two variables, we approximate the degrees of freedom as

\begin{equation}
\nu = \frac{\big(V_i/n_i + V_j/n_j\big)^2}{V_i^2/n_i^2 + V_j/n_j^2}
\end{equation}

following the method proposed by Welch, 1947.

§ NUMERICAL EVALUATION

Experiments were conducted to evaluate the tightness of the presented bounds and the accuracy of the error approximation.
We tested two problems: a simple multi-armed bandit and a grid-world MDP.
We calculated the bounds under various conditions using a simple Monte Carlo estimator for both problems.
For the MDP, we also tested MCTS with upper confidence tree (UCT) sampling.
All testing was done in the Julia language using the POMDPs.jl framework [Egorov et al., 2017].
The $\alpha$ were minimized with Nelder-Meade optimization [Nelder and Mead, 1965].

The simple Monte Carlo solver draws $n$ samples of actions from a discrete distribution over the complete action space.
Distribution weights for each action $w_i$ are calculated from the softmax over the value estimates as

\begin{equation}
w_i \gets \frac{\exp{\bar{Q}_i \tau^{-1}}}{\sum_i \exp{\bar{Q}_i \tau^{-1}}}
\end{equation}

where $\tau$ is a temperature parameter, which was set to 10 for experiments in this work.
Actions that have not yet been sampled have $q_i \gets \infty$.
In the MDP setting, after sampling the first action, the remaining actions are selected according to an $\epsilon$-greedy baseline policy, up to a specified depth $h$.
Source code for the experiments is provided in the supplementary materials.

§.§ Multi-Armed Bandit

The bandit problem had 10 arms, each with a Gaussian payout distribution.
The tests were run in episodes, with each episode having differently parameterized distributions.
The payout mean values were sampled from a uniform distribution on $[-1, 1]$.
Each payout distribution had the same variance $\sigma^2 = 0.5$ each episode.

The bandit expected payouts were estimated with the Monte Carlo solver for varying numbers of samples $n$.

For each episode, the error probability bound was calculated for the two actions with the highest estimated values.
The sub-optimality bound was calculated for the best action for each sampling level and variance.
The margin was $\epsilon = 0.1$ for both experiments.
We calculated the true rates of error and sub-optimality in excess of $\epsilon$ across episodes.
The average upper bound calculated for each trial and one standard error bounds are plotted along with the true rate curves in <ref>.
The x-axis plots the total number of samples for the the search, not necessarily the samples over the measured action.

Bandit experiment results. Each graph gives the calculated general upper bound, the CLT bound, the error rate estimation, and the true observed error rate for varying numbers of sample rollouts. Figure (a) gives the action error and figure (b) gives the value error.

As can be seen, both the error and sub-optimality bounds decay exponentially with the number of samples drawn.
The CLT bounds are tighter than the general bounds in both cases for all sampled points.
The error estimate is significantly closer than either bound, though it does under-estimate the observed rate a several points.

§.§ Grid World MDP

In the grid world MDP, an agent is required to navigate through a discrete, 2D world from a random initial position to a goal destination.
There are two goal locations generating positive rewards of 10 and 3, and two penalty locations giving negative rewards of -10 and -5.
The episode terminates when the agent reaches any of the four reward-generating locations.
Each step, the agent takes an action to attempt to move in one of the four cardinal directions.
The agent successfully moves in the intended direction with a probability of 0.7 and moves randomly otherwise.
The grid tested was $10 \times 10$.
The time discount was set to $\gamma = 0.95$.

§.§.§ Simple MC Solver

We used offline value iteration to learn a baseline policy and approximate optimal action values.
The offline solver terminated after an update step with a Bellman residual less than $10^{-6}$.
The simple Monte Carlo solver used an $\epsilon$-greedy baseline policy select actions after the initial step.
The bootstrap value estimates were generated by the baseline policy with additive uniform noise bounded between $-0.1$ and $0.1$.
We ran the tests with the solver searching to depths of 5, 10, and 25.
As with the bandit problem, we calculated the error and sub-optimality rates for the top performing actions.
The results for depth 10 search are shown in <ref>.

Grid world MC solver results. Each graph gives the calculated general Chernoff upper bound, the CLT bound, the error rate estimation, and the true observed error rate for varying numbers of sample rollouts. Figure (a) gives the action error and figure (b) gives the value error.

As with the bandit results, both the error and sub-optimality bounds decay exponentially with the number of samples drawn, and the CLT bounds hold more tightly.
At low samples, the bounds are tighter for the MDP than for the bandit problem, due to the higher true error rates.
Unlike the bandit problem, the t-distribution estimate does not underestimate the observed error rate at any point.
This is likely due to the fact that the t-distribution estimate does still introduce a conservative error assumption in the bounds used to calculate the leaf node bias $\zeta$ terms.
Using an estimate of leaf node error instead of an upper bound could result in more accurate estimates, though with more frequent under-approximation.

§.§.§ Monte Carlo Tree Search

We tested the tightness of the bounds on the popular MCTS-UCT algorithm, using the MCTS.jl implementation.
We tested the performance with max depth set to 25 steps for each search.
The UCT exploration constant was set to 1.0.
As with the simple solver, leaf node values were estimated by the baseline policy with additive uniform noise.
The results are shown in <ref>.

Gridworld MCTS solver results. Each graph gives the calculated general upper bound, the CLT bound, the error rate estimation, and the true observed error rate for varying numbers of sample rollouts. Figure (a) gives the action error and figure (b) gives the value error.

The same trends can be observed in the MCTS results as in the simple Monte Carlo search results, however, the bounds are significantly looser at high sample counts in the MCTS tests.
This is likely due to the influence of several shallow searches early in the search leading to higher magnitude $\zeta$ terms.
Unlike the simple solver, the UCT algorithm does not satisfy the convergence condition, as value estimates include results from early search results during exploration.
This leads to a well known underestimation problem, often corrected through post-search max-value backup.
This leads to the observed rate of over-estimation being zero throughout search.
Including early samples in the estimates also means that the samples are not from a stationary distribution, and that the CLT may not necessarily apply.
Despite this, the CLT bound does not underestimate the true value at any search setting.

These results suggest several adaptations to MCTS methods to improve expected error bounds.
For a given number of samples, tighter error bounds can be achieved when the sample counts of the compare actions are comparable.
For situations where accurate action selection from a discrete set is important, this suggests that increased exploration at the end of search may improve error performance.
Deeper searches tend to result in better estimates than shallow searches.
This suggests modifications to either increase search depth or to more heavily weight samples from longer trajectories in the value estimate.

§ CONCLUSIONS AND FUTURE WORK

In this work, we presented empirical bounds on the error probability for convergent Monte Carlo solvers on non-stationary bandit problems and Markov decision processes.
We presented bounds for general, bounded-value problems with convergent solver and tighter bounds applicable when the central limit theorem holds.
To provide a useful metric for anytime search, we also proposed a simple error approximation based on the t-distribution.
We demonstrated that these bounds hold for a multi-armed bandit and a discrete grid world MDP for both a simple Monte Carlo estimator and a more complex tree search.
Results showed that the bounds are loose for low sample numbers, but decay exponentially as samples increase.
The results also show that the CLT bounds tend to be tighter than the general bounds, and that the error estimate is a good approximation of the true rate in the studied problems.
We anticipate that these measures may be especially useful to methods that use an informed baseline such as offline policy improvement and performance verification of safety-critical systems.

The current work only tested these bounds for simple bandits and MDPs.
Future work shall explore their performance in other, more complex tasks such as MCTS for partially observable MDPs [Silver and Veness, 2010].
Even with a moderate number of samples, the bounds were sometimes loose.
Future work will look to specialize the general form presented in this work to produce tighter bounds for specific algorithms.
The current work also only considered arithmetic mean estimators.
Future work will extend the bounds to more general estimators.

[Aras et al., 2019]
Efe Aras, Kuan-Yun Lee, Ashwin Pananjady, and Thomas A. Courtade.
A family of Bayesian Cramér-Rao bounds, and consequences
for log-concave priors.
In IEEE International Symposium on Information Theory
(ISIT), pages 2699–2703, 2019.

[Arlotto and Steele, 2016]
Alessandro Arlotto and J. Michael Steele.
A central limit theorem for temporally nonhomogenous Markov chains
with applications to dynamic programming.
Mathematics of Operations Research, 410 (4):0
1448–1468, 2016.

[Audibert et al., 2009]
Jean-Yves Audibert, Rémi Munos, and Csaba Szepesvári.
Exploration-exploitation tradeoff using variance estimates in
multi-armed bandits.
Theoretical Computational Science, 4100 (19):0
1876–1902, 2009.

[Auer et al., 2002]
Peter Auer, Nicolò Cesa-Bianchi, and Paul Fischer.
Finite-time analysis of the multiarmed bandit problem.
Journal of Machine Learning Research, 470
(2–3):0 235–256, 2002.

[Bernstein, 1927]
Serge Bernstein.
Sur l'extension du théoréme limite du calcul des probabilités aux
sommes de quantités dépendantes.
Mathematische Annalen, 970 (1):0 1432–1807,

[Berry, 1941]
Andrew C Berry.
The accuracy of the Gaussian approximation to the sum of
independent variates.
Transactions of the American Mathematical Society, 490
(1):0 122–136, 1941.

[Beygelzimer et al., 2011]
Alina Beygelzimer, John Langford, Lihong Li, Lev Reyzin, and Robert E.
Contextual bandit algorithms with supervised learning guarantees.
In International Conference on Artificial Intelligence and
Statistics (AISTATS), volume 15, pages 19–26, 2011.

[Dobrushin, 1956]
Roland L Dobrushin.
Central limit theorem for nonstationary Markov chains.
Theory of Probability & Its Applications, 10
(1):0 65–80, 1956.

[Egorov et al., 2017]
Maxim Egorov, Zachary N. Sunberg, Edward Balaban, Tim Allan Wheeler, Jayesh K.
Gupta, and Mykel J. Kochenderfer.
POMDPs.jl: A framework for sequential decision making under
Journal of Machine Learning Research, 18:0 26:1–26:5,

[Esseen, 1942]
Carl-Gustav Esseen.
On the Liapunoff limit of error in the theory of probability.
Arkiv for matematik, astronomi och fysik, A: 1–19, 1942.

[Feldman and Domshlak, 2014]
Zohar Feldman and Carmel Domshlak.
Simple regret optimization in online planning for Markov decision
Journal of Artificial Intelligence Research, 51:0
165–205, 2014.

[Gosset, 1908]
William Sealy (Student) Gosset.
The probable error of a mean.
Biometrika, 60 (1):0 1–25, 1908.
ISSN 00063444.

[Grill et al., 2016]
Jean-Bastien Grill, Michal Valko, and Rémi Munos.
Blazing the trails before beating the path: Sample-efficient
Monte-Carlo planning.
In Advances in Neural Information Processing Systems
(NeurIPS), pages 4673–4681, 2016.

[Grill et al., 2019]
Jean-Bastien Grill, Omar Darwiche Domingues, Pierre Ménard, Rémi
Munos, and Michal Valko.
Planning in entropy-regularized Markov decision processes and
In Advances in Neural Information Processing Systems
(NeurIPS), pages 12383–12392, 2019.

[Hölder, 1889]
Otto Hölder.
Ueber einen Mittelwertsatz.
Gött. Nachr., 1889:0 38–47, 1889.

[Kaufmann and Koolen, 2017]
Emilie Kaufmann and Wouter M. Koolen.
Monte-Carlo tree search by best arm identification.
In Advances in Neural Information Processing Systems
(NeurIPS), pages 4897–4906, 2017.

[Kearns et al., 1999]
Michael J. Kearns, Yishay Mansour, and Andrew Y. Ng.
A sparse sampling algorithm for near-optimal planning in large
Msarkov decision processes.
In International Joint Conference on Artificial Intelligence
(IJCAI), pages 1324–1231, 1999.

[Kocsis and Szepesvári, 2006]
Levente Kocsis and Csaba Szepesvári.
Bandit based Monte-Carlo planning.
In European Conference on Machine Learning (ECML), 2006.

[Lai and Robbins, 1985]
Tze Leung Lai and Herbert Robbins.
Asymptotically efficient adaptive allocation rules.
Advances in Applied Mathematics, 60 (1):0
4–22, 1985.

[Markov, 1951]
Andrey A. Markov.
Investigations of general experiments connected by a Markov chain.
Collected Works, pages 465–509, 1951.

[Maurer and Pontil, 2009]
Andreas Maurer and Massimiliano Pontil.
Empirical Bernstein bounds and sample-variance penalization.
In Conference on Learning Theory (COLT), 2009.

[Nelder and Mead, 1965]
John A. Nelder and R. Mead.
A simplex method for function minimization.
Computer Journal, 70 (4):0 308–313, 1965.

[Peel et al., 2013]
Thomas Peel, Sandrine Anthoine, and Liva Ralaivola.
Empirical Bernstein inequality for martingales : Application to
online learning.
HAL, 2013.

[Shah et al., 2020]
Devavrat Shah, Qiaomin Xie, and Zhi Xu.
Non-asymptotic analysis of Monte Carlo tree search.
In Joint International Conference on Measurement and Modeling
of Computer Systems (SIGMETRICS), pages 31–32, 2020.

[Shevtsova, 2011]
I. Shevtsova.
On the absolute constants in the Berry-Esseen-type inequalities.
Doklady Mathematics, 89:0 378–381, 2011.

[Silver and Veness, 2010]
David Silver and Joel Veness.
Monte-Carlo planning in large POMDPs.
In Advances in Neural Information Processing Systems (NIPS),

[Silver et al., 2016]
David Silver, Aja Huang, Chris J. Maddison, Arthur Guez, Laurent Sifre, George
van den Driessche, Julian Schrittwieser, Ioannis Antonoglou, Vedavyas
Panneershelvam, Marc Lanctot, Sander Dieleman, Dominik Grewe, John Nham, Nal
Kalchbrenner, Ilya Sutskever, Timothy P. Lillicrap, Madeleine Leach, Koray
Kavukcuoglu, Thore Graepel, and Demis Hassabis.
Mastering the game of Go with deep neural networks and tree search.
Nature, 5290 (7587):0 484–489, 2016.

[Silver et al., 2017]
David Silver, Julian Schrittwieser, Karen Simonyan, Ioannis Antonoglou, Aja
Huang, Arthur Guez, Thomas Hubert, Lucas Baker, Matthew Lai, Adrian Bolton,
Yutian Chen, Timothy P. Lillicrap, Fan Hui, Laurent Sifre, George van den
Driessche, Thore Graepel, and Demis Hassabis.
Mastering the game of go without human knowledge.
Nature, 5500 (7676):0 354–359, 2017.

[Sinn and Chen, 2013]
Mathieu Sinn and Bei Chen.
Central limit theorems for conditional Markov chains.
In International Conference on Artificial Intelligence and
Statistics (AISTATS), pages 554–562, 2013.

[Sutton and Barto, 2018]
Richard S. Sutton and Andrew G. Barto.
Reinforcement Learning: An Introduction.
MIT Press, 2 edition, 2018.

[Tichavsky et al., 1998]
P. Tichavsky, C.H. Muravchik, and A. Nehorai.
Posterior Cramér-Rao bounds for discrete-time nonlinear
IEEE Transactions on Signal Processing, 460
(5):0 1386–1396, 1998.

[Welch, 1947]
Bernard Lewis Welch.
The Generalization of Student's problem when several different
population variances are involved.
Biometrika, 340 (1-2):0 28–35, 1947.
ISSN 0006-3444.

[Zhang et al., 2021]
Zihan Zhang, Jiaqi Yang, Xiangyang Ji, and Simon S. Du.
Variance-aware confidence set: Variance-dependent bound for linear
bandits and horizon-free bound for linear mixture MDP.
Computing Research Repository, 2021.

§ PROOFS

This section provides proofs the presented action error bounds and shows how the same procedure can be applied to prove the value error bound.
We first present several lemmas.
The first lemma presents a bound on the difference in expected error of action value estimates generated by a Monte Carlo search.
lemma]lem: exp error
Given a Monte Carlo estimate $\bar{Q}_i$ of true value $Q_i$, the expected error is

\begin{equation}
\mathbb{E}[\bar{Q}_i - Q_i] \leq \frac{1}{|\mathcal{C}_i|}\sum_k \gamma^{H_k}\epsilon^0_k
% \mathbb{E}[\bar{Q}_i - Q_i] - \mathbb{E}[\bar{Q}_j - Q_j] \leq \frac{1}{|\mathcal{C}_i|}\sum_k \gamma^{H_k}\epsilon^0_k + \frac{1}{|\mathcal{C}_j|} \big|\sum_l\gamma^{H_l}\epsilon^0_l\big|
\end{equation}

where $\epsilon_k^0$ is the expected error of the baseline estimate at the end of the trajectory $k$ of depth $H_k$.

A proof of <ref> is given in <ref>.
The next lemma shows that an appropriate upper bound on the conditional likelihood of an estimator of an unknown parameter is also an upper bound on a probability of the parameter.
lemma]lem: cramer
Let $\theta$ be a population parameter and let $\bar{\theta}$ be a potentially biased estimator of that parameter with variance $\bar{\sigma}^2$.
If the probability

\begin{align}
P(\bar{\theta} \geq a \mid  \theta = b) \leq c(\bar{\sigma}^2)
\end{align}

holds, where $c(\bar{\sigma}^2)$ is a function that increases monotonically with $\bar{\sigma}^2$, then the inequality

\begin{align}
P(\theta \geq b \mid \bar{\theta} = a) \leq c(\bar{\sigma}^2)
\end{align}

also holds.

<Ref> follows from the Bayesian Cramér-Rao inequality [Tichavsky et al., 1998, Aras et al., 2019], which, for appropriate priors, implies

\begin{equation}
\mathrm{Var}[\theta \mid X] \leq \mathrm{Var}[\bar{\theta} \mid X]
\end{equation}

for distributional parameter $\theta$, estimator $\bar{\theta}$, and data $X$.
The next lemma bounds the true variance of an expectation estimator around the sample variance.
lemma]lem: variance
Let $\bar{X} = \frac{1}{n}\sum_i^n X_i$ be an estimator of the mean over $n$ random variables with variance $\bar{\sigma}^2$.
Given a estimate $\bar{x} = \frac{1}{n}\sum_i^n x_i \sim X_i$ and sample variance $V_n = \frac{1}{n}\sum_i^n (x_i - \bar{x})^2$, the variance is bounded as

\begin{equation}
\bar{\sigma}^2 \leq \frac{b^2}{n}\Bigg(\frac{1}{b}V_n^\frac{1}{2} + \sqrt{\frac{-\ln{\alpha}}{n - 1}}  \Bigg)^2
\end{equation}

where $b = \max_i x_i - \min_i x_i$ is the range of the sample values and $\alpha \in (0, 1)$ is a significance level.

<Ref> follows from the sample variance penalty for empirical Bernstein bounds stated in theorem 10 of the work by Maurer and Pontil, 2009.
The final lemma provides an upper bound to the finite sample penalty of the Berry-Esseen theorem when the true second and third moments are not known.
lemma]lem: berry-esseen
Given a set of zero-mean $i.i.d.$ random variables $X_1, \dots, X_n$ with variance $\sigma^2$ and $Y = 1/n sum_i X_i$, the difference in between the true cumulative distribution function $F_n$ and the normal Gaussian CDF $\Phi$ is bounded as

\begin{equation}
\sup_{x \in \mathbb{R}}|F_n(x) - \Phi(x)| \leq \frac{C}{n}
\end{equation}

where for all $n > 0$ and some value of $C \leq 0.4748$.

§.§ Proof of <Ref>

The proof begins by constructing defining a bound as

\begin{align}
& P\big(\bar{Q}_i - Q_i \geq \epsilon \big)~\label{proof start 1} \\
= & P\big(\bar{Q}_i - Q_i - \mathbb{E}\big[\bar{Q}_i \big] \geq \epsilon - \mathbb{E}\big[\bar{Q}_i \big] \big) \\
= & P\big(\bar{Q}_i - \mathbb{E}\big[\bar{Q}_i \big] \geq \epsilon - \mathbb{E}\big[\bar{Q}_i - Q_i\big] \big) \\
\leq & P\big(\bar{Q}_i - \mathbb{E}\big[\bar{Q}_i \big] \geq \epsilon - \sum_{k \in Ch_i}\gamma^{H_k}\epsilon^0_k \big)~\label{error step 1} \\
\leq & P\big(\bar{Q}_i - \mathbb{E}\big[\bar{Q}_i \big] \geq \epsilon - \sum_{k \in Ch_i}\gamma^{H_k}\epsilon^U_k \big)~\label{error bound step 1}
\end{align}

where <ref> follows from <ref>.
The expected baseline estimator error $\epsilon^0$ is not generally known, though its magnitude can be limited by a term $\epsilon^U \geq |\epsilon^0|$.
Applying this inequality results in <ref>.

We can bound the probability of the zero-mean quantity in <ref> with a Chernoff-Hoeffding bound as

\begin{align}
P\big(\bar{Q}_i - Q_i \geq \epsilon \big)
\leq & \exp{\Big(\frac{-2n_i (|\epsilon - \zeta|^+)^2}{n_i\bar{\sigma}_i^2}\Big)}~\label{chernoff step 1}
\end{align}

where $\bar{\sigma}^2$ gives the variance of the value estimators.
Since the true variance is not generally known, we apply <ref> with a union bound to arrive at

\begin{align}
\leq & \exp{\Big(\frac{-2n_i (|\epsilon - \zeta|^+)^2}{n_i\hat{\sigma}_i^2}\Big)} + \alpha~\label{variance step}
\end{align}

which concludes the proof.

§.§ Proof of <Ref>

The proof begins by constructing a bound on the likelihood of search results as

\begin{align}
& P\big(\bar{Q}_i - \bar{Q}_j \geq \delta \mid Q^\pi_j - Q^\pi_i = \epsilon \big)~\label{proof start} \\
= & P\big(\bar{Q}_i - \bar{Q}_j - \mathbb{E}\big[\bar{Q}_i - \bar{Q}_j \big] \geq \delta - \mathbb{E}\big[\bar{Q}_i - \bar{Q}_j \big] \big) \\
= & P\big(\bar{Q}_i - \bar{Q}_j - \mathbb{E}\big[\bar{Q}_i - \bar{Q}_j \big] \geq \delta + \epsilon - \mathbb{E}\big[\bar{Q}_i - Q_i\big] + \mathbb{E}\big[\bar{Q}_j - Q_j\big] \big) \\
\leq & P\big(\bar{Q}_i - \bar{Q}_j - \mathbb{E}\big[\bar{Q}_i - \bar{Q}_j \big] \geq \delta + \epsilon - \sum_{k \in Ch_i}\gamma^{H_k}\epsilon^0_k + \sum_{l \in Ch_j}\gamma^{H_l}\epsilon^0_l \big)~\label{error step} \\
\leq & P\big(\bar{Q}_i - \bar{Q}_j - \mathbb{E}\big[\bar{Q}_i - \bar{Q}_j \big] \geq \delta + \epsilon - \sum_{k \in Ch_i}\gamma^{H_k}\epsilon^U_k - \sum_{l \in Ch_j}\gamma^{H_l}\epsilon^U_l \big)~\label{error bound step}
\end{align}

where the conditional dependence term is omitted after <ref> for clarity.
The  <ref> follows from <ref> and the convergence condition.
The expected baseline estimator error $\epsilon^0$ is not generally known, though its magnitude can be limited by a term $\epsilon^U \geq |\epsilon^0|$.
Applying this inequality results in <ref>.

We use this term with <ref> to construct a Chernoff-Hoeffding bound on the sub-optimality as

\begin{align}
P\big(Q_j - Q_i \geq \epsilon \mid \bar{Q}_i - \bar{Q}_j = \delta \big)
\leq & \exp{\Big(\frac{-2n_i n_j(|\delta + \epsilon - \zeta|^+)^2}{n_i\bar{\sigma}_i^2 + n_j\bar{\sigma}_j^2}\Big)}~\label{chernoff step}
\end{align}

where $\bar{\sigma}^2$ gives the variance of the value estimators.
Since the true variance is not generally known, we apply <ref> with a union bound to arrive at

\begin{align}
\leq & \exp{\Big(\frac{-2n_i n_j(|\delta + \epsilon - \zeta|^+)^2}{n_i\hat{\sigma}^{2, \alpha}_i + n_j\hat{\sigma}^{2, \alpha}_j}\Big)} + \alpha~\label{variance step 1}
\end{align}

which concludes the proof.

<Ref> can be proved for action $a_i$ by applying the same procedure while omitting the terms related to $a_j$.

§.§ Proof of <Ref>

We begin this proof at <ref> of the general action error proof.

\begin{align}
& P\big(\bar{Q}_i - \mathbb{E}\big[\bar{Q}_i \big] \geq \epsilon - \sum_{k \in Ch_i}\gamma^{H_k}\epsilon^U_k \big) \\
\leq & 1 - \Phi(\epsilon - \zeta_i, 0, \frac{\bar{\sigma}^2}{n_i}) + \frac{C}{n_i}~\label{eq: CLT} \\
\leq & 1 - \Phi(\epsilon - \zeta_i, 0, \frac{\hat{\sigma}^{2, \alpha}}{n_i}) + \alpha + \frac{C}{n_i}~\label{eq: CLT union}
\end{align}

The inequality of <ref> follows from the definition of the cumulative distribution function and application of <ref>.
The inequality of <ref> follows from application of the variance upper bound of <ref> with the union bound, concluding the proof.

§.§ Proof of <Ref>

We begin this proof at <ref> of the general action error proof.

\begin{align}
& P\big(\bar{Q}_i - \bar{Q}_j - \mathbb{E}\big[\bar{Q}_i - \bar{Q}_j \big] \geq \delta + \epsilon - \sum_{k \in Ch_i}\gamma^{H_k}\epsilon^U_k - \sum_{l \in Ch_j}\gamma^{H_l}\epsilon^U_l \big) \\
\leq & 1 - \Phi(\delta + \epsilon - \zeta_i - \zeta_j, 0, \frac{\bar{\sigma}^2}{n_i + n_j}) + \frac{C}{n_i + n_j}~\label{eq: CLT 1} \\
\leq & 1 - \Phi(\delta + \epsilon - \zeta_i - \zeta_j, 0, \frac{\hat{\sigma}^{2, \alpha}}{n_i + n_j}) + \alpha + \frac{C}{n_i + n_j}~\label{eq: CLT union 1}
\end{align}

The inequality of <ref> follows from the definition of the cumulative distribution function and application of <ref>.
The inequality of <ref> follows from application of the variance upper bound of <ref> with the union bound, concluding the proof.

§.§ Proof of <Ref>

Define $\phi$ to be a policy that follows the Monte Carlo sampling strategy until a depth $h$ and then follows the optimal policy.
Using this policy, we can define the expected error of an action value estimate to be

\begin{align}
\mathbb{E}[\bar{Q}_i - Q_i] &= \mathbb{E}\Big[\frac{1}{n}\sum_{j \in Ch} \bar{v}_j - Q_i\Big] \\
&= \mathbb{E}\Big[\frac{1}{n}\sum_{j \in Ch} \bar{v}_j - Q^\phi_i\Big] + Q^\phi_i - Q_i \\ % + \mathbb{E}[Q^\phi_i - Q_i] \\
&= \frac{1}{n}\sum_{j \in Ch}\gamma^{H_j} \mathbb{E}\Big[(\hat{v}_j - v_j)\Big] + (Q^\phi_i - Q_i) \\% + \mathbb{E}[Q^\phi_i - Q_i] \\
&= \frac{1}{n}\sum_{j \in Ch} \gamma^{H_j} \epsilon^0_j + (Q^\phi_i - Q_i)\\
&\leq \frac{1}{n}\sum_{j \in Ch} \gamma^{H_j} \epsilon^0_j~\label{err def}
\end{align}

where $v_j$ and $\hat{v}_j$ are the true and estimated values at the end of the sample trajectory, respectively.
Intuitively, this is a decomposition of the expected error to the error of the baseline estimator and the expected error caused by the Monte Carlo sampling procedure.
The inequality in <ref> follows from the fact the the value the optimal policy value $Q_i \geq Q^\phi_i$ for any policy $\phi$.
Under the Monte Carlo convergence conditions, $\mathbb{E}[Q^\phi_i] = Q^\pi_i$.

§.§ Proof of <Ref>

We begin this proof with the statement of the Berry-Esseen theorem with known distributional parameters, with the third moment denoted $\rho = \mathbb{E}[|X|^3]$.

\begin{align}
\sup_{x \in \mathbb{R}}|F_n(x) - \Phi(x)| & \leq \frac{C\rho}{\sqrt{n}\sigma^3} \\
& = \Big( \frac{C}{\sqrt{n}} \Big) \Big( \frac{(1/n) \|X\|_3^3}{(1/n)^{3/2} \|X\|_2^3}\Big) \\
& \leq \Big( \frac{C}{n} \Big)~\label{eq: norms}
\end{align}

The inequality in <ref> follows from the fact that $\|X\|_p \geq \|X\|_q$ for any $p \leq q$, which results from Hölder's inequality [Hölder, 1889].


1Indian Institute of Astrophysics, II Block, Koramangala, Bangalore 560034,
India.
3Indian Centre for Space Physics, 43 Chalantika, Garia Stn. Road, Kolkata
700084, India.

# Study of accretion flows around an ultraluminous X-ray source M82 X-1 using
NuSTAR data

Santanu Mondal1,* Biswaraj Palit1 and Sandip K. Chakrabarti2

###### Abstract

We study the spectral properties and accretion flow behavior of an
ultraluminous X-ray source M82 X-1 using NuSTAR observations. We use the
physical two component advective flow (TCAF) model to fit the data and to
derive the accretion flow properties of the source. From the model fitted
parameters, we found that M82 X-1 is harboring an intermediate mass black hole
at its centre, where the mass varies from $156.04^{+13.51}_{-15.30}$ to
$380.96^{28.38}_{-29.76}$ M⊙. The error weighted average mass of the black
hole is $273\pm 43$ M⊙, which accreted in nearly super-Eddington rate. The
Compton cloud was compact with a size of $\sim 13r_{g}$ and the shock
compression ratio had nearly intermediate values except for the epoch four.
These indicate a possible significant mass outflow from the inner region of
the disk. The quasi periodic oscillation (QPO) frequencies estimated from the
model fitted parameters can reproduce the observed QPOs. The robustness of the
model parameters is verified by drawing the confidence contours among them.

###### keywords:

accretion, accretion discs — black hole physics — stars: black holes —
X-rays:individual (M82 X-1)

<EMAIL_ADDRESS><EMAIL_ADDRESS>
#### 000 0000 1– 1

##  1 Introduction

Ultraluminous X-ray sources (ULXs) appear to be very luminous and non-nuclear
extra-galactic objects. Most of the ULXs are thought to be X-ray binaries
powered by accretion onto a central compact object which can be a stellar-mass
black hole (StMBH) or a massive black hole (MBH), and their luminosities are
comparable with or above the Eddington luminosity of stellar black holes. It
was also proposed that ULXs could be powered by accreting intermediate-mass
black holes (IMBHs: $\sim 100-10^{6}$ M⊙; Mezcua (2017) for a review) in the
sub-Eddington regime (Colbert & Mushotzky, 1999; Matsumoto & Tsuru, 2001). The
plausible ULX IMBH candidates might be those at the high-end of the high mass
X-ray binary (HMXB) luminosity distribution (see Mineo $et~{}al$., 2012), with
luminosity $\gtrsim$ $10^{41}$ erg s-1 in some cases. Some of these candidates
seem to show variability properties consistent with the Galactic black hole
binaries (GBHBs; Godet $et~{}al$., 2009; Pasham $et~{}al$., 2014) and the
presence of cold accretion disks (Godet $et~{}al$., 2012), suggesting masses
in the IMBH regime. However, it is not the case always, and it has been
evidenced that a large population of ULXs does not comply with canonical state
transition like GBHBs (Stobbart $et~{}al$., 2006; Gladstone $et~{}al$., 2009,
and references therein). Therefore, another possibility came up, which
suggests that accretion onto StMBH in the super-Eddington regime could power
these sources (Shakura & Sunyaev, 1973; King & Pounds, 2003; Poutanen
$et~{}al$., 2007). Apart from extra-galactic evidence, there is now evidence
for a Galactic ULX (Wilson-Hodge $et~{}al$., 2018). However, most of their
features remain unknown.

M82 X-1 is an ultraluminous X-ray source present in the cigar-shaped galaxy
M82, which is a star burst galaxy with a high star formation rate (SFR) $\sim$
10M⊙ yr-1 (Iwasawa, 2021). It is located at a distance of 3.6 Mpc (Freedman
$et~{}al$., 1994). Matsumoto & Tsuru (2001) reported the presence of a low
luminosity Active Galactic Nuclei (LLAGN) in M82-galaxy. In the proximity of
M82 X-1 source lies several other X-ray luminous sources (namely X-2, X-3 and
X-4 (see Fig. 1; Brightman $et~{}al$., 2020)). M82 X-1 can achieve an X-ray
luminosity as high as $\sim$ 1041 erg s-1, making it the brightest among all
other sources. Most of the research on ULXs is focused on identifying the
properties of the central compact object and the nature of accretion flow
which is capable of producing such high X-ray luminosities. Similarly, several
intriguing features of M82 X-1 have been employed to determine its mass and
accretion geometry. Analysis of Rossi X-Ray Timing Explorer (RXTE)
observations by Kaaret $et~{}al$. (2006b) had shown a 62.5 day X-ray orbital
period, while Qiu $et~{}al$. (2015) and more recently, Brightman $et~{}al$.
(2019) claimed that there is degeneracy regarding the source of this orbital
flux modulation and instead the period may be super-orbital in nature
originating from M82 X-2 which is now a confirmed pulsar (Bachetti $et~{}al$.,
2021). M82 X-2, which is the second most luminous source in this galaxy, is
only separated from M82 X-1 by a mere $5^{\prime\prime}$ angular separation.
These two factors contribute to an enormous amount of contamination in
observed properties of the source. Currently, only Chandra satellite has the
resolving power required to distinguish X-1 from X-2. In order to constrain
the amount of contamination, spectral studies were required to be conducted
above 10 keV.

Further, it was claimed that, M82 X-1 might be located in a young star
cluster; MGG-11 (McCrady $et~{}al$., 2003). Patruno $et~{}al$. (2006)
simulated the binary evolution of this star cluster, with compact object mass
range between 10-5000 M⊙ and donor star mass range 18-26 M⊙. The results from
their simulations agreed with observations from MGG-11 and required that this
dense stellar cluster hosted highly luminous BH binary systems with donor mass
$\sim$ 22 M⊙ and the BH mass range of 200-5000 M⊙. In addition, they also
argued that even with a lower IMBH range (200-400 M⊙), it requires either
beaming of disk outflow or mild super-Eddington accretion to explain its high
luminosity. Hence, if M82 X-1 belonged to this young dense cluster, then
evidences favour an IMBH scenario.

The bolometric luminosity of M82 X-1 exceeds the Eddington luminosity for a 10
M⊙ BH by at least two orders of magnitude (Dewangan $et~{}al$., 2006). If
indeed M82 X-1 is an IMBH accreting close to super-Eddington rate, then the
mass is $\geq$ 700 M⊙ (Okajima $et~{}al$., 2006). Considering a rapidly
rotating BH at moderate super-Eddington accretion rates, Feng & Kaaret (2010)
derived a mass range of $200-800$ M⊙. Their fitting (kerrbb model; Li
$et~{}al$., 2005) results agreed with L $\propto$ T4 variation which indicated
a blackbody emission from the disk, contrasting the results of Matsumoto
$et~{}al$. (2003) who had reported no correlation among these quantities.
However, Okajima $et~{}al$. (2006) pointed out that within the slim-disk
formulation (used to explain the high disk temperature), even a massive
stellar mass BH (MStMBH) can explain spectral curvature beyond 3 keV. Using
high quality XMM-Newton observations, they reported a p-value of 0.61 (where p
is the exponent in the relation T $\propto$ r-p). Radiation hydrodynamic
simulations supported the formation of disk accreting in super-Eddington rate
(Ohsuga $et~{}al$., 2005).

A joint analysis of XMM-Newton and INTEGRAL satellite, conducted by Sazonov
$et~{}al$. (2014), indicated a steep energy rollover at $\sim$ 10 keV. This is
very typical of ULXs. They also suggested the existence of cold gas along the
line of site of M82 X-1. Recently, Brightman $et~{}al$. (2020) carried out an
extensive spectral analysis of source M82 X-1 using simultaneous Chandra,
NuSTAR, and Swift/XRT data. Results from their spectral fitting showed that
DISKPBB*SIMPL model gave the best fits to spectra of M82 X-1 when the
contamination from M82 X-2 was less (when it was assumed to be in an off
state). In an earlier work, Brightman $et~{}al$. (2016) calculated the mass of
BH in M82 X-1 from normalization-Radius relation of DISKPBB model fit and
concluded that super-Eddington accreting StMBH and estimated the mass as 26 M⊙
and 125 M⊙ respectively. All estimates pointed towards a MStMBH. Overall, it
was established that a puffy disk most accurately explained the spectral
features of M82 X-1, which in turn led to the conclusion that this ULX is
undergoing super-Eddington accretion. Thus, two very important features of ULX
have been identified, one being a high energy cutoff and other super-Eddington
accretion.

Along with the spectral variability, M82 X-1 also showed variability in its
lightcurve. Quasi-periodic oscillations (QPOs) in M82 X-1 were first reported
by Strohmayer & Mushotzky (2003) in the range of very-low frequency (LF)
$\sim$ 54 mHz. Another QPO at $\sim$ 100 mHz was found by Fiorito & Titarchuk
(2004) associated with high/soft states ($\Gamma$ $\sim$ 2.1-2.7). Using the
LFQPO frequency-Spectral Index correlation technique, they suggested that M82
X-1 harbours intermediate mass black holes (IMBHs) $\sim$ 1000 M⊙. Mucciarelli
$et~{}al$. (2006) showed that mHz range QPOs varied overtime, ranging from 50
mHz to 170 mHz. This corresponded to mass ranging from 10 M⊙ to 1000 M⊙. The
mass estimation of 25-520 M⊙ done by Dewangan $et~{}al$. (2006) using the
photon index v/s QPO frequency relation falls in the above range. More
recently, “twin-peaked” stable 3:2 ratio QPOs; 3.32 Hz and 5.07 Hz were
detected for this source by Pasham $et~{}al$. (2014); yielded mass roughly
$\sim$ 400 M⊙.

Since the proposition of mass outflow from inner regions of an advected disk
in ULXs by Poutanen $et~{}al$. (2007), several ULXs (NGC 1313 X-1 Pinto
$et~{}al$. (2016), NGC 5408 X-1 Pinto $et~{}al$. (2016), and NGC 6946 X-1
Kosec $et~{}al$. (2018)) have been found to exhibit soft-excess below 2 keV,
which indicates a wind dominated accretion disk. However, since the soft X-ray
emission from M82 X-1 coincides with the diffuse plasma emission (Miller,
2005), spectral information below 2 keV cannot be trusted. Hence, it has been
difficult to study M82 X-1 in the context of wind launching from accretion
disks. Only recently Brightman $et~{}al$. (2020) has hinted at a possible
positive correlation between mass outflow rate and absorption by galactic
column density. So far, we noticed that investigating temporal properties have
mostly favoured IMBH while simple multicoloured DISK+POWERLAW continuum models
or hybrid-disk models DISKPBB have favoured MStMBHs. The accretion onto these
objects is exhaustively studied, and several phenomenological models (see
Kaaret $et~{}al$., 2017, for a review) are proposed in the literature,
however, we know very little about the disk formation, origin of inner hot
region, its optical depth, geometry, etc. Therefore, we require a self-
consistent models that are able to explain the harder energy range spectra,
which are less contaminated.

Despite the fact that the source has been studied in different energy bands,
we know very little about the underlying physical processes that created such
variable spectral features. Also, the literature estimated a broad range of
mass of the compact object, leaving a range of questions: (1) does this system
harbour a StMBH or an IMBH? (2) what is the nature of accretion behind such a
high luminosity? (3) is the high column density along the line of sight (LOS)
due to the wind from the disk, or due to some moving cloud along the LOS? In
this paper, we aim to address different observational aspects after performing
spectral fitting of NuSTAR  data using two component advective flow (TCAF)
model (Chakrabarti & Titarchuk, 1995). The paper is organized as follows: In
Section 2, we discuss the observation log, data analysis procedures, and
modelling. In Section 3, we present our spectral fitting results and different
estimations. Finally, in Section 4, we make our concluding remarks.

##  2 Observation and data analysis

For the spectral analysis of M82 X-1 (located at RA=148.958 and Dec=+69.679 in
the galaxy M82) of each epoch of observation, we used the data of NuSTAR
(Harrison $et~{}al$., 2013) during 2016 in the energy range of 3$-$30 keV.
NuSTAR data was extracted using the standard NuSTAR Data Analysis Software
(NuSTARDAS v1.3.1) 111https://heasarc.gsfc.nasa.gov/docs/nustar/analysis/
software. The exposure time and the observation log are given in Table 1. We
ran nupipeline task to produce cleaned event lists and nuproducts task to
generate the source and background spectra with the standard filters. Source
events were selected from a circular region of 50′′. The background region
were selected from larger circular source-free region (of 60′′) to avoid
contamination. The data was grouped by grppha command, with a minimum of 25
counts in each bin. To avoid contamination of the source from X-2, we
considered those observations which have much higher flux compared to X-2,
following the data selection discussed in Brightman $et~{}al$. (2020).

We used XSPEC222https://heasarc.gsfc.nasa.gov/xanadu/xspec/ (Arnaud, 1996)
version 12.11.0 for spectral analysis. Each epoch of observation was fitted
using both phenomenological and physical models, which is discussed in the
later sections. For all observations, we used the absorption model tbabs
(Wilms $et~{}al$., 2000), for the Galactic absorption, keeping the hydrogen
column density ($N_{H}$) fixed to $0.9\times 10^{22}$cm-2 (Brightman
$et~{}al$., 2020; Kalberla $et~{}al$., 2005, the Leiden/Argentine/Bonn (LAB)
survey) during the fitting. In addition to the above $N_{H}$, the host galaxy
$N_{H}$ is also taken into account using the partial covering ionization
(ZXIPCF) model (Reeves $et~{}al$., 2008), which is discussed later. We used
chi-square statistics for the goodness of the fitting.

Table 1: NuSTAR observations used in this paper Epoch | OBSID | Date | MJD | Exposure
---|---|---|---|---
E1 | 80202020002 | 2016-01-26 | 57413 | 36 ksec
E2 | 30202022002 | 2016-06-03 | 57542 | 39 ksec
E3 | 30202022008 | 2016-07-29 | 57598 | 42 ksec
E4 | 90202038002 | 2016-10-07 | 57688 | 45 ksec

###  2.1 Modelling with TCAF

There are several models available in the literature which have been used to
study the source M82 X-1. Alternatively, a more physical model, namely a TCAF
model is also present in the literature and this provides a natural
explanation of two phase configuration. For the spectral fitting, we used TCAF
model along with a multiplicative ZXIPCF model.

According to TCAF, matter falling onto a BH has two components, one is the
high angular momentum Keplerian flow sitting at the equatorial plane and the
second component is the hot, low angular momentum, sub-Keplerian flow, which
forms the inner hot region after forming a shock due after satisfying Rankine-
Hugoniot shock conditions (see Chakrabarti, 1989). The boundary layer of this
region has physical origin rather than phenomenological, which is called
CENtrifugal pressure supported BOundary Layer (CENBOL). This model not only
considers accretion, but also launches jet in a self-consistent way; therefore
it can explain accretion-ejection or disk-jet connection in BH systems (see
Chakrabarti, 1999; Mondal & Chakrabarti, 2021).

The model requires five parameters (if the mass is unknown and is considered
to be a free parameter) which are (i) mass of the BH ($M_{\rm BH}$), (ii) disk
accretion rate ($\dot{m}_{\rm d}$), (iii) halo accretion rate ($\dot{m}_{\rm
h}$), (iv) size of the CENBOL or the location of the shock ($X_{\rm s}$ in
$r_{g}=2GM_{\rm BH}/c^{2}$ unit), and (v) shock compression ratio ($R$). In
this model two-temperature equations for the electron and protons are solved
to compute the inverse-Comptonized spectra. The emergent spectrum from TCAF
has mainly three components, (1) multicolour blackbody spectrum which is
coming from the disc, (2) the hard radiation; from the upscattering of the
soft photons from the disk by the hot corona, and (3) the scattering of the
hard radiation by the cold disc; which is so-called reflection component. Till
date, TCAF model has been successfully applied to study the outburst
properties of StMBHs (Debnath $et~{}al$., 2014; Mondal $et~{}al$., 2014; Jana
$et~{}al$., 2016, and references therein), AGNs (Nandi $et~{}al$., 2019;
Mondal & Stalin, 2021), and estimated the mass of the central BHs (Molla
$et~{}al$., 2017; Nandi $et~{}al$., 2019). In the timing domain, the
oscillation of the shock in this model self-consistently explains the
occurrence of QPOs, and this interpretation in the context of thick disks was
proposed by Molteni $et~{}al$. (1996). Later, the shock oscillation model was
verified and applied extensively to several GBHBs (Chakrabarti $et~{}al$.,
2015; Mondal, 2020; Singh $et~{}al$., 2021, and references therein).
Considering the success of this model in explaining the observed spectral and
temporal properties of GBHBs, we now employ it to the ULX M82 X-1 as ULXs have
been compared with GBHBs. M82 X-1 also showed two Type-C QPO frequency, $\sim$
3.32 Hz and 5.07 Hz (Pasham $et~{}al$., 2014). All these observational
findings and evidence have prompted us to use this model and study the nature
of accretion flow of M82 X-1 which gives rise to such a high luminosity.

##  3 Results

###  3.1 Spectral fit

In Figure 1, we show the TCAF model fitted 3-30 keV spectra for all four
epochs. We have obtained acceptable fit statistics with reduced $\chi^{2}$
values ranging from $0.97$ to $1.15$, mentioned with the figures. The data
quality is noisy above 25 keV. The corresponding model fitted parameters are
provided in Table 2. It can be seen that the mass of the BH varied in a range
between 156 M⊙ and 381 M⊙, the error weighed value is $273\pm 43M_{\odot}$.
This implies that M82 is harboring an IMBH at its centre. The mass accretion
rates we obtained are in the super-Eddington regime with values of
$\dot{m}_{\rm d}$ from 2.09 to 2.53 $\dot{M}_{\rm Edd}$ and $\dot{m}_{\rm h}$
from 2.88 to 2.98 $\dot{M}_{\rm Edd}$. This agrees with the previous claims of
M82 X-1 being a super-Eddington accretor. The CENBOL was compact, varies in a
narrow range between $11$ and 15 $r_{g}$. The intermediate values of $R$,
indicate a possible mass outflow from the system (Chakrabarti, 1999). In
Figure 2, variation of all model parameters with MJD are shown. The top five
panels show the variation of the mass of the BH, disk and halo accretion
rates, location of the shock and shock compression ratio respectively, which
are variable during the observation period. The bottom panel shows covering
fraction from ZXIPCF model. Considering epochs E3 and E4, which are separated
by 68 days, we notice a significant change in $\dot{m}_{\rm d}$ from 2.42 to
2.09 $\dot{M}_{\rm Edd}$, and $R$ from 3.77 to 5.07. However, $X_{s}$ remained
more or less the same, which can be due to decrease in $\dot{m}_{\rm h}$. On
the contrary, in case of epochs E2 and E3, separated by 56 days, increase in
$\dot{m}_{\rm d}$ from 2.42 to 2.53 $\dot{M}_{\rm Edd}$ shifts the shock
location towards the BH, when the $\dot{m}_{\rm h}$ remained the same.

It should be noted that the ZXIPCF model fitted NH along the LOS and covering
factor ($C_{f}$) are significantly variable and high. The super-Eddington mass
accretion naturally explains that the wind was launched from the disk, which
blocked the central continuum radiation, therefore NH is also high (col. 8 in
Table 2). The anti-correlation between NH and $\dot{m}_{\rm d}$ can explain
the continuum luminosity driven outflow from the disk. The increase in
$\dot{m}_{\rm d}$, increases the luminosity of the continuum. Assuming that
the photon momentum is transferred to the wind from the disk, therefore it
increases the momentum of the wind outflow as well. A higher momentum flow can
escape the gravitational barrier of the disk and give a low value of NH along
the LOS, whereas, the effect is opposite when $\dot{m}_{\rm d}$ is low (see an
active galactic nucleus case study in Mondal $et~{}al$., 2022). However,
unlike many other ULXs, it has not been able to study M82 X-1 in the context
of wind launching from the inner disk from absorption line variability (Pinto
$et~{}al$., 2016), mainly due to diffuse emission of the background.

Figure 1: TCAF model fitted 3$-$30.0 keV spectra for all four epochs of M82
X-1 during NuSTAR era.  Table 2: Best fitted TCAF model parameters are
provided. Here, NH and Cf are hydrogen column density and dimensionless
covering fraction for the partial covering model respectively.

Obs Ids. | Epoch | $M_{\rm BH}$ | $\dot{m}_{\rm d}$ | $\dot{m}_{\rm h}$ | $X_{\rm s}$ | R | $N_{H}$ | Cf | $\chi^{2}/dof$
---|---|---|---|---|---|---|---|---|---
|  | $M_{\odot}$ | $\dot{M}_{\rm Edd}$ | $\dot{M}_{\rm Edd}$ | $r_{\rm g}$ |  |  $\times 10^{22}$ cm-2 |  |
80202020002 | E1 | $380.96^{+28.38}_{-29.76}$ | $2.494^{+0.142}_{-0.102}$ | $2.876^{+0.194}_{-0.192}$ | $14.37^{+1.35}_{-1.82}$ | $2.65^{+0.36}_{-0.30}$ | $21.29^{+1.16}_{-1.38}$ | $0.73^{+0.08}_{-0.05}$ | 317.6/310
30202022002 | E2 | $156.04^{+13.51}_{-15.30}$ | $2.417^{+0.03}_{-0.03}$ | $2.977^{+0.08}_{-0.07}$ | $14.87^{+1.61}_{-1.66}$ | $3.16^{+0.13}_{-0.21}$ | $23.77^{+1.13}_{-1.34}$ | $0.76^{+0.06}_{-0.10}$ | 334.4/299
30202022008 | E3 | $222.57^{+10.76}_{-13.14}$ | $2.527^{+0.08}_{-0.07}$ | $2.974^{+0.021}_{-0.106}$ | $11.06^{+1.33}_{-1.27}$ | $3.77^{+0.24}_{-0.12}$ | $18.23^{+0.91}_{-1.12}$ | $0.84^{+0.03}_{-0.04}$ | 395.1/342
90202038002 | E4 | $195.75^{+4.91}_{-4.75}$ | $2.090^{+0.057}_{-0.056}$ | $2.899^{+0.041}_{-0.042}$ | $11.42^{+1.48}_{-1.92}$ | $5.07^{+0.23}_{-0.24}$ | $28.47^{+0.61}_{-0.59}$ | $0.81^{+0.04}_{-0.04}$ | 292.5/301

Figure 2: Variation of the best fitted TCAF model parameters with MJD are
shown. The bottom panel shows the partial covering fraction parameter of
ZXIPCF model.

The BH mass proposed for M82 X-1 varies in a significantly broad range between
$\sim$10 to $10^{3}$ M⊙. However, TCAF model fitted BH mass varies in a narrow
range. To compare the TCAF model fitted BH mass with that present in the
literature, we plot all possible estimated masses in Figure 3. As discussed in
section 1, the powering of ULXs can be explained depending on the accretor
mass and its mass accretion rates at different regimes. In addition to other
possibilities, our spectral study shows that an IMBH can also accrete in the
super-Eddington limit. This estimate agrees with the claims in the literature
that super-Eddington accretion may not be limited to stellar mass BHs (Kaaret
$et~{}al$., 2006a). It is evident that super-massive BHs can accrete in super-
Eddington limit (Du $et~{}al$., 2015; Liu $et~{}al$., 2021, and references
therein) and the physical processes in low mass BHs to the supermassive BHs
can be scaled by their mass (McHardy $et~{}al$., 2006). Furthermore, in
principle, advective flows can accrete in the super-Eddington regime.
Observational evidence of the IMBH AGN source-RX J1140.1+0307 suggests it is
accreting at $\sim$ 10 times its Eddington limit (Jin $et~{}al$., 2016). One
of the often discussed models for growth of BHs into super-massive BHs (SMBHs)
is through rapid mass accretion in $\sim$ 100 M⊙ BHs (Greene $et~{}al$., 2020;
Toyouchi $et~{}al$., 2021). Such IMBHs may be formed after the death of the
earliest known Pop-III stars. Therefore, super-Eddington accretion onto IMBHs
is a necessary requirement to form SMBHs within appropriate timescales (Greene
$et~{}al$., 2020, and references therein). To further strengthen our claim we
have drawn confindence contour of the models parameters in the next section.

Figure 3: All possible estimates of the mass of M82 X-1 obtained from
different analyses in the literature. Red, green, and blue colored data points
represent mass estimates above 1000M⊙, between 1000 M⊙-100 M⊙ and less than
100 M⊙ respectively. Authors who have only reported a range of MBHs, for them
we have only plotted the lower limit and extended an arrow upwards towards the
upper limit where the length of the arrow is proportional to range of the mass
estimation. X-axis indicates the models used by the authors. The points inside
the rectangle represent the mass estimated in this paper.

###  3.2 Robustness of model fitted parameters

Figure 4 shows the confidence contours of $M_{BH}$ with other TCAF model
parameters. We have used steppar task in XSPEC to generate these plots. The
top panels show the confidence contours of $M_{BH}$ with disk and halo
accretion rates, and the bottom panels show the correlation with the location
of the shock and its compression ratio. Different model parameters are
labelled in x-axis. Three different contour colors (red, green, and blue)
correspond to one, two, and three sigma confidence levels. It is also the same
as $\Delta\chi^{2}$ fit statistics of 2.3, 4.61, and 9.21. The correlation of
$M_{BH}$ with other model parameters are clearly visible from the contour
plots. The model parameters assess better precision. The correlations are
clear, which show that there is no obvious degeneracy between parameters and
the solution is robust. This additional test reasonably agrees our new
possibility of powering the ULXs.

Figure 4: Confidence contours of $M_{\rm BH}$ with other model parameters for
the epoch E1 are shown. Three different contour colors show the confidence
level of one (red), two (green), and three sigma (blue). Here, x-labels follow
the parameter notations used in Table 2, e.g. mdotd, mdoth, and $Xs$
represents $\dot{m}_{d}$, $\dot{m}_{h}$, and $X_{s}$ respectively.

###  3.3 QPO frequency from model parameters

As explained by the propagatory shock oscillation model (Chakrabarti
$et~{}al$., 2005, 2008), the QPO frequencies follow the relation below:

$\nu_{\rm QPO}=\frac{c}{r_{\rm g}\;M_{\rm BH}\;RX_{\rm s}^{1.5}}$ (1)

Considering model fitted parameters from Table 2, we have estimated the QPO
frequencies for all four epochs, the $\nu_{\rm QPO}$ values range from 1.83 Hz
to 3.53 Hz (see, Table 3). The error weighed $\nu_{\rm QPO}$ is 2.93 $\pm$
0.32 Hz. In epochs E2 and E3, the values are closer to the lower twin-peaked
QPO (3.32 Hz) estimated by Pasham $et~{}al$. (2014) from power density
spectra. Therefore, our model fitted parameters are also able to reproduce the
observed QPOs. Conversely, it can be also be inferred that BH mass estimated
from the QPO frequency is along the line of our model fitted mass. As evident
from Figure 3, mass values of M82 X-1 dominate the 100M⊙ \- 400M⊙ region,
consistent with our estimate.

Table 3: Estimated QPO frequencies against those reported by Pasham $et~{}al$. (2014) Epoch | $\nu_{\rm QPO}$ | $\nu_{\rm QPO}$
---|---|---
| Estimated | Observed
| Hz | Hz
E1 |  1.83$\begin{subarray}{c}+0.38\\\ -0.42\end{subarray}$ |
E2 |  3.53$\begin{subarray}{c}+0.66\\\ -0.72\end{subarray}$ |  3.32 $\begin{subarray}{c}+0.06\\\ -0.06\end{subarray}$
E3 |  3.24 $\begin{subarray}{c}+0.63\\\ -0.59\end{subarray}$ |  5.07 $\begin{subarray}{c}+0.06\\\ -0.06\end{subarray}$
E4 |  2.61 $\begin{subarray}{c}+0.52\\\ -0.67\end{subarray}$ |

##  4 Conclusions

In this paper we study the accretion properties of the ultraluminous X-ray
source M82 X-1 using NuSTAR data with the TCAF model. Our main conclusions are
follows:

  * •

TCAF model can successfully fit and explain the ULXs data in the same way as
it did for Galactic black hole binaries and active galactic nuclei.

  * •

The model fitted mass falls in the intermediate BH mass range between 156 M⊙
and 381 M⊙.

  * •

Our estimates suggest that both the Keplerian and sub-Keplerian mass accretion
rates are nearly in the super-Eddington range.

  * •

The CENBOL was found to be compact around 13 rg during the observation period.

  * •

The shock compression ratio ($R$) falls in the intermediate range indicates
that the mass outflow from the CENBOL region was significant (Chakrabarti,
1999).

  * •

The confidence contours of mass with other parameters points to the robustness
of the estimation of the mass and the other model parameters.

  * •

The super-Eddington mass accretion can be the reason of disk wind as observed
in the literature.

  * •

The QPO frequencies estimated from the TCAF model fitted parameters agree well
with the observations in the literature.

This study simultaneously explains both spectral and temporal properties of
the ULX M82 X-1 using the same model parameters which have been obtained from
spectral fitting.

## Acknowledgements

We thank the referee for making constructive comments and suggestions. SM
thanks M. Brightman for discussions on resolving the source. SM and BP
acknowledge Ramanujan Fellowship research grant (File #RJF/2020/000113) by
DST-SERB, Govt. of India for this research. This research has made use of the
NuSTAR Data Analysis Software (nustardas) jointly developed by the ASI Science
Data Center (ASDC), Italy and the California Institute of Technology
(Caltech), USA. This research has also made use of data obtained through the
High Energy Astrophysics Science Archive Research Center Online Service,
provided by NASA/Goddard Space Flight Center.

## References

  * Arnaud (1996) Arnaud, K. A. 1996, Astronomical Society of the Pacific Conference Series, Vol. 101, XSPEC: The First Ten Years, ed. G. H. Jacoby & J. Barnes, 17
  * Bachetti $et~{}al$. (2021) Bachetti, M., Heida, M., Maccarone, T., $et~{}al$. 2021, arXiv e-prints, arXiv:2112.00339
  * Brightman $et~{}al$. (2020) Brightman, M., Walton, D. J., Xu, Y., $et~{}al$. 2020, ApJ, 889, 71
  * Brightman $et~{}al$. (2016) Brightman, M., Harrison, F. A., Barret, D., $et~{}al$. 2016, ApJ, 829, 28
  * Brightman $et~{}al$. (2019) Brightman, M., Harrison, F. A., Bachetti, M., $et~{}al$. 2019, ApJ, 873, 115
  * Chakrabarti & Titarchuk (1995) Chakrabarti, S., & Titarchuk, L. G. 1995, ApJ, 455, 623
  * Chakrabarti (1989) Chakrabarti, S. K. 1989, ApJ, 347, 365
  * Chakrabarti (1999) —. 1999, A&A, 351, 185
  * Chakrabarti $et~{}al$. (2008) Chakrabarti, S. K., Debnath, D., Nandi, A., & Pal, P. S. 2008, A&A, 489, L41
  * Chakrabarti $et~{}al$. (2015) Chakrabarti, S. K., Mondal, S., & Debnath, D. 2015, MNRAS, 452, 3451
  * Chakrabarti $et~{}al$. (2005) Chakrabarti, S. K., Nandi, A., Debnath, D., Sarkar, R., & Datta, B. G. 2005, arXiv e-prints, astro
  * Colbert & Mushotzky (1999) Colbert, E. J. M., & Mushotzky, R. F. 1999, ApJ, 519, 89
  * Debnath $et~{}al$. (2014) Debnath, D., Chakrabarti, S. K., & Mondal, S. 2014, MNRAS, 440, L121
  * Dewangan $et~{}al$. (2006) Dewangan, G. C., Titarchuk, L., & Griffiths, R. E. 2006, ApJL, 637, L21
  * Du $et~{}al$. (2015) Du, P., Hu, C., Lu, K.-X., $et~{}al$. 2015, ApJ, 806, 22
  * Feng & Kaaret (2010) Feng, H., & Kaaret, P. 2010, ApJL, 712, L169
  * Fiorito & Titarchuk (2004) Fiorito, R., & Titarchuk, L. 2004, ApJL, 614, L113
  * Freedman $et~{}al$. (1994) Freedman, W. L., Hughes, S. M., Madore, B. F., $et~{}al$. 1994, ApJ, 427, 628
  * Gladstone $et~{}al$. (2009) Gladstone, J. C., Roberts, T. P., & Done, C. 2009, MNRAS, 397, 1836
  * Godet $et~{}al$. (2009) Godet, O., Barret, D., Webb, N. A., Farrell, S. A., & Gehrels, N. 2009, ApJL, 705, L109
  * Godet $et~{}al$. (2012) Godet, O., Plazolles, B., Kawaguchi, T., $et~{}al$. 2012, ApJ, 752, 34
  * Greene $et~{}al$. (2020) Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257
  * Harrison $et~{}al$. (2013) Harrison, F. A., Craig, W. W., Christensen, F. E., $et~{}al$. 2013, ApJ, 770, 103
  * Iwasawa (2021) Iwasawa, K. 2021, A&A, 652, A18
  * Jana $et~{}al$. (2016) Jana, A., Debnath, D., Chakrabarti, S. K., Mondal, S., & Molla, A. A. 2016, ApJ, 819, 107
  * Jin $et~{}al$. (2016) Jin, C., Done, C., & Ward, M. 2016, MNRAS, 455, 691
  * Kaaret $et~{}al$. (2017) Kaaret, P., Feng, H., & Roberts, T. P. 2017, ARA&A, 55, 303
  * Kaaret $et~{}al$. (2006a) Kaaret, P., Simet, M. G., & Lang, C. C. 2006a, ApJ, 646, 174
  * Kaaret $et~{}al$. (2006b) —. 2006b, Science, 311, 491
  * Kalberla $et~{}al$. (2005) Kalberla, P. M. W., Burton, W. B., Hartmann, D., $et~{}al$. 2005, A&A, 440, 775
  * King & Pounds (2003) King, A. R., & Pounds, K. A. 2003, MNRAS, 345, 657
  * Kosec $et~{}al$. (2018) Kosec, P., Pinto, C., Fabian, A. C., & Walton, D. J. 2018, MNRAS, 473, 5680
  * Li $et~{}al$. (2005) Li, L.-X., Zimmerman, E. R., Narayan, R., & McClintock, J. E. 2005, ApJS, 157, 335
  * Liu $et~{}al$. (2021) Liu, H., Luo, B., Brandt, W. N., $et~{}al$. 2021, ApJ, 910, 103
  * Matsumoto & Tsuru (2001) Matsumoto, H., & Tsuru, T. G. 2001, in American Institute of Physics Conference Series, Vol. 599, X-ray Astronomy: Stellar Endpoints, AGN, and the Diffuse X-ray Background, ed. N. E. White, G. Malaguti, & G. G. C. Palumbo, 758–761
  * Matsumoto $et~{}al$. (2003) Matsumoto, H., Tsuru, T. G., Watari, K., Mineshige, S., & Matsushita, S. 2003, in Astronomical Society of the Pacific Conference Series, Vol. 289, The Proceedings of the IAU 8th Asian-Pacific Regional Meeting, Volume 1, ed. S. Ikeuchi, J. Hearnshaw, & T. Hanawa, 291–294
  * McCrady $et~{}al$. (2003) McCrady, N., Gilbert, A. M., & Graham, J. R. 2003, ApJ, 596, 240
  * McHardy $et~{}al$. (2006) McHardy, I. M., Koerding, E., Knigge, C., Uttley, P., & Fender, R. P. 2006, Nature, 444, 730
  * Mezcua (2017) Mezcua, M. 2017, International Journal of Modern Physics D, 26, 1730021
  * Miller (2005) Miller, J. M. 2005, Ap&SS, 300, 227
  * Mineo $et~{}al$. (2012) Mineo, S., Gilfanov, M., & Sunyaev, R. 2012, MNRAS, 419, 2095
  * Molla $et~{}al$. (2017) Molla, A. A., Chakrabarti, S. K., Debnath, D., & Mondal, S. 2017, ApJ, 834, 88
  * Molteni $et~{}al$. (1996) Molteni, D., Sponholz, H., & Chakrabarti, S. K. 1996, ApJ, 457, 805
  * Mondal (2020) Mondal, S. 2020, MNRAS, 492, 804
  * Mondal $et~{}al$. (2022) Mondal, S., Adhikari, T. P., Hryniewicz, K., Stalin, C. S., & Pandey, A. 2022, arXiv e-prints, arXiv:2204.02133
  * Mondal & Chakrabarti (2021) Mondal, S., & Chakrabarti, S. K. 2021, The Astrophysical Journal, 920, 41
  * Mondal $et~{}al$. (2014) Mondal, S., Debnath, D., & Chakrabarti, S. K. 2014, ApJ, 786, 4
  * Mondal & Stalin (2021) Mondal, S., & Stalin, C. S. 2021, Galaxies, 9, 21
  * Mucciarelli $et~{}al$. (2006) Mucciarelli, P., Casella, P., Belloni, T., Zampieri, L., & Ranalli, P. 2006, MNRAS, 365, 1123
  * Nandi $et~{}al$. (2019) Nandi, P., Chakrabarti, S. K., & Mondal, S. 2019, ApJ, 877, 65
  * Ohsuga $et~{}al$. (2005) Ohsuga, K., Mori, M., Nakamoto, T., & Mineshige, S. 2005, ApJ, 628, 368
  * Okajima $et~{}al$. (2006) Okajima, T., Ebisawa, K., & Kawaguchi, T. 2006, ApJL, 652, L105
  * Pasham $et~{}al$. (2014) Pasham, D. R., Strohmayer, T. E., & Mushotzky, R. F. 2014, Nature, 513, 74
  * Patruno $et~{}al$. (2006) Patruno, A., Portegies Zwart, S., Dewi, J., & Hopman, C. 2006, MNRAS, 370, L6
  * Pinto $et~{}al$. (2016) Pinto, C., Middleton, M. J., & Fabian, A. C. 2016, Nature, 533, 64
  * Poutanen $et~{}al$. (2007) Poutanen, J., Lipunova, G., Fabrika, S., Butkevich, A. G., & Abolmasov, P. 2007, MNRAS, 377, 1187
  * Qiu $et~{}al$. (2015) Qiu, Y., Liu, J., Guo, J., & Wang, J. 2015, ApJL, 809, L28
  * Reeves $et~{}al$. (2008) Reeves, J., Done, C., Pounds, K., $et~{}al$. 2008, MNRAS, 385, L108
  * Sazonov $et~{}al$. (2014) Sazonov, S. Y., Lutovinov, A. A., & Krivonos, R. A. 2014, Astronomy Letters, 40, 65
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Singh $et~{}al$. (2021) Singh, C. B., Mondal, S., & Garofalo, D. 2021, Monthly Notices of the Royal Astronomical Society, 510, 807–814
  * Stobbart $et~{}al$. (2006) Stobbart, A. M., Roberts, T. P., & Wilms, J. 2006, MNRAS, 368, 397
  * Strohmayer & Mushotzky (2003) Strohmayer, T. E., & Mushotzky, R. F. 2003, ApJL, 586, L61
  * Toyouchi $et~{}al$. (2021) Toyouchi, D., Inayoshi, K., Hosokawa, T., & Kuiper, R. 2021, ApJ, 907, 74
  * Wilms $et~{}al$. (2000) Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914
  * Wilson-Hodge $et~{}al$. (2018) Wilson-Hodge, C. A., Malacaria, C., Jenke, P. A., $et~{}al$. 2018, ApJ, 863, 9