{'fraction_non_alphanumeric': 0.029378273741744478, 'fraction_numerical': 0.01203218704926744, 'mean_word_length': 4.426673532440783, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Bayesian Approach to Modelling Multi-Messenger Emission from Blazars using Lepto-Hadronic Kinetic Equations

Abstract Blazar TXS 0506+056 is the main candidate for a coincident neutrino and gamma-ray flare event. In this paper, we compose a detailed kinetic lepto-hadronic emission model capable of producing a photon and neutrino spectrum given a set of parameters. Our model includes a range of large-scale geometries and both dynamical and steady-state injection models for electrons and protons. We link this model with a Markov Chain Monte Carlo sampler to obtain a powerful statistical tool that allows us to both fit the Spectral Energy Distribution and study the probability density functions and correlations of the parameters. Assuming a fiducial neutrino flux, we demonstrate how multi-messenger observations can be modelled jointly in a Bayesian framework. We find the best parameters for each of the variants of the model tested and report on their cross-correlations. Additionally, we confirm that reproducing the neutrino flux of TXS 0506+056 requires an extreme proton to electron ratio either in the local acceleration process or from external injection.


INTRODUCTION
The modelling of blazar emission and blazar flares is a topic with a long history (see for example Böttcher 2019), yet a number of questions still remain open. Firstly, the nature of the hot emitting plasma remains debated. Early models were mostly leptonic (Jones 1968), in the sense that only electrons and positrons were considered to play a role in the emission. More recently, lepto-hadronic and hadronic models have started to gain a wide acceptance (Mastichiadis & Kirk 1995;Gao et al. 2017). Secondly, the number of relevant emitting plasma zones remains unknown. Currently, the main models are based on the assumptions of a single emission zone (a spherical blob) (Böttcher et al. 2013;Cerruti et al. 2018;Keivani et al. 2018) or two main emission zones (jet and sheath models) (MacDonald et al. 2015;Potter 2018), although some multi zone emission models exist (Graff et al. 2008;Böttcher & Dermer 2010). Thirdly, the plasma heating mechanisms remain poorly understood. Proposed mechanisms include Magnetic Reconnection (Morris et al. 2019;Petropoulou et al. 2019) and Internal Shocks (Spada et al. 2001;Kino et al. 2004;Böttcher & Dermer 2010;Pe'er et al. 2017). Recently, a mayor milestone was achieved in addressing the first problem with the detection of directionally coincident neutrinos and gamma-rays from blazar TXS 0506+056 (IceCube et al. 2018), as only models with a hadronic component are able to explain the presence of neutrinos. However, this in turn has raised an issue with the amount of energy required for heating and keeping the particles hot, as hadronic models require much more energy than leptonic ones. The other questions also remain open. In fact, we can see that at least two of the studies explaining the emission of TXS 0506+056 use a one zone model (Cerruti et al. 2018;Keivani et al. 2018) whereas others use a nested leaky box model (Gao et al. 2019). The nature of the plasma heating mechanism is addressed by none of the studies. Instead, either they inject an already heated plasma or they set a steady population of particles, whose parameters are given based on observational and energetic constrains or determined in a multi-dimensional parameter search.
The methods that these studies use allow them to obtain the best fit parameters according to a grid search or a fit by eye. But they do not allow for a full statistical analysis on the parameters, their distributions and their correlation. This work is an attempt at achieving this objective. This work also builds upon a lepto-hadronic model for all the emission processes inside the plasma with a single zone of emission. We choose that the main particles (electrons and protons) can be in a steady state by default or can be injected. In addition, we also investigate the effects of the geometry of the emitting plasma on the modelling by simulating a sphere and a disk.
On top of the emission model we add a Markov Chain Monte Carlo (MCMC) method that allows us to quickly fit the Spectral Energy Distribution (SED) of TXS 0506+056. This method of obtaining a fit for the SED gives us important statistical information about the parameters and their correlations, encoded in the posterior distribution functions for each of the parameters and correlation plots for each pair of parameters.
This study starts with an overview of the model, the volumes we simulate, the kinetic equations we solve and the approach we take for the initial particles in section §2. We continue with a brief presentation of how we fit the models in section §3. Subsequently, we report our results reproducing TXS 0506+056 in section §4. And finally we draw conclusions in section §5. As cosmological parameters, we have chosen those of Collaboration et al. (2019): H0 = 67.66, Ωm = 0.3111 and ΩΛ = 0.6889.

METHOD
Our approach to simulating SEDs follows four steps. First, we define the geometry of the emitting volume, which we will employ to set the escape timescale and the relation between the photon and neutrino populations and the flux an observer would receive. Second, we set the type of population behaviour for protons and electrons, choosing between leaving them with fixed populations, or injecting particles from outside. Third, we simulate the evolution of the plasma in the volume until all the particles in it reach a steady state.
As a final step, our model takes the population number density of photons (nγ) and neutrinos (nν ) and transforms it to the flux (νFν ) that an observer would receive.
For photons, this transformation is given by (adding a prime to variables in the plasma frame Böttcher & Baring 2019, equation 2.3): where ε is the energy of a photon normalized to the rest mass of an electron, δ is the Doppler boosting factor, DL is the luminosity distance, V is the volume that we are simulating, τ esc is the escape timescale for photons and c is the speed of light.
To simplify the notation, we define S f as the factor between ε 2 n γ (ε ) and νFν : The neutrino transformation from the obtained population to the theoretically observed flux is similar to that of photons, but we need to take into account the flavour of the produced neutrinos and neutrino oscillations. In order to do this, we assume that the neutrinos have had time to be distributed equally among the three flavours.
Note that to simplify subsequent formulas, and because we usually work in the plasma frame, we will drop the primes.

Simulated volumes
We consider two emission region geometries in our modelling. A sphere and a thin disk. The former is a first approximation to an unknown shape of the emitting plasma whereas the latter tries to take into account a possible shape derived from some heating mechanisms, such as magnetic reconnection and shock acceleration. In figure 1 we show a simple sketch of both volumes. It should be noted that one of the basic assumptions of blazar modelling in general is that the emission zone is homogeneous and isotropic as all the equations for particle interactions assume isotropy at the interaction point. Although the effect is smaller for a sphere than a disc model, any choice of geometry will introduce an inherent anisotropy, but it is beyond the scope of this work to study this effect..

Sphere
The Spherical model (figure 1 left) is defined by one characteristic length: the radius of the blob R. The free escape timescale, that is, the time a particle that does not interact in the volume takes to escape, (τ f ree,esc ) is then set as the average length that a particle would have to traverse to escape from a sphere divided by c: and the relation between ε 2 nγ (ε) and νFν can be then obtained as:

Thin Disk
The Thin Disk model (figure 1 right) is defined by two characteristic lengths: the radius of the disk (R) and the height of the disk (h). The first order approximation to the escape timescale is: and the relation between ε 2 nγ (ε) and νFν can be then obtained as: It is interesting to notice that S f does not depend on the height of the shell, as that dependency goes away when calculating V τesc .

Kinetic equation approach to emission
The model we employ solves the coupled kinetic equations for 12 different species of particles: protons, neutrons, electrons, charged and neutral pions, muons and antimuons, and electron and muon neutrinos and antineutrinos. The general form of a kinetic equation, omitting the linear term follows: As shown in Equation (7), the variation of the population of a species of particle is due to: external injection Qe; internal injection Qi, which represents particles produced by processes internal to the simulated volume; loss of particles L due to interactions with other particles; an escape term where τesc represents the timescale for escaping the simulated volume, which we take as a multiple of the free timescale (τesc = στ f ree,esc ); a decay term where τ dec represents the decay timescale; and a continuous loss of energy where the energy lost is proportional to γ 2 , such as synchrotron losses for charged particles and the Inverse Compton (IC) process for electrons.
Usually, the expressions for photon injection and losses are given in terms of emission and absorption coefficients (jν and αν ). For charged particle losses, they are given in terms of an energy loss. These are related according to: where is the reduced Plank constant.
In our model, all the charged particles create synchrotron emission, have synchrotron losses and cause synchrotron absorption. For this, we model the emission as: where Pν (γ) is the energy emitted by a particle with energy γ at frequency ν and has the form (Crusius & Schlickeiser 1986): where q and m are the charge and the mass of the emitting particle, B is the magnetic field strength, x is the frequency normalized to a critical frequency x = ν/νc with νc := ν0γ 2 := 3 4π and CS (x) is a function given in terms of Whittaker functions: The loss of energy of the emitting particles is modelled as: where σT is the Thomson cross-section and uB is the energy density of the magnetic field. The absorption of photons is modelled as: Photons interact with energetic electrons due to the inverse Compton process, which we model as: where ε is the energy of a photon normalized with respect to the energy of a rest electron (ε is the initial energy of the photon and εs is the energy after the scattering), ne is the electron population, nγ is the photon population and g u,d (γ, ε, εs) are two functions referring to the upscattering and downscattering of photons which where calculated by Jones (Jones 1968, equations 40 and 44): if qu ≤ 1 ≤ εs ε , εs < γ, The loss term for photons can be obtained as: where R (γe, ε) is the reaction rate of electrons of energy γe with photons of energy ε and can be found from integrating the Klein-Nishina cross section. where: and Li2 the dilogarithm. Following Böttcher et al. (2012), we use the simplified prescription for the electron energy losses in inverse comptom scattering: where uγ is the energy density of the photon field: Photons can interact among themselves creating electron and positron pairs, which we model according to (Böttcher & Schlickeiser 1997). Photon losses are modelled according to the fully integrated cross section in a similar way to the photon loss term for IC.
Hadrons (protons and neutrons) interact with photons producing hadronic resonances which ultimately decay in the form of hadrons and pions. To model this we follow the approximations of Hümmer et al (Hümmer et al. 2010) in their SimB model. We diverge though in the calculations of the loss and self-injection term for hadrons as, instead of applying a cooling term, we model every process as a loss and where applicable, we reinject the corresponding hadron with reduced energies.
Neutrons, charged pions and muons decay, producing ultimately protons, electrons, positrons and neutrinos. To model these processes, we follow Hümmer et al (Hümmer et al. 2010) for the decay of pions and Lipari et al (Lipari et al. 2007) for the decay of muons and the production of neutrinos.
Finally, we do not model other processes such as protonproton interaction, heating of charged particles when absorbing photons and heating of electrons when downscattering photons because their impact is negligible relative to that of the processes included in our model.

Steady State Models
To show a simple example of the capabilities of our model, we first use it in a steady state model to fit the SED of TXS 0506+056. In our steady state model the populations of both electrons and protons do not evolve over time. This also covers a scenario where both are injected constantly over time with a fixed distribution that would lead to the same steady state distribution that we here fix from the outset.
The initial distribution of electrons is considered to be a broken power law and for protons a simple power law: if γe,min < γe < γ e,break Bγ −p 2 e if γ e,break < γe < γe,max , where A, B and C are factors that normalize the distributions and can be calculated from the conditions: γp,max γ p,min and we have: We can assume that the ratio of protons to electrons is given by η. This means that, given a combined mass density ρ of accelerated electrons and protons, the steady populations for electrons and protons are: if γe,min < γe < γ e,break γ p 2 −p 1 e,break γ −p 2 e if γ e,break < γe < γe,max , Initially, in the simulated volume there are only electrons and protons, so we need to evolve the plasma until it reaches a 'true' steady state for all the particles. This is detected by calculating a statistic related to the relative change in population of particles with time: where n is the kind of particles we are interested in. In our case, we test three particles species for convergence: photons, electrons and neutrinos. We consider that the system has reached the photonic steady state when the photon statistic has a value lower than 10 −8 , and the final steady state when all three have a value lower than 10 −8 . The value 10 −8 has been chosen as a good compromise between the accuracy of the steady state and the simulation time, given that the system only asymptotically approaches a true steady state. We check the photon, neutrino and electron populations for convergence for the following reasons. The photon population directly shapes the SED, the electron population has a long covergence timescale especially at low energies due to cooling, and neutrinos are the last particles to be produced and therefore formally the last to achieve steady state.
The steady state for protons and electrons is implemented by choosing not to update their populations during the evolution of the system. But we can use this fact to have an a priori value for the energy density in protons and electrons: From the energy densities, we can also derive an estimate for the luminosity of the plasma as:

Dynamic Models
Dynamic models break the assumption of a steady state population of electrons and protons by allowing them to evolve over time. To arrive at a steady state population of particles we have two options: we can inject particles arriving from outside the volume with a given distribution, or we can keep a steady population of cold particles that are accelerated to higher energies within the volume. In this work we focus on the first option. There are various reasons for this. One is that acceleration introduces various parameters, the acceleration timescales and mechanisms, which are not well understood. Second, the accelerated particles have to come from a population of cold particles, which we do not model. Third, acceleration requires a bit of tuning to obtain the right slope and enough maximum energy for the particle populations. Fourth, one can argue that acceleration is only going to happen at the boundaries of the simulated volume (e.g. shock fronts), and thus, its effect is to inject energetic particles with a given distribution to the rest of the volume. Particle injection is characterized by a luminosity (L) and a normalized particle distribution f (γ) such that: where V is the volume of the plasma where the injection is happening, m is the mass of the injected particle and Q0 is the number density of injected particles. For a known particle distribution, we can solve for Q0 as: where γ is the average energy of the injected particles. In our case, we inject electrons and protons. Assuming that the ratio of injected protons to electrons is η, we have: and thus: Once a steady state is achieved and assuming that the injection follows a broken power law, simplified versions of the kinetic equations can be applied to obtain an approximation to the final electron and proton populations.
For low energies, we have that the cooling time due to synchrotron and Inverse Compton processes is longer than the escape timescale of the particles. This means that we can approximate the population evolution as: By taking the derivative equal to 0, we obtain that in the steady state limit the population is: If a steady state is achieved from dynamically injected protons and electrons, the dynamic model and the steady state model relate according to: For higher energies where we can not ignore the effects of synchrotron cooling and Inverse Compton losses, we can approximate the population with a power law with a given index p. We obtain: In the steady state, we obtain: For (2 − p) γSτesc 1 we recover the previous result where we ignored the contributions of Synchrotron and Inverse Compton cooling. When (2 − p) γSτesc 1 we obtain that the main driver of the population is cooling, and we get: where we see that the effect that cooling produces is to increase the index of the injected distribution of particles by 1.
Particles with energies lower than the minimum energy of the injected energy can appear due to cooling. If the turnover energy in a charged particle population due to cooling lies below the level associated with the lower injection cut-off Lorentz factor, an intermediate regime appears in the particle spectrum. In the absence of additional physical processes, this would lead to a power law slope of 2. Below the turnover energy, escape takes over from cooling as the main driver, which leads to a drop in the population.

BAYESIAN AND MARKOV CHAIN MONTE CARLO ANALYSIS
Model fitting can be done in various ways. An easy and straightforward way of doing it is by eye (Keivani et al. 2018), where the goodness of the fit is estimated by seeing if the simulated spectrum fits the observations. The disadvantages of this approach are that the validity of the fit can not be validated numerically and that, depending on the parameters, finding a good fit might be a complicated task. Another way is by performing grid searches (Gao et al. 2019), where the parameter space is divided into a grid, and, for every point in this grid you assign a goodness of fit value, taking the best parameters as the final ones. This has the advantage of being fully automated but has the disadvantage of requiring a high amount of computations, as quite a lot of time is spent calculating the likelihood of bad parameters.
A third method, similar to the latter, is to perform an MCMC search over the parameter space (see for example Qin et al. (2018) for a study that employed MCMC to fit simple leptonic models to blazars). In this case, instead of dividing the parameter space into a grid and checking every point, a series of 'walkers' traverse the parameter space checking the likelihood of the fit at the point where they are. Their advance is biased towards improving the likelihood so as to guide their search towards 'good' parameters. Nevertheless, the search is ergodic, which means that with time, the whole parameter space will be traversed. Just like the previous method, MCMC searches are fully automated, but they have the advantage of being focused on the best parameters. This allows us to obtain two results besides the best parameters. On the one hand, it allows us to obtain a measure of the error in each of the parameters. On the other hand, a study of the behaviour of the MCMC walkers allows us to identify correlations between the parameters. This in turn allows us to redefine our parameters so as to eliminate the correlations between them and improve the fits. It should be noted that not all MCMC sampling algorithms perform worse under correlated parameters. Nevertheless, finding good 'fit parameters' a priori can give a very useful physical insight in the interconnectedness of the model parameters.
Therefore, we will refer to two sets of parameters: 'model parameters', or simply parameters, which are the input of our model; and 'fitted parameters', which are the set that MCMC fits. Both sets are related by invertible transformations, so it is easy to go from one to other.
MCMC requires some information from our part to perform its search. Just as grid search needs a grid to be defined, MCMC needs information about the distribution of the parameters it has to fit. The priors on the fitted parameters are derived by starting from flat priors in either linear or logarithmic space of the model parameters. The distributions are chosen so that parameters that enter the model as multiplicative terms are defined in logarithmic space and parameters that enter as exponents are defined in linear space. The limits of the distributions are taken very loosely, even to the point where in some cases we allow unphysical values. The only hard limits are the relations γe,min < γ e,break < γe,max and p1 < p2, and in the case of thin shells, h < R.
The log-likelihood of the model is calculated as: where yn are the observations,ŷn are the predictions of our model and σn are the errors of the observations, which we assume to be a few percent of the observations, depending on the observed band. In our case, we take that for optical frequencies, we have σn = 0.1 yn; for X-rays σn = 0.2 yn; and for gamma-rays σn = 0.5 yn.
The log-likelihood of the 'fitted parameters' is then defined as a combination of the log-likelihood of the model plus the a-priori log-likelihood of the 'fitted parameters'. Denoting the 'model parameters' as P the 'fitted parameters' as F and the a-priori log-likelihood of a set of parameters as p (X), we obtain: where the term |dF/dP | represents the Jacobian of the transformation between parameters, and p (P ) = 0 by construction.
In our implementation, we use an affine invariant sampler (Goodman & Weare 2010) implemented in the Python package emcee (Foreman-Mackey et al. 2013) using 128 walkers. The total numbers of samples taken differ between runs and are reported along with the outcomes of a run.
The first data points of each walker correspond to the period of 'burn in', where the walkers spread from their initial conditions until they explore the log-likelihood function in an uncorrelated way. During each of the runs, the data is binned and the values of the average, median, standard deviation and 16 and 84 percentiles are checked for trends. After they remain consistent for several bins we consider that the run has found the final distribution and the run is stopped.
The bins where the aforementioned values drift are considered the 'burn-in' phase and their data discarded for the final study.
As an additional consistency check for the remaining data. We rebin the data and plot the PDFs to verify that the features of the distributions remain essentially unchanged.
While it is never formally possible to rule out trends with an extremely long correlation time, we are therefore confident that the results do not contain biases.
All the results reported in tables correspond to the median and the 16 and 84 percentiles.

APPLYING OUR MODELS IN A BAYESIAN FRAMEWORK
Before trying to fit a data set, we do an a priori analysis of the physics of the emission processes in order to obtain constraints and correlations for our parameters. These relations are then used to find good initial parameters for the MCMC analysis and to interpret the outcome of the fitting process.

Initial Assumptions and Constraints
Figure 2 shows a simple schematic view of an SED, which we can use to obtain a series of constraints and first guesses for our search: • We can assume that the first peak of the SED is near the optical band, and that it is caused by the synchrotron emission of the electrons at the break point γ e,break . This gives us the constraint: • We know that we have to obtain photons of a certain energy, and that those come from the Inverse Compton process with high energy electrons. So, we can derive a lower limit for the maximum value of γe,max: ν obs,max γe,max mec 2 2π • A second lower limit for the value of γe,max, albeit much less limiting, can be obtained from the X-ray valley. We can safely assume that the emission for the low energy part of the valley comes almost entirely from the synchrotron emission of high energy electrons. This gives us the constraint: • We can fit a line between the optical and the X-ray emission. As it is above the synchrotron peak, we can assume that this emission comes from electrons with γe > γ e,break and, as the theoretical slope of synchrotron emission is known, we can find the constraint: • We can fit a line between the X-ray emission and Gamma-ray emission, and we can assume that the emission in this range is going to be caused by the upscattering of synchrotron photons by the Inverse Compton process. The inverse Compton process increases the energy of photons as ε af ter ≈ γ 2 e ε bef ore . For photons from the synchrotron peak, upscattered by electrons at γ e,break , their energy is above the fit. So, these photons should come from synchrotron photons emitted by electrons with γe < γ e,break . As the Inverse Compton process does not modify the slope of the emission, we obtain: • The radio observations can be considered as upper limits for our model, as the flux can come from the extended emission in the jet (see for example Cerruti et al. 2018). From this fact, we can derive another constraint. Assuming, as before, that the optical observations are the peak of the synchrotron spectrum, and that they are caused by electrons with energy γ e,break , we can extrapolate the observations down to the frequency associated with γe,min: so: At lower frequencies, this can be extrapolated further using the index 4/3. And we obtain: so: (νFν ) radio = (νFν ) break νmin ν break α 1 ν radio νmin Now, all the frequencies can be transformed into associated gamma factors for electrons in the fluid frame: (νFν ) radio = (νFν ) break γe,min γ e,break 2α 1 γ e,radio γe,min Note that γ e,radio is not a real value for the gamma factor of electrons, as it can be < 1. Also note that it might be the case that the synchrotron self absorption frequency is higher than the frequency associated with γe,min. In this case, by using γe,min we would be overestimating the approximation of (νFν ) radio . The constraint is that the observed flux in radio has to be less than the calculated, so: (νFν ) radio,obs (νFν ) break,obs γe,min γ e,break 2α 1 γ e,radio γe,min • The energy of the photons at the optical band is enough so that the blob can be considered thin to the radiation. In this case, the main loss term for photons is escape. Under the simplification of a δ-approximation to synchrotron emission (Böttcher et al. 2012, Chapter 3) 1 , which we only employ here in order to make an estimate, we have: The emission coefficient is related to the rate of injection of photons as: The δ-approximation also allows us to have that ν/ν0 = γe, and thus, we can relate ν, ε and γe as: which means that: As the main losses for photons is the escape term we have: In the steady state, ∂n/∂t = 0, so, we can solve for the equilibrium population of photons as: The electron population is related to the density of the blob and we have, assuming that γe,min γ e,break : for γe ≤ γ e,break .
The transformation from population to observed flux requires us to multiply the previous expression by ε 2 , which we do using the relation with γe: Doing the final transformation we have: It is important to remember that this is only valid for energies where the blob is not thick to the radiation, and only for γ < γ e,break .
• Following the previous point, we can also derive the flux for the low energy part of the x-ray band, as we assume that those photons are synchrotron emission from high energy electrons. The population of high energy electrons is: Then, the photon population for the high energy part of the synchrotron spectrum is: for γe > γ e,break . So, the X-Ray flux is: ∝ V δ 4 B 2 γ 3 e,x−ray ρ (p1 − 1) γ p 1 −1 e,min γ p 2 −p 1 e,break γ −p 2 e,x−ray . Following one of our previous constraints, γe,x−ray γe,max, so: Dividing our previous prediction for optical flux and X-rays flux, we have: • Photons from synchrotron processes get upscattered by high energy electrons due to Inverse Compton processes. As an approximation, and just so we can derive a functional form for the photon population at high energies, we can say that the emission coefficient for IC photons depends on the population of photons at low energy and the population of electrons at high energy. As before, we assume that the energy of a photon is increased as ε af ter ≈ γ 2 e ε bef ore .
Transforming to the rate of injection of photons, we obtain: Again, we assume that the main losses are the escape process, so we can calculate the population of gamma ray photons in our blob as: To compare with the optical flux, we just multiply by ε 2 : We assume that the photons that get upscattered to γrays are those from the synchrotron peak, and that the electrons with γ e,break are the ones mainly responsible for the upscattering. This means that εoptγ 2 e,break = εγ−ray and nγ (εopt) = nγ εγ−ray/γ 2 e,break . This means that the ratio of fluxes is: ∝ 1 δB τescρ (p1 − 1) γ p 1 −1 e,min γ −p 1 e,break .
• Dividing the optical flux reduced by δBγ 2 • And reducing once more by δBγ 2 e,break : The estimations and relations that we have derived in the previous subsection, and the numerical estimates that we have obtained in this can then be used in two ways. First, the relations can guide our search for parameters that are usually together and group them into the 'fitted parameters' for the MCMC fitting procedure. And second, the estimates allow us to give good estimation of the initial value of a fitted parameter, if we happen to chose one of them, or at least provide boundaries for combinations of the initial parameters. At the start of our fitting process, we initialize our walkers in a small sphere, taken to be much smaller than the allowed range for all the parameters, around the estimates.

Other Settings and Results
Our model includes several parameters that set the physics and numerics of the emitting plasma. Some of these are fitted for, and are summarized in table 1 along with their allowed ranges. Among these parameters, the electron distribution is taken to be a broken power law (p1 and p2 refer to the slopes across this break). Other parameters are kept fixed and are discussed below.
Following other authors (Gao et al. 2019;Cerruti et al. 2018), we assume for the accelerated proton population an injection with a single power law in energy with decay slope of 2. The lower cut-off value of the proton population γp,min does not have a noticeable impact on the final fit results and has been fixed at 10. For the value of γp,max we have chosen a reasonable energy that allows for the creation of neutrinos with the energy detected by Ice Cube. Above the threshold of 5 TeV, the estimated energy deposited in a pion due to photon-proton pion production in the best case of multi pion production is ∼ 0.6 of that of the parent proton (Hümmer et al. 2010). Then, through the muon decay channel, the neutrino can obtain almost all the energy of the parent pion. The energy of the detected neutrino is ≈ 290 TeV (Ice-Cube et al. 2018), which gives a value for γp,max ≈ 10 5.71 . In order to allow for other neutrino production channels to contribute, we choose to increase this value to γp,max = 10 6 , although as we will show in section 4.6, this value does not affect the predicted neutrino flux at 290 TeV.
The escape time for charged particles will be longer than that for uncharged particles. Following Gao et al. (2019), we take this ratio of escape times, σ (see section 2.2), to be equal to 300.
In our dynamic models where the proton and electron populations are not initially present, we set both number densities to numerically small values at the start.
We set the range of the numerical grid that stores the electron distribution broad enough to safely encompass the emergent electron distribution from injection, using a lower cut-off Lorentz factor of 2 and an upper cut-off value that is set to max γe,max, mp me γp,max , where γe,max and γp,max are the upper cut-off values respectively of the electron and proton injections. The lower and upper boundaries of the photon population grid have also been chosen generously to account for all possible photon production channels. These have respectively been set at 10 −12 and equal to the electron upper grid boundary.
In principle we can sidestep the issue of charge neutrality in the plasma by allowing the low-temperature pools of electrons and protons to compensate for any charge imbalance in the injection terms. However, assuming the injection to be either acceleration of these low-temperature populations (e.g. at boundary shocks), or to originate from the same external region, we have no obvious physical mechanism that suggests a larger number of protons than electrons to be accelerated (if anything, the opposite), nor for the occurence of a larger pool of protons than electrons in the source of the injection, and our default approach therefore assumes η = 1. We explore the implications of relaxing this constraint in section 4.6.
This still leaves us with all the parameters of Table 1 to fit. Figure 3 is a combination of the best SED fits (tables 2, 3, 4 left and right) that we can obtain with our model given the parameters that we set a priori. It can be appreciated that all the models can reproduce the SED very well, but they fall short by ∼ 3 orders of magnitude in the neutrino flux. We discuss this point in the next subsection of our work. Figures C1, C2, C3 and C4 show the corner plots that we obtain from the MCMC fitting. The blue lines denote the parameters that give the best SED fit.

Spherical Steady State
The spherical steady state model has 9 parameters: (γe,min, γ e,break , γe,max, p1, p2, B, δ, R and ρ), but we know from our analysis of section 4.1 that some of them are always grouped. Thus, we choose to redefine our parameters in the following way: ν br := δBγ 2 e,break , ρ := ρ, As we showed in equation (59), we need the Jacobian of the transformation to properly calculate the log-likelihood of the model parameters. This is given by: The best SED that we can achieve is shown in figure 3 where we can see that it fits quite well to the available data.
The corner plot of the fit can be seen in figure C1.

Thin Disk Steady State
The thin disk steady state model has 10 parameters: (γe,min, γ e,break , γe,max, p1, p2, B, δ, h, R and ρ). And   again, we choose to redefine our parameters in the following way:  Again, we need the Jacobian of the transformation, which is given by:
Unfortunately, due to equation (54), the parameter that naturally would take the place of N1 in the steady state models depends on the value of S, which is not known a priori. Even if it were known, it separates two different behaviours for the particles in two energetic regimes. This makes obtaining simple relations for the synchrotron, x-ray and gamma-ray fluxes a not trivial problem that, although solvable, does not help the fitting. Because of this, we have chosen to not redefine our parameters and fit directly the model parameters.

Thin Disk Injection
The spherical injection model has 10 parameters: (γe,min, γ e,break , γe,max, p1, p2, B, δ, h, R and L). Again, for the same reason as for the case of spherical injection, we choose to not redefine our parameters.

First Interpretation of Results
First, we check the validity of some of our initial estimates for the initial fitted parameters. For the value of log 10 δBγ 2 e,break we had estimated 8.50. Figure 4 reflects that, with the exception of the sphere injection model, where we obtain a value ∼ 3 times bigger, our estimate was correct. A possible explanation for this minor discrepancy lies in equation (54). This equation sets the energy which separates the regimes where particles mainly escape from the regime where particles mainly cool. Spherical models have typically higher values of R than disks for h, which in turns means higher values for τesc. Then, everything else being equal, we know that γ e,break,sphere < γ e,break,disk .
Second, we check the validity of our initial estimate for log 10 δ 2 V τesc γ −3 e,break , where we had obtained a value of 35.39. Note that the factor V /τesc depends only on R 2 regardless of the geometry. Figure 5 shows that we systematically obtain results below the estimate by a factor ∼ 3. Part of this discrepancy can be explained by looking at the gamma-ray peak of figure 3. We can see that our model obtains a gamma-ray flux a bit above the data points, although well within the error bars. So, the modelled flux is above that that we used for the estimate, and thus its value should be reduced. The other part comes from the physics used in the estimation, which is based on a crude approximation of inverse Compton upscattering. Gao et al. (2019) give an estimate for the Eddington luminosity of the black hole at the center of TXS 0506+056 assuming a mass similar to that of M87: 10 47.8 erg · s −1 . As one of our parameters for the injection models is precisely injected luminosity, we can directly compare with this value. Figure 6 shows that the distribution of values we obtain for spherical injection are well below the estimated one, whereas the values for disk injection are above it.

Interpretation of Posterior Distribution Results and Cross-Correlations
As discussed in 3 when we presented the procedure of fitting with MCMC, one of the products that we obtain is the posterior density function for each of the parameters. And that, in turn, can provide interesting information, for example about how individual parameters impact a model or whether or not bimodalities occur across the range of good fits. One example of the former are the posterior density functions for γe,max. We have not found examples of bimodal distributions. The four γe,max posterior distributions clearly show a lower cut off at low energies and a flat shape at the rest of the energies in the explored parameter space. The lower cut-off can be explained by the constraints of equations (61) and (62), whereas the flat shape can be explained by realizing that our model is not constraining at this range. In fact, an increase of γe,max requires a number of particles so small that it can barely modify ρ or L.    Another of the products of using MCMC are the correlation plots between the parameters, where we can clearly see the relations between the parameters and how changing one affects others. Good examples are the triad γ e,break , δ, R or the correlation between ρ and R or h depending on the geometry of the emitting volume for steady state models.
The correlations between γ e,break , δ and R can be explained by equation (101). V /τesc is a function of R 2 regardless of the geometry of the model, which is the cause of the cross-correlations.
The correlation between ρ and R or h involves the es-cape timescale τesc. As stated in the previous paragraph, V /τesc is only a function of R 2 , which means that in (2) the only geometric factor is R 2 . But, from our exploration of both synchrotron emission and inverse Compton emission, we found that the population of photons in the rest frame of the fluid depends on τesc. This means that the population depends directly on R or h depending on whether the geometry is an sphere or a disk.
A more detailed analysis of the posterior density functions for the case of the spherical injection model reveals  something interesting. In the case of the distribution for p1, we find that its distribution is much more extended than the ones for the other three models. Moreover, looking at the distributions for γe,min and γ e,break we find that their ranges overlap and, as can be seen in figure 7, if we plot the distribution for both combined we find that we obtain a single peak. We infer that the first part of the broken power law is not relevant for the fit, as it is very narrow and its slope can vary wildly.
To check the validity of this reasoning. We have chosen to also test the case of the spherical injection but with a simple power law.

Spherical Injection with a Power Law
The spherical injection with a power law model has 7 parameters: (γe,min, γe,max, p, B, δ, R and L) where we again have chosen to not redefine them. Table 5 corresponds to the parameter values that we have found, figure 8 shows the SED corresponding to the best parameters, and figure C5 shows the corresponding corner plot of the model.
As we can see in figure 8, the fit that we obtain is very similar to that of the broken power law model. We conclude that a simple power law is indeed sufficient to reproduce the observed SED for an injection model in spherical geometry. Figure 3 shows the comparison between the SED we obtained with the four models and the observed fluxes. For η = 1 (equal injection numbers for protons and electrons), the SED can be produced very well, but the neutrino flux remains orders of magnitude below that allowed by the tentative neutrino detection. Below we discuss two mechanisms to increase the neutrino flux.

Large proton to electron ratio injection
One option is to shift the value of gamma p,max. As discussed in section 4.3, a value of 10 6 is enough to produce neutrinos with the appropriate energy through more than one channel. Nevertheless, higher energies allow for a greater energy range for protons, pions and muons to produce neutrinos with ∼ 290 TeV. Additionally, an increase of the maximum energy for protons also increases the total number of protons that have energies above the threshold to produce neutrinos at 290 TeV, thus leading to an increase of the total neutrino flux. For example Keivani et al. (2018) uses γp,max = 1.6 × 10 6 − 3 × 10 9 and Cerruti et al. (2018) has γp,max = 6 × 10 7 − 2.5 × 10 9 .
The second option is to increase the ratio of protons to electrons, both in the steady state population and in the injected population. The main consequence that a higher value of η would have is an increase of the number of protons at all energies, including those with the appropriate one to produce neutrinos. Thus, we will obtain a higher neutrino flux at all energies. For example Gao et al. (2019) uses η ≈ 8606 (as we show in Appendix B) and Cerruti et al. (2018) has η ≈ 200 Appendix A. Figure 9 shows a simple test of these two hypotheses with our model. As a base we have taken the best parameters for the spherical injection model. Then, to check the first hypothesis we have increased the value of γp,max to 10 9 and to check the second we have increased the value of η by   a factor 1750 and modified the value of L according to equation (45) so that the value of injected electrons is the same in the two models. As we predicted in the previous paragraphs, in the first case we see an increase in the total flux of neutrinos, specially those with energies higher than 500 TeV, whereas in the second case the increase is restricted only to those of energies around 290 TeV. Following figure 9, and combined with the simplicity of the spherical injection with a single power law model, we choose to go further in our modelling and allow the value of η to vary. The neutrino flux is an upper limit, so we do not have an actual data point to compare our model results, but for illustrative purposes, we have chosen to act as if the upper limit was an actual detection and use its value in the calculation of the log-likelihood.
The model that we fit has 8 parameters: (γe,min, γe,max, p, B, δ, R, η and L) Table 6 corresponds to the parameter values that we have found, figure 10 shows the fit that we obtain with the best parameters and figure C6 shows the corresponding corner plot of the model. Figure 11 plots the resulting proba-bility density function for η, showing a broad distribution peaking around log 10 η = 3.48. Below η of ∼ 10 3 , the neutrino flux will be underpredicted, whereas for eta approaching 10 4 it becomes difficult to avoid overpredicting the X-ray flux in the SED valley.
It can be seen from figure 10 that both the SED and the neutrino flux can be fitted at the same time without overwhelming the x-ray flux. But this requires a value the ratio of protons to electrons of η ≈ 3388, which we can not explain.

SUMMARY AND CONCLUSIONS
We present a detailed kinetic lepto-hadronic emission model capable of producing photon and neutrino spectra for a range of parameters, plasma geometries and electrons and protons either injected or at steady state. This model has been developed independently of other approaches (e.g. Cerruti et al. (2015) or Gao et al. (2017)), and can be used as a cross-check on their work. We have coupled our model   to a Markov-Chain Monte Carlo method, and discuss the characteristics of the parameter space that this allows us to probe in a novel manner. The method is applied to the multi-messenger case of TXS 0506+056.
MCMC sampling reveals the full probability distribution of the fit parameters as well as their cross-correlations and degeneracies. We demonstrate that the six iterations of our model (between sphere and disc geometry, steady state and injection, simple power-law and varying η) are all able to fit TXS 0506+056. The MCMC method shows that a single power law injection of electrons is sufficient in spherical geometry, that γe,max is poorly constrained and that some of the fitted variables are correlated.
To reproduce a high level of neutrino flux comparable to the peak value allowed by the tentative detection by Ice-Cube requires a large proton to electron ratio. Our MCMC method quantifies the probability distribution of this ratio eta, and we find log 10 η = 3.48 +0.26 −0.28 . If the accelerated protons and electrons are not taken to be from different origins altogether, this poses a challenge to models of charged particle acceleration to explain the origin of an overwhelmingly larger number of accelerated protons than electrons.

APPENDIX A: NOTES ON REPRODUCING THE WORK OF CERRUTI ET AL.
As a way of checking our model, we try to reproduce the results of Cerruti et al. (2018). Their study has both hadronic and lepto-hadronic models and their parameters and the goodness of the fit are available in the extra materials. We choose to reproduce the parameters for a lepto hadronic model that give the best χ 2 value.
Cerruti et al.'s model is based upon a stationary particle distribution with the shape of a broken power-law with an exponential drop which produces the primary photon emission. Although for this work, they set γ e,break = γe,max, obtaining a simple power law. In addition, they calculate the lepton population that would result from photon-photon pair production and bethe-heitler cascades and subsequently add their own emission to the primary emission. For the neutrino emission, they calculate the resulting emission from pion and muon decays.
This hybrid approach to the calculation of the lepton population poses a problem for our model as we can either have a single steady state population or a fully dynamic one. Moreover, as we will see when we derive the corresponding parameters for our model, these can easily cause that one of our assumptions, that of the lepton injection due to pair production and bethe-heitler processes being small, to be broken. With this in mind, we choose to employ the spherical steady state version of our model, knowing that it will not be able to fully reproduce the shape of the observed SED.
As said, their model defines the steady populations for electrons and protons as: Ne,p (γ) = K e,p γ −αe,p exp − γ γe,p,max .
Their population parameter K e,p is given in term of the normalizations for the electron distribution at γ = 1 and the ratio with protons (r). This means that: So: Our model defines the population parameters differently but we can transform between the two as: So: where A and C are normalization constants for the particles distributions. Table A1 shows the results of the translation of the parameters between Cerruti et al.'s model and ours. There are three parameters that will have a very marked effect on the resulting SED: γe,max, γp,max and η. The first (γe,max) effectively sets the maximum energy of synchrotron and Inverse Compton photons, and as it is two orders of magnitude lower than ours, we can predict that we will not be able to fit neither the high energy end of the synchrotron peak, nor that of the Inverse Compton peak. The second (γp,max) contributes to increasing the amount of protons available to interact with photons in both, pion production and bethe-heitler pair production. This will cause the amount of injected leptons and neutrinos to increase, and their energies will be bigger. The third (η) has a similar effect to γp,max, but increases the number of protons at all energies. This means that not only the pion and lepton production will be higher, but also the synchrotron emission. Figure A1 shows the results obtained by our model using Cerruti et al.'s parameters. We have chosen to separate the simulated SED in the total range and a trusted range because, as we have explained, we can not be sure of the emission that would come from electrons with energies above that of injection. Still, for the trusted range, we obtain an acceptable fit with the observed flux levels. Although the neutrino flux is still lower than the estimated by Ice Cube at the relevant energies.

APPENDIX B: NOTES ON REPRODUCING THE WORK OF GAO ET AL.
A second check that we can do for our model is to try to reproduce the SED of TXS 0506+056 employing the parameters of the lepto-hadronic model of Gao et al.'s hybrid model. It must be noted that in their model, they separate the luminosities of protons and electrons, whereas in our model we combine them into one luminosity and regulate the particle populations by means of the ratio of protons to electrons parameter.
By separating Ke,pγ −αe,p into a normalized distribution and a constant, as we did in section 2.4 we easily find: Qe,p,0 = Le,p γe,p me,pc 2 .
The ratio of injected protons to electrons is then: Another difference between our model and theirs is that they take that the charged particles do not escape in the free escape timescale, as they assume that protons, electrons and positrons escape with a velocity of c/300. We have chosen to simulate this by increasing the free escape timescale of charged particles by 300 which has the same effect. Another difference lies in the definition of the free escape timescales,  where we have included a factor 3/4 due to geometrical reasons (3): tesc,our = 3 4 tesc,gao.
In table B1 we list the original parameters and our choices.
B1 shows the resulting SED and neutrino flux that we obtain applying our model with the parameters of (Gao et al. 2019). It is clear that, although the shape is correct, we are not able to correctly reproduce the flux level, being our prediction lower by ∼ 1 order of magnitude. It should be noted that their model is, in fact, a two zone model, where a sphere is inside another one, and here we have chosen to reproduce the emission of the inner one. Additionally, it is important to remember that our model neglects the effects of pair production and Bethe-Heitler processes, as for the parameters we employ they are not important. But, just as in the previous case, Gao et al.'s parameters require the ratio of protons to electrons to be high, which breaks our assumption about these two processes. Moreover, this injection of electrons and positrons could be above that from the external injection, making the whole SED to increase.