\section{Introduction}
\label{sec:intro}



Central to radiative modelling of pulsar wind nebulae (PWN) is the numerically solution the diffusion-loss equation for electron-positron pairs,
\begin{equation}
\label{eq:difflos}
\frac{\partial N(\gamma,t)}{\partial t}=-\frac{\partial}{\partial \gamma}\left[\dot{\gamma}(\gamma)N(\gamma,t) \right]-\frac{N(\gamma,t)}{\tau(\gamma,t)}+Q(\gamma,t),
\end{equation}
where $N(\gamma,t)$ is the electron-positron distribution function with energy $\gamma$ at time $t$.
The full version of the equation is found in \citet{ginzburg1964}, where the radial coordinate and other second order terms are taken into account.
Whereas the lack of radial coordinate $r$ in the equation usually adopted hampers the implementation of a detailed accounting of diffusion of particles and of the magnetic field distribution, it makes the problem more tractable. 
Despite this caveat, this kind of simplified models have been very successful in order to fit the spectra of most of the known PWNe and has taught us as a lot about the underlying PWN physics. Among the most important pieces of information gathered we now know (see e.g., \cite{Aharonian1997,zhang2008, gelfand2009, dejager2009, tanaka2011, martin2012, Bucciantini2011, martin2012, torres2013, vorster2013, Torres2014,torres2018,zhu2021}: 

\begin{itemize}

\item PWNe are particle dominated: magnetization of a few percent are common, although there might be some observational bias embedded into this fact.

\item PWNe are found having very similar injection parameters (high energy slope: 2.2-2.8; low-energy slope: 1.0-1.6)  and break energies (Lorenz factor around $10^{5 \div 6}$), but there is hardly any correlation among parameters.

\item PWNe have high multiplicity (order $10^{4 \div 6}$).

\item Comparing SEDs of the PWNe as observed today mixes pulsars of different spin-down power and age; generates a variety of distributions. A normalized comparison of the SEDs reduces the dispersion leaving no significant outliers. 

\item From the known population of PWNe, only Crab is self-synchrotron Compton (SSC) dominated at high energies (which is understood since we can derive the conditions when this can plausibly happen and prove that essentially only Crab fulfills them). Most others are dominated by the comptonization of the far infrared photon density, with the  
SSC being less relevant. 

\item  The usually adopted GALPROP code may tend to underpredict local values of FIR and NIR needed for PWNe to shine up in TeV as they do. 
GALPROP bins could be coarse for particular PWNe studies.

\item We need to use a more complex dynamics in aged PWNe. Reverberation is very important, but not clearly understood. Reverberation can produce stages of super-efficiency, where, as a result of contraction, e.g., the luminosity in X-rays $L_x > L_{sd}$

\item One key point related to Eq. \ref{eq:difflos} is that, as demonstrated in \citet{martin2012},  all terms of the equation are relevant when studying time evolution, since their relative significance change along time.
Relative differences in flux can reach 100\% between time dependent models once we let them evolve, meaning that the parameters obtained for even a young (worse for a middle-aged PWN) can differ significantly from one model to another depending on the approximations.

\end{itemize}


The recent 3DMHD numerical simulations of the Crab PWN \citep{Porth:2014,Olmi:2016}, despite they are yet limited in time, 
showed that the onset of kink instabilities leads to a large amount of dissipation with respect to that observed at lower dimensionality.
Such dissipation renders the total pressure nearly constant across the whole PWN, as well as to a relatively limited 
magnetic field gradient.
Given that, by definition, one-zone models have zero gradients, this may explain the success of radiative models to explain the spectral energy distribution of PWNe \citep{Gelfand2017}.




In this context, TIDE is a radiative model prepared to compute the spectral and dynamical evolution of young and middle-aged pulsar wind nebulae (PWNe). It solves the diffusion-loss equation for the electron-positron pairs injected into the nebula from the termination shock taking into account synchrotron, inverse Compton (IC), adiabatic, and Bremmstrahlung energy losses, evolving in time the magnetic field, the radius of the PWN
and other dynamical variables, as 
the shock trajectories of the SNR, allowing to follow the evolution further than the beginning of the reverberation phase.
The TIDE code has been used in several occasions to model PWNe observations and study its physics (e.g., see the latest in 
\cite{torres2018,bandiera2020} and references therein. 
Here, we are presenting the latest significant upgrades, including parallelization of some of its components and automatic fitting in a fully time-dependent scenario, and using these upgrades to tackle the issue of possible degeneracy of solutions.
Real fitting of PWNe multi-frequency data is not usually done (despite it is common to find the word {\it fit} in the literature (see the formerly quoted papers, for instance), it rather refer to visually acceptable descriptions when time-dependent, complex models are involved. 
It is thus relevant to briefly clarify what is incorporated in the fitting model and where is room for future improvements.


Using this automatic fitting tool we can more properly address whether degeneracy exists. 
Given a relatively complete set of observational data, how unique are PWN models? 
Can a significantly different set of model parameters lead to similar spectral energy distributions? 
We also use such tool to consider, for the first time to our knowledge, whether a significant set of PWN multi-frequency spectral energy distribution (SED) data can actually imposes a constraint over the age of the system that emits it.
The rest of this work goes directly into presenting the fitting technique and the results obtained, 
leaving a more detailed accounting of the physical components to an Appendix.
















\section{Fitting method}
\label{sec:min}

\subsection{Identifying free parameters}

Time-dependent PWNe models, despite all simplifying assumptions,
are still subject to a significant degree of freedom. 
Table \ref{tab:free-para} identifies all (in principle) free parameters that enter into our model.
The meaning and formulae associated with this set of parameters can be read in Appendix~\ref{sec:appa}.
The parameters in Table \ref{tab:free-para} --unless in obvious cases such as the age of the PWN-- are additionally assumed to be 
constant along time, which is in itself an approximation.
There could be occasions (e.g., spin-down rate transitions, glitches, bursts, etc.) where the
injection or the energy distribution or the magnetization could be subject to sudden changes.




\begin{table*}
\scriptsize
\centering
\begin{threeparttable}
\caption{Identification of the whole set of possible free parameters for a time-dependent automatic fitting. All of these parameters can be changed directly in the code, and searched in the minimization procedure when sufficient data warrants it.}
\label{tab:free-para}
\begin{tabular}{lll}
\hline
 & Name & Comment\\
\hline
\hline
$\delta$ & Data systematic uncertainty & Accounts for the systematic uncertainty between data points obtained from different detectors \\
$t_\text{age}$ & Age of the PWN & Relevant for the particle accumulation, dynamics, energetics. \\
$v_\text{psr}$ & Proper motion of the pulsar & Activates/deactivates the injection of new particles depending on the pulsar position with respect to the PWN shell \\
$n$\tnote{$^\dagger$} & Braking index & Affects the evolution of the spin-down luminosity \\
$L_0$\tnote{$^\dagger$} & Initial spin-down luminosity & Affects the evolution of the spin-down luminosity \\
$\tau_0$\tnote{$^\dagger$} & Initial spin-down age & Affects the evolution of the spin-down luminosity \\
$d$ & Distance & Relevant for the flux normalization \\
$\gamma_\text{min}$ & Minimum energy at injection & Relevant for the spectrum normalization and can affect to the luminosity of the spectrum at radio frequencies \\
$\gamma_\text{b}$ & Energy break & Relevant for the spectrum population and the position of the spectrum break \\
$\alpha_l$ & low-energy index & Relevant for the spectral definition at low energies.  \\
$\alpha_h$ & high-energy index  & Relevant for the spectral definition at high energies. \\
$\varepsilon$ & Containment factor & Relevant for fixing the synchrotron cut-off \\
$E_\text{sn}$ & Supernova explosion energy & Relevant for the evolution of the PWN radius \\
$M_\text{ej}$ & Ejected mass & Relevant for the evolution of the PWN radius \\
$\rho_\text{ism}$ & Interstellar medium density & Relevant for the evolution of the PWN radiu \\
$\delta_\text{core}$ & SNR core density index & Relevant for the evolution of the PWN radius \\
$\omega$ & SNR envelope density index & Relevant for the evolution of the PWN radius \\
$\eta_\text{B}$ & magnetic fraction & Controls synchrotron radiation and magnetic field. \\
$\eta_\text{other}$ & Other energy releases fraction & Controls spectrum normalization and magnetic field. \\
$w_\text{i}$ & Photon background energy density & Affects inverse Compton luminosity. Can come from observations or environmental modelling. \\
$T_\text{i}$ & Photon background energy density & Affects inverse Compton luminosity. Can come from observations or environmental modelling. \\
$n_\text{He}/n_\text{H}$ & He/H density ratio & Relevant for Bremsstrahlung radiation spectrum. \\
\hline
\end{tabular}
\begin{tablenotes}
\item [$\dagger$] Their values are taken consistently 
with the age of the PWN is the latter is let to vary. 
\end{tablenotes}
\end{threeparttable}
\end{table*}

\subsection{Fitting }

Our model now incorporates a new routine package in order to perform fits of multi-wavelength PWNe 
spectra from radio to VHE, in a multi-dimensional space, along time. We refer to this as TIDEfit.
It uses the Downhill Simplex Method (also called Nelder-Mead method, see \cite{numerical_recipes}) where the initial set of parameters is represented in a polygon (called simplex) in a $N$-dimensional parameter space, being $N$ the number of free parameters, with $N+1$ vertices.
The Nelder-Mead method prompts a sort of transformations to this polygon, which provide new solutions that are tested, and accepted or rejected according to its nearness to the data, in such a way that the polygon surrounds the final solution up to the desired accuracy for final convergence.
It is known that such a method suffer of a great loss of efficiency when the number of parameters is high enough \citep{gao2012}.
However, these authors proposed a modified version of the method to improve the efficiency in these situations, which is already implemented in the {\it minimize} function of the {\it Scipy} library for Python.


Getting a good fit via methods such as the Monte Carlo Markov Chain implies realizing a high number of iterations  (of the order of tens of thousands or even more) in which the diffusion-loss equation must be solved at every time. 
Due to the complexity that implies solving this equation in detail and the potentially high number of free parameters, other typical methods (e.g., minimum least squared) which require the calculations of first partial derivatives (and probably second partial derivatives too) with respect of all the free parameters are not preferred either, due to their time-consumption.
Instead, the simplex method avoid the use of first partial derivatives and can typically converge in $\sim 1000$ iterations for 6-8 free parameters if the initial choice of the simplex is not extremely far from the fitting one.
See Table 1 for typical search ranges.

\subsection{Data treatment}


The multi-wavelength data for a generic PWN come from different detectors, and from different observation times (sometimes decades apart). Thus, in addition of the observational uncertainty of each detector per se and the complexities arising from treating data with vastly different relative error bars, there exist a cross-calibration uncertainty that must be taken into account.
This issue --as well as some misconceptions arising when fitting data in astrophysics-- are discussed by \cite{hogg2010}.
We pose that the most convenient prescription is to introduce an extra free parameter in our fit such that the real uncertainty for each point is systematically enlarged in proportion to the flux calculated from the model. 
That is (see \cite{hogg2010}),
\begin{equation}
s_i^2 = \sigma_i^2 + \delta^2 F_i^2,
\end{equation}
where $\sigma_i$ is the observational uncertainty, $\delta$ is the free parameter that accounts for the systematic uncertainty and $F_i$ is the spectral flux obtained from the model.
As discussed in \cite{hogg2010}, the $\chi^2$ function typically used in astrophysics is only a part of the full formula for the likelihood $\chi^2$ function, because the coefficient depending on the uncertainty is neglected. This approximation commonly done in the literature is right if we do not consider systematic errors in our data, so they are constant independently of the values of the model parameters. Here, instead of minimizing the typical $\chi^2$ function, we maximize the likelihood function \citep{hogg2010}
\begin{eqnarray}
\label{eq:likefunc}
\ln p(y | \nu, \sigma, \vec{x}, \delta) = - \frac{1}{2} 
\sum_i \left \{ \frac{[y_i(\nu) - F(\nu,\vec{x})]^2}{s_i^2(\nu,\sigma_i,\delta)} - \right.  \nonumber \\  \left. 
\ln 2 \pi s_i^2(\nu,\sigma_i,\delta) \right \},
\end{eqnarray}
where $\nu$ is the observed frequency and $\vec{x}$ is our set of free parameters.
Note that Eq.~\ref{eq:likefunc} uses the logarithm of $p$
to have smoother control over the function's gradient and the variation tolerance before consider convergence of the numerical scheme.

If the likelihood function has a complicated shape in the parameter phase-space, it is possible --in the context of the Nelder-Mead method-- to fall into local minima of the likelihood function 
depending on the first guess of the parameters. 
To avoid this, we noticed that a good practise is to fix first the systematic uncertainty at an a-priori seemingly reasonable value (10 - 20\%), 
then fit the rest of parameters, 
in order to use the results as an initial guess for a second refined fit now leaving the systematic uncertainty as a free parameter too.
All of this is automatically implemented in our code (see Appendix B for some further notes on implementation and numerical cost).




\section{Applications to two canonical PWNe with varying data quality: How
unique are their fits?}
\label{sec:results}

We now fit the spectra of two canonical PWNe with TIDE: the Crab Nebula and 3C 58. Fortunately, most of the parameters can be fixed for observations in both nebulae, as we can see in Table \ref{tab:fixed}. However, the number of data points and their uncertainty varies.
For instance, whereas the age of Crab is certain, the age of 3C 58 is assumed. 
A total of 8 parameters has been left free to fit both PWNe.
Apart the systematic uncertainty, these are the energy break and the spectral slopes of the injection function, the ejected mass of the SNR, the magnetic fraction and the FIR and NIR energy densities.

A special note goes for the containment factor. 
The value taken for Crab comes from a previous fit where this parameter is free, but since it only affects the synchrotron cut-off where the number of data points is small, the uncertainty in this parameter is high because the $\chi^2$ changes little when we look for the confidence level. 
Thus we have taken it fixed to minimize the number of parameters.
We come back to this in Appendix~\ref{sec:appa}, and show in detail how different containment factors affects the spectra.
The fitting ranges and the results obtained for each nebula are compiled in Table \ref{tab:results}.

\begin{table}
\scriptsize
\centering
\caption{Data from observations -first panel-, computed from the first panel -second panel-, or assumed -third panel- in the models of the Crab nebula and 3C 58.}
\label{tab:fixed}
\begin{tabular}{llcc}
\hline
\hline
Period (ms) & $P$ & 33.4 & 65.7\\
Period derivative (s s$^{-1}$) & $\dot{P}$ & $4.2 \times 10^{-13}$ & $1.93 \times 10^{-13}$\\
Braking index & $n$ & 2.509 & 3\\
Observed radius (pc) & $R_\text{pwn}$ & 1.8 & 2.3\\
Distance (kpc) & $d$ & 2 & 2\\
Pulsed radiation fraction & $\eta_\text{other}$ & 0.0142 & 0.0018\\
\hline
Characteristic age (yr) & $\tau_\text{c}$ & 1296 & 5398\\
Age (yr) & $t_\text{age}$ & 968 & 2500\\
Spin-down  (erg s$^{-1}$) & $L$ & $4.5 \times 10^{38}$ & $2.7 \times 10^{37}$\\
Initial spin-down  (erg s$^{-1}$) & $L_0$ & $3.1 \times 10^{39}$ & $9.3 \times 10^{37}$\\
Initial spin-down age (yr) & $\tau_0$ & 750 & 2878\\
\hline
SN energy explosion (erg) & $E_\text{sn}$ & $10^{51}$ & $10^{51}$\\
ISM density (cm$^{-3})$ & $\rho_\text{ism}$ & 0.5 & 0.1\\
FIR temperature (K) & $T_\text{cmb}$ & 70 & 25\\
NIR temperature (K) & $T_\text{cmb}$ & 5000 & 2900\\
Containment factor & $\varepsilon$ & 0.3 & 0.5\\
\hline
\end{tabular}
\end{table}

\subsection{Crab Nebula (with LHASSO data)}
\label{sec:crab}


\begin{figure*}
\centering
\includegraphics[width=0.45\textwidth]{crab_sed-eps-converted-to.pdf}
\includegraphics[width=0.45\textwidth]{crab_freeage_sed-eps-converted-to.pdf}
\\
\includegraphics[width=0.45\textwidth]{3c58_sed-eps-converted-to.pdf}
\includegraphics[width=0.45\textwidth]{3c58_freeage_sed-eps-converted-to.pdf}
\caption{Models fitted for the Crab Nebula (top) and 3C 58 (bottom). On the left column, the age is fixed and on the right, it is a free parameter (see results in Table \ref{tab:results}) The data points for Crab are taken from \cite{baldwin1971,maciasperez2010}, (compilations of radio observations), \cite{buehler2014} (from IR to $\gamma$-rays and VHE) and the new points offered by the LHAASO collaboration \citep{lhaaso2021}. For 3C 58, they come from \cite{green1986,morsi1987,salter1989} (radio), \cite{torii2000} (infrared), \cite{green1994,slane2008} (X-rays), \cite{abdo2013,ackermann2013,li2018} (gamma-rays) and \cite{aleksic2014} (VHE).}
\label{fig:sed}
\end{figure*}



The Crab Nebula is the most studied PWN, and has precise measurements of its flux at all frequencies from radio to VHE.
Regarding the flux data, we included the compilation shown in \cite{buehler2014} and the recently VHE data points by LHAASO \citep{lhaaso2021}.
The resulting spectral energy distribution (SED) from the automatic fitting process is shown in Figure \ref{fig:sed}.
The automatic fit has reduced $\chi^2$ of 1.013, and the 1$\sigma$ band (depicted in grey in the plot) is almost invisible at all energies in this scale (compare with the companion panel for 3C 58, where this band is more prominently seen).
The fit reproduces the data, particularly in critical regions such as in the synchrotron decay or in the gamma-ray peak.
The Crab PWN is SSC dominated at the gamma-ray peak, with a significant contribution of comptonized CMB photons at larger energies.
Despite the generous range of the parameters search, there are no surprises in the fit, neither in the radiative process dominance nor in the free parameters values.
There is no alternative set of solutions outside the range quoted for the free parameters providing a similarly good fit. 
We have tested this thoroughly, starting from dozens of very different initial conditions in the Nelder-Mead simplex:
We performed $\sim 150$ fits with randomly selected initial simplex 
obtaining convergence to the solution quoted in more than 73\% of the cases. 
In the remaining ones, visual inspection or the value of the reduced $\chi^2$ would immediately qualify the fit as non-compliant.
These non-convergent cases happen when the values of the energy break $\gamma_b$ and/or the magnetic fraction $\eta_\text{B}$, and/or the initial systematic uncertainty
are close to one of the extreme values inside its range of variation.
The Crab nebula spectra is thus uniquely fit in the surrounding of the parameters described, and in the surrounding of no others. 
There seems to be little room for degeneracy in our PWNe models when such significant set of data is considered.



It is interesting to remark that the best model actually slightly overpredicts the flux at TeV energies, especially around 100 TeV (see Figure \ref{fig:sedzoom} to zoom in the details). 
\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{crab_sed_zoom-eps-converted-to.pdf}
\caption{Zoom of the Crab Nebula spectrum shown in Figure \ref{fig:sed} between 100 GeV and 1 PeV.}
\label{fig:sedzoom}
\end{figure}

Given that the PWN radius in the model is just an average value representing a more complex physical situation, one can think that the inverse Compton yield at high energies could be reduced by if the radius of the nebula is larger, at least within the corresponding uncertainty.
If so, the magnetic energy density would be lower and, in consequence the magnetic field.
In turn, the magnetic fraction could be larger in order to power up the magnetic field sufficiently, decreasing the normalization of the particle distribution function to compensate.
Despite of the reduction of the normalization, the synchrotron luminosity could be fitted because it is proportional to $\sim B^2$.
A lower number of pairs would decrease the IC flux helping to fit the IC flux in both components (IC-CMB and SSC).
This strategy may however be considered disadvantageous since it would require an ejected mass of about $6M_{\odot}$, which seems too low for current models of supernovae \cite{Smartt:2009}. 
\cite{Owen2015} state
that the total mass of gas plus dust in the Crab Nebula is $7.2 \pm 0.5 M_\odot$, consistent with a progenitor star mass of $\sim
9 M_\odot$.
However, 
we note on the side that disregarding this fact, and allowing for smaller ejecta masses than the range we adopt in Table 
\ref{tab:fixed}, we obtain even better fits than the ones shown.


Another option to decrease the inverse Compton yield at high energies is that part of the energy invested in the acceleration of electron-positron pairs, $(1-\eta_\text{B})$, is used to accelerate an hadronic component, or to emit pulsed emission, or other multi-messenger emission such as neutrinos or gravitational waves.
That is, that part of the energy is simply not available due to it being spent in other processes that are or may be, in fact, naturally occuring.
We have taken this into account. We take the fraction of pulsed emission integrating the data from \cite{Torres2019}, which corresponds to 1.42\% of the spin-down luminosity, and we have actually discounted such a quantity from our contribution to particles (as explained in the Appendix~\ref{sec:appa}). However,  Figure \ref{fig:sed} shows that this imply little change.

Figure \ref{fig:crabchi} shows the shape of the reduced chi-squared function 
around the solution given by TIDE for each parameter pair. 
Most of the uncertainty in the model is introduced by the subdominant FIR and NIR energy densities: the 1, 2 and 3$\sigma$ ellipsoids are elongated and cover the entire possible value parameter range. This occurs as a result of the fact that none of them is actually relevant for the fit. Interestingly, however, the automatically selected FIR inverse Compton contribution is larger than the NIR one. 
The rest of the parameters, those related with the injection function (energy break and low and high energy injection slopes) together with the magnetic fraction and the ejected mass are well determined in the $\chi^2_r$ maps, confirming that the degeneracy is minimal.





\subsubsection{The automatic fit in the context of other models}

The automatically chosen parameters are compatible with earlier studies using previous incarnations of the code (e.g., \cite{Torres2014, aleksic2015}). 
Differences introduced by the improved treatment of the dynamics explain 10-20\% changes in the resulting radius and ejected mass (we now discount magnetospheric pulsed emission, use the new SNR/PWN shock positions as in \cite{bandiera2021}, time-dependence of the spin-down luminosity in the radius evolution equations as in \cite{martin2016}, and an improved treatment for the energy conservation in the dynamical evolution as in \cite{bandiera2020}).
The small variation in radius also reflects in the magnetic fraction, while the magnetic field is similar.
Regarding the FIR and NIR energy densities, the fitting algorithm confirms that the IC-FIR and IC-NIR contributions are negligible in order to explain the IC flux.

Apart from TIDE, other models use the new data from LHAASO to update their corresponding fit, but it is difficult to compare them directly due to the features of each model being different.
Here, we will compare our results with three different models. One of them is previous to the LHAASO data, and based on the model of \cite{meyer2010} and published with the MAGIC data points released in \cite{aleksic2015}.
The other two models include the new data: the paper of the LHAASO collaboration itself \citep{lhaaso2021} and the one proposed in \cite{nie2022}.
\cite{meyer2010} (hereinafter MEY) fits the Crab spectrum with a steady one zone model composed by two pair populations: one which accounts for the radio emitting pairs and another one for the wind accelerated particles. These two populations are not calculated by evolving the diffusion-loss equation but  as a fit of the particle spectrum observed today. 
The model uses 15 parameters to fit the spectrum.
These are thus significant differences, which makes a direct comparison unfeasible.
The model used in \citep{lhaaso2021} (hereinafter LHA) is also a steady one zone model, but more simplified in comparison with MEY, needing only three parameters, in particular,  the spectral index at high energies, the magnetic field, and the super exponential energy cut-off.
Finally, TIDE and the model of \cite{nie2022} (hereinafter NIE) assume a broken power-law injection of particles and evolve them in time taking into account adiabatic and only synchrotron radiative energy losses.
From the right-hand side of Equation \ref{eq:difflos}, NIE neglects the second term and radiative losses are included in the escape term (third term of Equation \ref{eq:difflos}).
As it was shown in \citealt{martin2012}, this kind of approximations has a significant impact in the results along the time evolution, obtaining differences in the calculated flux of more than a factor 2 by using the same parameters.
Qualitatively bearing in mind account all said differences, selected model values are not far from one another.
For instance, the maximum energy at injection for the models MEY, LHA and NIE is, respectively, 3.7, 2.15 (we take the super-exponential energy cut-off as a reference) and 5.649 PeV (this is the initial maximum energy at injection, so we should take it as an upper limit).
In the case of TIDE, the maximum energy gives 4.7 PeV (see Table \ref{tab:results}) which is an intermediate value of those obtained by these models.
We have seen in \cite{lhaaso2021,nie2022} some trials to include an hadronic component into the spectrum to explain the enhance of the VHE flux at $~$1 PeV in the Crab Nebula.
Introducing this component is certainly possible, but for the current set of data, there is not enough data to really favour a lepto-hadronic model over a leptonic one: it would introduce more free parameters (at least two) to fit only one additional data point.
However, the appearance of such enhancement may require this inclusion in the future.

\subsubsection{A lower-energy explosion}
\label{crab-low}

It has been earlier suggested that SN 1054 was an atypical, low-energy explosion (e.g., \cite{Frail1995,Yang2015}).
Here, we shall also explore this possibility by considering, on the one hand, a model with $E_{sn} = 10^{50}$ erg, and on the other
by additionally consider $E_{sn}$ as a free parameter within a range $10^{49}$ --  $10^{51}$ erg. 

For the $E_\text{sn}$ case, the model clearly worsens and requires an ejected mass lower than 7M$_\odot$ to be competitive. 
If we prevent this, then we cannot fit the radius of the nebula correctly. 
For the latter case, when we allow $E_\text{sn}$ to vary, we obtain a very similar result as show in Table~\ref{tab:results}: $E_\text{sn}=8.3^{+0.5}_{-0.4} \times 10^{50}$ erg, which is close to the a priori canonical value of $10^{51}$ erg as initially assumed.







\subsection{3C 58}
\label{sec:3c58}

3C 58 is a PWN associated with the pulsar PSR J0205+6449 \citep{murray2002}.
Regarding the morphology, this PWN shows a faint jet and a toroidal structure perpendicular to it.
3C 58 has been detected in radio, IR, X-rays, and gamma-rays (see the references in the caption of Figure \ref{fig:sed}).
This makes 3C 58 another good example to test the automatic fitting code, albeit
there is still an ongoing discussion on its age and distance in comparison to the more certain values of the Crab nebula. 
Regarding the age, \cite{stephenson2002} associate the nebula to the supernova SN 1181, resulting a current age of 841 yr, but other estimates give an age of up to 7 kyr (see e.g., \citealt{fesen2008} for a compilation).
For our model, we use the dynamical age of 2.5 kyr given by \cite{chevalier2004,chevalier2005} and come back to this discussion below, with a new approach.
We also note that the required ejected mass in order to fit the radius using an age of 841 yr would also be low ($< 8M_\odot$).
On the other hand, the distance ranges from 2 kpc \citep{kothes2010} to 3.2 kpc \citep{roberts1993}.
We take the newest determination of 2 kpc, where new H I data is used.

We obtain a fit with a reduced chi-squared of 0.74 and a systematic uncertainty of 0.01, which is on the limit considered for this parameter.
In terms of the reduced $\chi^2$, the fit is a bit farther from ideal, but we understand that this may be strongly influenced by the number of data points considered in X-rays (189 of 232 in total) and radio to a lesser extent (24 of the remaining 43 taking out the X-ray points), where the fit is really good. The data outside these ranges have errors larger than an order of magnitude, then their contribution to the reduced $\chi^2$ becomes really low.
Regarding the solution convergence (see Figure \ref{fig:3c58chi}), if we compare with the solution for the Crab Nebula, the ellipsoids in the $\chi^2$ maps seem more elongated. Therefore the relative errors of the parameters are higher, but in terms of degeneracy, the ellipsoids are still closed and the solution is not degenerate.

As in the case of the Crab Nebula,  we also discount from the energy input for particles the pulsed radiation emitted ($1.8 \times 10^{-3}$, see Table~\ref{tab:fixed} and \citealt{li2018}), with no noticeable effect on the spectrum and the parameters. The modelled flux surpasses the data points at energies higher than 100 GeV. For Crab, this situation gives a chance to speculate on possible corrections to the model due to the high quality of the data, but in this case, we have less data points and their dispersion makes this difference no-significant in comparison with the variations of the model in other parts of the IC spectrum.  

Finally, regarding convergence of the solution, and as in the previous case, we have also performed $\sim$150 simulations with random sets of initial parameters converging to the final solution in more than 70\% of the cases. The initial sets of parameters for the non-convergent cases have in general values close to the limit range for the energy break $\gamma_\text{b}$ and magnetic fraction $\eta_\text{B}$. Thus, a very similar scenario and results as in the case of the Crab Nebula appears.




\subsubsection{The fit in the context of other models}

There are not many models of 3C 58 apart from \cite{tanaka2013} and 
the studies we did earlier with TIDE in the framework of different observations
\cite{torres2013,li2018}.
As expected from the outcome of Crab, the differences with the previous fit published in \cite{li2018} are also few, and are mainly explained by the improvement on the expansion model and the mathematical fitting used in the version presented in this paper, as commented before for the case of Crab.


Regarding the model of \cite{tanaka2013}, we recall that the nebula was detected one year later at gamma-rays and VHE \cite{abdo2013,ackermann2013,aleksic2014} and thus these data has not been used then. Despite of this, there is a general good agreement with our obtained slopes for the low and high energy injection power-laws (1.0 \& 3.0, respectively) and the location of the energy break in the two alternative models that they present ($5-9 \times 10^4$).
The expansion model is a ballistic expansion of the nebula where the parameter to be fitted is the expansion velocity.
In this case, the two models are indeed different, from 780 to 2000 km s$^{-1}$.
In our model, the expansion velocity at 2.5 kyr is 1056 km s$^{-1}$, which is roughly in agreement with mentioned models too, but in our case we have to take into account that it is not constant and changes with time like $\sim t^{1/5}$. 
Finally, the magnetic field obtained by us is 15.6 $\mu$G, which is close but lower than the minimum value given by \cite{tanaka2013} where the magnetic field lies between 17 and 40 $\mu$G. This small deviation in the magnetic field can be produced by the disagreement of the models on the considered radius of the PWN. While \cite{tanaka2013} consider a radius of 2 pc for 3C 58, our model predicts a slightly larger radius of 2.31, which is closer to what is observed. The magnetic fraction in \cite{tanaka2013} is also a factor 3 lower than what we obtain from the fit. 
However it is very difficult to compare both models directly since the evolution in time is computed differently.
In fact, we remark that comparing models based on different theoretical priors provides not more than a qualitative feeling of the stability of the solution with respect to the assumptions.



\begin{table*}
\begin{threeparttable}
\center
\scriptsize
\caption{Parameters obtained for the Crab Nebula and 3C 58 from TIDEfit. The parameters where a searching range is defined are those used to do the fit. The quantities inside the parenthesis are the confidence interval at 95\%.}
\label{tab:results}
\begin{tabular}{llccccc}
\hline
Parameter & Symbol & Range & \multicolumn{2}{c}{Crab Nebula} & \multicolumn{2}{c}{3C 58}\\
 & & & Age fixed & Age free & Age fixed & Age free\\
\hline
\hline
Free parameters \\
\hline
Age (yr) & $t_\text{age}$ & [300,1500 / 2000,5000] & 968 & 914 (901,931) & 2500 & 2108 (2075,2140)\\
Magnetic fraction ($\times 10^{-2}$) & $\eta_\text{B}$ & [0.01,50] & 2.048 (2.027,2.053) & 2.2 (2.1,2.4) & 1.06 (1.03,1.09) & 1.77 (1.70,1.84)\\
Ejected mass ($M_\odot$) & $M_\text{ej}$ & [7,12 / 12,25] & 7.9 (7.6,8.1) & 7.0 (7,7.2) & 17.2 (16.8,17.6) & 10.8 (10.5,11.0)\\
Energy break ($\times 10^5$) & $\gamma_{b}$ & [0.1,100] & 4.9 (4.7,5.3) & 5.8 (5.4,6.2) & 0.88 (0.86,0.90) & 0.76 (0.74,0.78)\\
Low energy index & $\alpha_{l}$ & [1,4] & 1.46 (1.44,1.48) & 1.51 (1.49,1.53) & 1.00 (1,1.03) & 1.00 (1,1.04)\\
High energy index & $\alpha_{h}$ & [1,4] & 2.470 (2.468,2.473) & 2.478 (2.476,2.481) & 3.012 (3.009,3.015) & 3.012 (3.009,3.015)\\
FIR energy density (eV cm$^{-3}$) & $w_\text{fir}$ & [0.1,5] & 0.10 (0.1,0.15) & 0.10 (0.1,0.18) & 0.22 (0.1,0.42) & 0.37 (0.1,0.66)\\
NIR energy density (eV cm$^{-3}$) & $w_\text{nir}$ & [0.1,5] & 0.1 (0.1,1.38) & 0.3 (0.1,2.4) & 0.41 (0.1,0.99) & 0.75 (0.1,1.56)\\
Systematic uncertainty & f & [0.01,0.5] & 0.106 & 0.103 & 0.010 & 0.010\\
\hline
Derived values\\
\hline
Magnetic field ($\mu$G) & $B$ & \ldots &  86.7 (84.4,89.1) & 88.5 (86.6,91.4) & 15.6 (15.3,15.9) & 18.0 (17.6,18.4)\\
Maximum energy at injection ($\times 10^9$) & $\gamma_{max}(t)$ & \ldots & 9.25 (9.21,9.26) & 9.6 (9.4,10.1) & 2.70 (2.66,2.74) & 3.51 (3.44,3.58)\\
PWN Radius (pc) & $R_\text{pwn}$ & \ldots & 1.82 (1.79,1.85) & 1.78 (1.73,1.83) & 2.31 (2.28,2.33) & 2.30 (2.25,2.35)\\
\hline 
Statistics\\
\hline
Degrees of freedom\tnote{$^\dagger$} & DOF & \ldots & 580 & 579 & 456 & 455\\
Reduced chi-squared & $\chi^2_\text{r}$ & \ldots & 1.013 & 0.987 & 0.740 & 0.718\\
\hline
\end{tabular}
\begin{tablenotes}
\item [$\dagger$] Age free models have one more free parameter, but also one more data point, the spin-down luminosity. It is necessary to fix correctly the initial spin-down luminosity and the initial spin-down age.
\end{tablenotes}
\end{threeparttable}
\end{table*}

\begin{figure*}
\raggedleft
\includegraphics[width=0.16\textwidth]{crab_gridmap_8_9.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_8_10.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_8_13.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_8_17.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_8_22.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_8_24.jpg}
\\
\includegraphics[width=0.16\textwidth]{crab_gridmap_9_10.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_9_13.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_9_17.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_9_22.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_9_24.jpg}
\\
\includegraphics[width=0.16\textwidth]{crab_gridmap_10_13.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_10_17.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_10_22.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_10_24.jpg}
\\
\includegraphics[width=0.16\textwidth]{crab_gridmap_13_17.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_13_22.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_13_24.jpg}
\\
\includegraphics[width=0.16\textwidth]{crab_gridmap_17_22.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_17_24.jpg}
\\
\includegraphics[width=0.16\textwidth]{crab_gridmap_22_24.jpg}
\caption{Reduced $\chi^2$ maps for Crab Nebula. The dark blue, turquoise and yellow contours correspond to the 1, 2 and 3$\sigma$ confidence levels.}
\label{fig:crabchi}
\end{figure*}

\begin{figure*}
\raggedleft
\includegraphics[width=0.16\textwidth]{3c58_gridmap_8_9.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_8_10.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_8_13.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_8_17.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_8_21.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_8_23.jpg}
\\
\includegraphics[width=0.16\textwidth]{3c58_gridmap_9_10.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_9_13.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_9_17.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_9_21.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_9_23.jpg}
\\
\includegraphics[width=0.16\textwidth]{3c58_gridmap_10_13.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_10_17.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_10_21.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_10_23.jpg}
\\
\includegraphics[width=0.16\textwidth]{3c58_gridmap_13_17.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_13_21.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_13_23.jpg}
\\
\includegraphics[width=0.16\textwidth]{3c58_gridmap_17_21.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_17_23.jpg}
\\
\includegraphics[width=0.16\textwidth]{3c58_gridmap_21_23.jpg}
\caption{Reduced $\chi^2$ maps for 3C 58. The contours have the same correspondence as in Figure \ref{fig:crabchi}.}
\label{fig:3c58chi}
\end{figure*}

\section{Using detailed observations and PWN models to determine the system's age}

The quantity of data available and its completeness in the entire electromagnetic spectrum in both cases allow us to explore using 
TIDEfit to estimate the system's age directly from the fit, adding it to the rest of the free parameters.
In Table \ref{tab:results}, we show the results obtained for both nebulae when such approach is adopted.
In terms of reduced $\chi^2$, the quality of the fits are quite similar when comparing with those where the age is fixed, and the age obtained is very close to the expected.
In the case of Crab, the true age falls inside the 2$\sigma$ range obtained from the model.
For 3C 58, for which a real age is not known, we find an age of 2059 yr, which implies a $\sim$440 yr younger PWN than assumed in the previous model.
The change in the age has concurrently changed the values of the magnetic fraction, the ejected mass, and the FIR and NIR energy densities.
The rest of parameters remain vary similar, included derived parameters as the magnetic field, the radius and the maximum energy at injection.
A difference of only 440 years is in fact already pointing towards the larger values in the age dichotomoy stated earlier (where the 
discussion was between 871 and 2500 years).
This is the first time that a radiative code is used to fit the age of PWNe.
Of course, it is crucial that a complete set of data, from radio to VHE exists,  to avoid 
possible degeneracies when the age is fitted.


\section{Concluding remarks}
\label{sec:conc}

\begin{itemize}


\item We have considered how to take into account that the multi-wavelength data for a generic PWN come from different
detectors, and from different observation times (sometimes separated decades
apart). Using the prescriptions given in \citet{hogg2010}, we proposed that a cross-calibration
uncertainty must be taken into account when fitting models, and posed that introducing an extra parameter to deal with it is the best option. 

\item We have successfully incorporated a fitting technique into a fully-time dependent numerical model for PWN. 
The appendix~\ref{sec:appa} accounts for the physical components of the model.
They are treated in much the same way as in earlier incarnations of TIDE, prior to automatic fitting,  
including some improvements like those related to the dynamics, or the accounting of the energy reservoir spent in non-PWN related emission, and others mentioned. 
The difference with former versions of the code, and in general with other codes, is 
that all of these components are now included into an automatic fitting technique. This is a result of an improved single-model computational cost (see Appendix~\ref{sec:imple})
as well as the introduction of a relatively fast algorithm technique that avoids the use of derivatives. 
The algorithm introduced in the code TIDEfit is based on the Nelder-Mead search for solutions, which was also used recently in 
the search of best models for magnetospheric emission from pulsars (see \cite{Iniguez2022}).
TIDEfit is able to cope with variations and minimization of a significant number of parameters (see Table \ref{tab:free-para}).

\item Our model maximizes the likelihood  $\chi^2$ function to search for a solution  instead of minimizing the typical expression for $\chi^2$ function. This is done to take into account systematic uncertainties that appear due to the diverse origin of the data, which is generally ignored in the literature regarding radiative models for PWNe. Our model fits the multi-frequency data once the free parameters of the model
are selected. This is done in a modular way, such that a different number of parameters (and which are those) can be chosen according to the problem at hand
and the number and quality of the data points at our disposal. 

\item This represents an improvement over other PWNe models, where the `fits' are in fact qualitatively obtained without a quality assessment. 
In addition, it allows us to tackle whether there is any significant degeneracy where very different sets of parameters would lead to equally good fits.
We performed convergence tests with random sets of initial conditions for the cases studied in this paper (Crab Nebula and 3C 58). The model converges to the final solution in more than 70\% of the cases for both nebulae independently of the initial set of parameters.
When not, parameters selected are on the verge of the search ranges and the fits are non-compliant.
There is no significant degeneracy since we find that neither Crab nor 3C58 data can be encompasses by any other set of parameters leading to an equally good solution.

\item Another application we introduced was the determination of the system's age by assuming it as a free parameter of the model. 
We find that when a significant dataset is considered, the age can be fixed from the spectrum at a relatively high accuracy, for instance, $<50$ years in the case of Crab.
We have used this possibility to actually provide a new estimation of the age of 3C58, for which a discussion is yet ongoing. We obtain an estimation 
of its age at $2108_{-33}^{+32}$ years, within 95\% confidence level. 


\end{itemize}



\section*{Data Availability}

No new observational data is herein presented. The dataset used in all figures is publicly available, and references are given correspondingly.
Any additional theoretical detail required is available from the authors on reasonable request. 


\section*{Acknowledgements}

We thank 
Elena Amato,
Rino Bandiera,
Niccolo Bucciantini,
Agnibha de Sarkar, 
Barbara Olmi,
and
Wei Zhang
for comments, discussions, encouragement, and tests. 
This work has been supported by the grants INAF grants MAINSTREAM 2018, SKA-CTA, PRIN-INAF 2019, by ASI-INAF n.2017-14-H.O, and the Spanish grant
PID2021-124581OB-I00 and  
by the program Unidad de Excelencia María de
Maeztu CEX2020-001058-M. 









\bibliographystyle{elsarticle-harv} 
  