\section{Introduction}

\subsection{Scrambling} \label{subsect_scr}

Every physical system carries information. Due to the time evolution of the system, the state of the information within it also evolves. Consider a quantum system of $K$ degrees of freedom, prepared in a pure state and evolved unitarily in time. Even though the state of the system remains pure throughout its unitary evolution, nevertheless, the system thermalizes after some time. This thermalization takes place in the following sense: although the time evolution is unitary, over time, the initial state becomes entangled with the other states in the Hilbert space of the system. The information initially contained only in the original state is distributed over the other states and is therefore mixed up. Thus, the system ``scrambles''~\cite{Hayden_2007,Sekino_2008} information. The characteristic time scale associated with this process is called the \emph{scrambling time} $t_s$~\cite{Hayden_2007,Sekino_2008}.

Before characterizing the scrambling time, we need to first address scrambling itself. Following~\cite{Hayden_2007,Sekino_2008,Susskind_2011}, which introduced the concept of scrambling, we call a system \emph{scrambled} if the information within it (i.e.\ the state of the system) is sufficiently distributed over the entire Hilbert space of the system with respect to some chosen measure.

Note that we keep the definition general on purpose: we do not limit the formulation to a specific measure of the uniformity of the state distribution. There exist multiple distinct such measures defined in the literature. Some of the first explicit definitions include the Haar measure in application to the choice of the mixing unitary transformation~\cite{Hayden_2007}, or to states (``Haar-scrambled'')~\cite{Sekino_2008}. Following~\cite{Page_1993}, for a model of $K$ qubits in~\cite{Sekino_2008}, the authors also define a system to be ``Page-scrambled'' when the entanglement entropy of any subsystem of $n<K/2$ qubits satisfies $S_n = n - O(\exp[2n-K])$. In general, a detailed picture of scrambling, especially for black holes, requires an understanding of scrambling at a microscopic level. In particular, the microscopic picture of scrambling, introduced in~\cite{Dvali_2013_1,Dvali_2014,
Dvali_2015_1,Dvali_2013_2,Dvali_2016_2} (see also further references therein), establishes the concept of ``memory modes'' (see the discussion below). Moreover, it discusses the concepts and the corresponding time-scales of ``one-particle entanglement''~\cite{Dvali_2013_2} and ``maximal entanglement''~\cite{Dvali_2016_2} for microscopic models of black holes.

However, in the present work we introduce a novel regime, which takes place on time-scales smaller than the various aforementioned scrambling time-scales. We thus establish a new pre-scrambling stage of a system's evolution, which requires quantification. In addition, we note that the results of the present paper are independent of the choice of a specific measure for scrambling.

What is the shortest possible scrambling time that a system can have? Following~\cite{Hayden_2007,Sekino_2008,Susskind_2011}, systems that scramble information in a time logarithmic in the number of degrees of freedom are the quickest scramblers. As in~\cite{Hayden_2007,Sekino_2008,Susskind_2011}, we refer to such systems as \emph{fast scramblers}. We broaden the above definition to also include the number of sites, modes, or generic quantum labels in a system, in addition to the number of degrees of freedom. Furthermore, to keep the definition general, we do not include any temperature dependence or any bound due to the latter, as is, for example, done in~\cite{Sekino_2008,Shenker_2014,Kitaev_2015,Maldacena_2016} among many other works. The results of the present paper apply in general, and require only that fast scramblers have a scrambling time that is logarithmic in the number of degrees of freedom $K$. The temperature dependent scrambling time for the model considered in section~\ref{sect_model} will be analyzed elsewhere~\cite{scr_to_appear}. In this paper we investigate the behavior of the above model at times preceding its scrambling time.

\subsection{Enhanced memory capacity}

Before considering the model of section~\ref{sect_model} we first need to introduce the concepts that are at its foundation. A sub-direction of the larger program of~\cite{Dvali_2018_1,Dvali_2019_1,Dvali_2013_1,Dvali_2014,
Dvali_2015_1,Dvali_2015_2,Dvali_2016_1,Dvali_2017,
Dvali_2018_2,Dvali_2018_3,Dvali_2019_2,Dvali_2019_3,
Dvali_2019_4,Dvali_2020_1,Dvali_2020_2} (see also further references therein) investigates general phenomena of systems that possess states with a high capacity to store information. As in~\cite{Dvali_2018_1,Dvali_2019_1,Dvali_2014,
Dvali_2015_1,Dvali_2015_2,Dvali_2016_1,Dvali_2017,
Dvali_2018_2,Dvali_2018_3,Dvali_2019_2,Dvali_2019_3,
Dvali_2019_4,Dvali_2020_1,Dvali_2020_2}, we refer to such systems as those with \emph{enhanced memory capacity}. One of the above phenomena is the effect of ``memory burden''~\cite{Dvali_2018_1,Dvali_2019_1}, which we briefly summarize below along with other terminology. For this, we largely rely on the recent detailed review included in~\cite{Dvali_2020_2}.

The degrees of freedom of a physical system are commonly described as quantum oscillators. The states of the system are then labeled by distinct sequences of the $K$ oscillators' occupation numbers $\ket{n_1, \dotsc, n_K}$. We refer to each such sequence as a \emph{memory pattern} which stores quantum information. The number of distinct patterns that can be stored within a microscopically narrow energy gap are a measure for the memory capacity of a system~\cite{Dvali_2018_2,Dvali_2019_2}. If the $\mathcal{N}$-many microstates $\ket{n_1, \dotsc, n_K}$ describing distinct patterns are degenerate in energy, they contribute to the microstate entropy $S = \ln (\mathcal{N})$. If a memory pattern contains a large amount of quantum information, it stabilizes the system in the state of enhanced memory capacity and slows down the system's evolution. Following~\cite{Dvali_2018_1,Dvali_2019_1}, we refer to this as the \emph{memory burden} effect.

A system can reach a state of enhanced memory capacity by the effect of \emph{assisted gaplessness}~\cite{Dvali_2019_2}. In its essence, assisted gaplessness occurs when a highly occupied \emph{master mode} interacts attractively with a set of \emph{memory modes}, lowering their energy gaps. The memory modes become effectively gapless and degenerate in energy, and can store large amounts of information. In their own turn, the memory modes backreact on the master mode via the memory burden effect~\cite{Dvali_2018_1,Dvali_2019_1} and slow down the change in its occupation number. As proposed by~\cite{Dvali_2018_1} and investigated in~\cite{Dvali_2020_2}, the memory burden effect can be avoided if the system possesses another set of modes $K'$ to which it can rewrite the information stored in the modes $K$.

\section{A prototype model} \label{sect_model}

We consider a specific prototype model of a system with enhanced memory capacity given in Eq.\ (34) of~\cite{Dvali_2020_2}, which is constructed in correspondence to the black hole's quantum $N$-portrait~\cite{Dvali_2013_1}. The model is given in Eq.\ (\ref{eqn_model}) with slight changes in notation for the sake of brevity. It contains two sets of bosonic memory modes $K$ and $K'$, with the corresponding creation and annihilation operators $\hat{a}_k^{\dagger}$, $\hat{a}_k$ and $\hat{a}_{k'}^{\dagger}$, $\hat{a}_{k'}$, respectively, for $k^{(\prime)}=1, \dotsc, K^{(\prime)}$. The operators $\hat{a}_{k^{(\prime)}}^{(\dagger)}$ obey the standard commutation relations (here and throughout we set $\hbar \equiv 1$)
\begin{equation} \label{eqn_comm_rel}
\left[ \hat{a}_{j^{(\prime)}}, \hat{a}^{\dagger}_{k^{(\prime)}} \right] = \delta_{j^{(\prime)}k^{(\prime)}}, \quad \left[ \hat{a}_{j^{(\prime)}}, \hat{a}_{k^{(\prime)}} \right] = 0, \quad \left[ \hat{a}^{\dagger}_{j^{(\prime)}}, \hat{a}^{\dagger}_{k^{(\prime)}} \right] = 0.
\end{equation}
The corresponding occupation number operators are given by $\hat{n}_{k^{(\prime)}}=\hat{a}_{k^{(\prime)}}^{\dagger}\hat{a}_{k^{(\prime)}}$ with eigenvalues $n_{k^{(\prime)}}$ and eigenstates $\ket{n_{k^{(\prime)}}}$. Additionally, the model contains two more bosonic modes $\hat{n}_a$ (the master mode) and $\hat{n}_b$, with creation and annihilation operators $\hat{a}^{\dagger}$, $\hat{a}$ and $\hat{b}^{\dagger}$, $\hat{b}$, respectively. These operators satisfy commutation relations analogous to Eq.\ (\ref{eqn_comm_rel}) and allow the exchange of occupation number between the modes $\hat{n}_a$ and $\hat{n}_b$. The couplings $C_m$ and $C_b$ parametrize the interaction strength among the memory modes $\hat{n}_{k^{(\prime)}}$, and between the modes $\hat{n}_a$ and $\hat{n}_b$, respectively.

The effective energy gaps of the $\hat{n}_k$ and $\hat{n}_{k^{(\prime)}}$ modes are
\begin{equation}
\varepsilon_k \equiv \varepsilon \left( 1- \dfrac{\hat{n}_a}{N} \right) \quad \text{and} \quad \varepsilon_{k'} \equiv \varepsilon \left( 1- \dfrac{\hat{n}_a}{N-\Delta} \right),
\end{equation}
correspondingly. Therefore, the modes $K$ and $K'$ become effectively gapless for $n_a=N$ and $n_a=N-\Delta$, respectively. In the beginning of the system's time evolution, $n_a=N$ and the $K$ modes are gapless, while the $K'$ modes are not. The $K$ modes initially contain all of the information in terms of excitations, whereas $K'$ modes are unoccupied. However, as $n_a$ diminishes from $N$ to $N-\Delta$, a state with $n_a=N-\Delta$, $n_b=\Delta$ and where the excitations have been transferred to the $K'$ modes, becomes energetically available.

We denote the total occupation number of the two memory sectors by
\begin{equation}
N_m \equiv \sum\limits_{k=1}^K n_k + \sum\limits_{k'=1}^{K'} n_{k'}.
\end{equation}
Note that it is conserved. Following~\cite{Dvali_2020_2}, we additionally truncate all memory modes to qubits. We denote the states of the system as
\begin{equation}
\ket{n_a, n_b, n_1, \dotsc, n_K, n_{1'}, \dotsc, n_{K'}} \equiv \ket{n_a} \otimes \ket{n_b} \bigotimes_{k=1}^K \ket{n_k} \bigotimes_{k'=1}^{K'} \ket{n_{k'}}.
\end{equation}
The system is evolved from an initial state
\begin{equation} \label{eqn_in_st}
\ket{in} = \ket{\underbrace{N}_{a}, \underbrace{0}_{b}, \underbrace{\overbrace{1, 1}^{=N_m}, 0, \dotsc, 0}_{K}, \underbrace{0, \dotsc, 0}_{K'}}
\end{equation}
where $\hat{n}_a$ and $\hat{n}_b$ have occupation numbers $N$ and $0$, respectively, and only the first $N_m$ modes of the initially gapless $K$ memory sector contain one particle each. Following~\cite{Dvali_2020_2}, we set the basic energy unit $e \equiv 1$. The Hamiltonian of the model is then
\begin{equation} \label{eqn_model}
\begin{aligned}
\hat{H} &= \varepsilon \left( 1 - \dfrac{\hat{n}_a}{N} \right) \sum\limits_{k=1}^{K} \hat{n}_k + \varepsilon \left( 1 - \dfrac{\hat{n}_a}{N-\Delta} \right) \sum\limits_{k'=1}^{K'} \hat{n}_{k'} + C_b \left( \hat{a}^{\dagger} \hat{b} + \text{H.c.} \right) + \\
&+ C_m \left\lbrace \sum\limits_{k=1}^{K} \sum\limits_{k'=1}^{K'} f_1(k,k') \left( \hat{a}_k^{\dagger} \hat{a}_{k'} + \text{H.c.} \right) + \sum\limits_{k=1}^{K} \sum\limits_{l=k+1}^{K} f_2(k,l) \left( \hat{a}_k^{\dagger} \hat{a}_l + \text{H.c.} \right) + \right. \\
&+ \left. \sum\limits_{k'=1}^{K'} \sum\limits_{l'=k'+1}^{K'} f_3(k',l') \left( \hat{a}_{k'}^{\dagger} \hat{a}_{l'} + \text{H.c.} \right) \right\rbrace,
\end{aligned}
\end{equation}
with
\begin{equation}
f_i(k,l) = \begin{cases}
		   F_i(k,l)-1, & F_i(k,l) < 0.5 \\
		   F_i(k,l),   & F_i(k,l) \geq 0.5
		   \end{cases}
\end{equation}
and
\begin{equation}
F_i(k,l) = \left( \sqrt{2}(k+\Delta k_i)^3 + \sqrt{7}(l+\Delta l_i)^5 \right) \text{ mod } 1,
\end{equation}
where
\begin{equation}
\Delta k_1 = \Delta k_2 = 1, \quad \Delta k_3 = K+1, \quad \Delta l_1 = \Delta l_3 = K+1, \quad \Delta l_2 = 1.
\end{equation}

For a spherically symmetric system, such as a black hole, the free energy gap $\varepsilon$ of the memory modes can be found by labeling the states by quantum numbers $(l,m)$ of spherical harmonics $Y_l^m(\theta,\phi)$~\cite{Dvali_2020_2}. From the number of independent solutions for one degree $l$ of a given $Y_l^m(\theta,\phi)$, the energy of the highest mode and hence the free gap of the memory modes can be estimated as~\cite{Dvali_2020_2}
\begin{equation}
\varepsilon = \sqrt{K}.
\end{equation}
Additionally, we require that in the vicinity of the initial state in Eq.\ (\ref{eqn_in_st}), the $K'$ memory sector is not gapless. This results in the constraint $|\varepsilon_{k'}| \gg 1/\sqrt{N_m}$~\cite{Dvali_2020_2} or, equivalently,
\begin{equation}
\Delta \gg \dfrac{N}{1+\sqrt{N_m}\sqrt{K}}.
\end{equation}
Due to the considerable duration of the numerical simulations, for the majority of parameter scans in section~\ref{subsect_indiv} we choose a value for $\Delta$ that does not satisfy this constraint. However, from the scan over $\Delta$ of the sought-after quantities, we are, firstly, able to extract a clear dependence over the majority of the domain, and secondly, we observe that one of the analyzed quantities acquires a dissimilar value only when $\Delta$ is in the immediate vicinity of $N$.
The requirement that the memory modes remain effectively gapless imposes the constraints on the couplings~\cite{Dvali_2020_2}
\begin{equation} \label{eqn_Cb_Cm_cond}
C_b \lesssim \dfrac{1}{\sqrt{N}}, \quad C_m \lesssim \dfrac{1}{\sqrt{N_m}\sqrt{K}}.
\end{equation}
Furthermore, as in~\cite{Dvali_2020_2}, we set
\begin{equation}
\begin{aligned}
K' &= K.
\end{aligned}
\end{equation}
This leaves a set of six free and constrained parameters
\begin{equation} \label{eqn_param}
N, K, N_m, \Delta, C_b, C_m.
\end{equation}
In terms of the black hole's quantum $N$-portrait~\cite{Dvali_2013_1}, for a black hole of entropy $S \gg 1$, the parameters are further constrained as $N=S, K=S, N_m=S/2$~\cite{Dvali_2020_2}. This scaling is considered in section~\ref{subsect_BH}.

\section{Two effects}

Based on the black hole's quantum $N$-portrait~\cite{Dvali_2013_1}, a microscopic framework to describe scrambling was established in~\cite{Dvali_2013_2,Dvali_2013_3,Dvali_2015_1,Dvali_2015_2,
Dvali_2016_1,Dvali_2016_2} among other works (see further references therein). However, as stated in section~\ref{subsect_scr}, our results are independent of the choice of a specific measure of the uniformity of the state distribution, which is necessary for a precise definition of scrambling. Instead, the only criterion we require a scrambled system to have, is for the system's state to be sufficiently distributed (with respect to some measure) over the \emph{entire} Hilbert space of the system. This is a property shared by all definitions of scrambling. Given this broad characterization, we expect two effects to occur not later than the instant in time the system becomes scrambled:
\begin{itemize}
\item[(1)] The system's wavefunction $\ket{\psi(t)}$, being initially in only one basis state $\ket{in}$, spreads over the entire basis $\{ \ket{v_i} \}$, given a minimum state probability threshold.
\item[(2)] The most probable basis state that the system can be found in, changes from $\ket{in}$ to some other basis state $\ket{v_i} \neq \ket{in}$.
\end{itemize}

Let us specify these effects and the times in which they occur more precisely. Namely, over some finite, non-zero time period, the wavefunction $\ket{\psi(t)}$ diffuses over the entire particle basis $\{ \ket{v_i} \}$, given a finite, fixed and sufficiently small precision $p$ to which the probabilities $|C_i|^2 = |\braket{\psi(t)|v_i}|^2$ of the individual basis states are rounded. For this arbitrarily small but fixed minimum state probability threshold $p < 1/\text{dim} \left( \{ \ket{v_i} \} \right)$, we denote the fraction of the basis that $\ket{\psi(t)}$ spreads over by $f$, and the time when $\ket{\psi(t)}$ spreads over the entire basis, i.e.\ the earliest time when $f=1$, by $t_f$. Concretely, we define $f$ as
\begin{equation}
f \equiv \dfrac{1}{\mathcal{N}} \sum\limits_{i=1}^{\mathcal{N}} H(|C_i|^2 - p)~, \quad
H(x) = \begin{cases}
	   0, & x \leq 0 \\
	   1, & x > 0
	   \end{cases},
\end{equation}
where $\mathcal{N} \equiv \text{dim} \left( \{ \ket{v_i} \} \right)$ is the number of basis states, i.e.\ the dimension of the Hilbert space of the system, and $H(x)$ is the left-continuous Heaviside step function. Clearly, for the model in Eq.\ (\ref{eqn_model}), $p$ has to be non-zero since the moment the system begins evolving, $\ket{\psi(t)}$ diffuses instantaneously over the entire basis, but with very small contributions of some states. However, at the time $t_f$, all of these probabilities increase above the given threshold of $p$.

Additionally, we denote by $t_c$ the time when the most probable state of the system changes from the initial basis state $\ket{in}$ to some other basis state $\ket{v_i} \neq \ket{in}$, i.e.\ the earliest time when $\underset{i}{\max} (|C_i|^2) \neq |C_{in}|^2$.

The above two effects are investigated numerically in section~\ref{sect_num} for the model in Eq.\ (\ref{eqn_model}). Our main goal is to understand how the two effects and their respective times $t_f$ and $t_c$ depend on the system's parameters in Eq.\ (\ref{eqn_param}). Note that the above definitions for the two effects of spreading involve quantitative measures of the quantum state's distribution in terms of probabilities. Due to this, our definitions differ somewhat to diffusion described in~\cite{Susskind_2011}. Furthermore, our definitions also differ from that of \emph{minimal scrambling} in~\cite{Dvali_2016_2}. Specifically, minimal scrambling for a system of $N$ qubits in~\cite{Dvali_2016_2} is defined as when the state of that system is a uniform superposition of $N$ basis states, each with $1/N$ probability, as opposed to a superposition of all $2^N$ basis states.

\section{Numerical results} \label{sect_num}

The numerical simulations are performed with QuSpin version $0.3.4$~\cite{QuSpin_1,QuSpin_2,QuSpin_3}. Unless stated otherwise, and except for the parameter(s) varied, the parameter values we use are
\begin{equation}
p=5 \cdot 10^{-21}, \quad N=4, \quad K=4, \quad N_m=2, \quad \Delta=1, \quad C_b=0.1, \quad C_m=0.1.
\end{equation}
When varying the parameters, we restrict their possible ranges to respect the constraints in Eq.\ (\ref{eqn_Cb_Cm_cond}). By the value chosen for $p$ we imply that we round each individual $|C_i|^2$ to $20$ decimal places, and if the result is $>0$, we add its contribution to $f$. The reason that the precision $p$ needs be sufficiently small is that for $N$ and/or $K$ large enough, $t_f$ quickly diverges from growing logarithmically and rapidly increases. However, we expect that for an arbitrary but fixed range of values for $N$ and/or $K$ there exists a sufficiently small precision $p$ such that $t_f$ behaves logarithmically on the given parameter range. Figure~\ref{fig_typ} depicts an exemplary plot of the two effects for the chosen parameter values.
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.75]{typ}
  \caption{\small Left (blue): fraction of basis $f$ that the state $\ket{\psi(t)}$ is spread over. Basis includes $140$ states in total. Right (red): maximum contribution $|C_i|^2$ of any single basis state $\ket{v_i}$ with $i \in [1,140]$, and (green, dot-dashed): contribution of the initial state $|C_{in}|^2$. The time step is $10^{-2}$. Note that the time scale is logarithmic. In this example $t_f=0.32$ and $t_c=6.25$.}
  \label{fig_typ}
\end{figure}

In several simulations we observe that after reaching $1$, the value of $f$ can decrease back below $1$ for a period of time, but then increases to $1$ and remains at this value. However, in most simulations, once $f=1$ is reached, it tends to stay at $1$. The maximal contribution $\underset{i}{\max} (|C_i|^2)$ can also come back to being equal to $|C_{in}|^2$, but then changes again to that of some other state. We also observe that the time $t_f$ does not always precede $t_c$; this is clear as the value of $f$ strongly depends on the chosen precision $p$.

We expect large values of $N$, $K$, $N_m$ (i.e.\ large systems with large Hilbert spaces) and $\Delta$ (i.e.\ large effective energy gap difference), and small values of $C_b$, $C_m$ (i.e.\ weak couplings) to slow down the spread of the wavefunction $\ket{\psi(t)}$. The goal of the numerical simulations is first, to verify and second, to quantify this behavior. Moreover, we expect the diffusion of $\ket{\psi(t)}$ as defined above in terms of the fraction $f$ of the basis to be a crucial pre-requisite for scrambling. This stems from a simple fact that before the system becomes scrambled over the entire basis with respect to any measure, its state has to firstly spread over the entire basis with some small non-zero contribution of every basis state. Furthermore, we expect there to be no ``diffusers'', ``spreaders'', ``distributors'' or ``pre-scramblers'' that are faster than the model in Eq.\ (\ref{eqn_model}), since in this system every mode of each sector is directly coupled to every other mode in that sector. It remains to investigate, how fast this ``fast pre-scrambler'' is.

\subsection{Individual parameter scaling} \label{subsect_indiv}

We first vary the parameters of the system individually. The results are shown in Fig.\ \ref{fig_indiv} with the respective best obtained fit functions. The corresponding fit parameters are listed in Tab.\ \ref{tab} in the Appendix. We observe that indeed, large values of $N$, $K$, $N_m$, $\Delta$ and small values of $C_b$, $C_m$ slow down the diffusion in terms of $t_f$. Moreover, we obtain new quantitative information from the dependence of $t_f$ and $t_c$ on the model parameters. However, for the dependence of $t_c$ on $N_m$ no behavior could be extracted due to irregularity of the data points. Furthermore, note that for the dependence of $t_f$ on $N_m$ the standard errors of the free fit-model parameters are comparable in magnitude with the values of the parameters themselves, and that the adjusted coefficient of determination $\overline{R}^2$ and the adjusted root-mean-square error are the lowest and highest, respectively, when compared to all other fits. We therefore can not consider this to be a good fit to the data. Unfortunately, due to numerical limitations, we weren't able to obtain more data in a reasonable amount of time.

However, for higher values of $N_m$ we expect $t_f$ to increase logarithmically with $N_m$ for the following two reasons. Firstly, note that $t_f$ doesn't increase logarithmically also for small values of $N$. The logarithmic behavior becomes evident only for higher values of $N$. Secondly, there are only three parameters in the model, namely $N$, $K$ and $N_m$, that affect the dimension of the Hilbert space. We expect the size of the Hilbert space to be the main factor that controls the speed at which the wavefunction spreads, and hence also the behavior of $t_f$. We therefore expect $t_f$ to depend on $N$, $K$ and $N_m$ in identical manner.
\begin{figure}[htbp]
    \centering
    \begin{subfigure}[b]{0.45\textwidth}
    	\centering
        \includegraphics[width=\textwidth]{N}
        \caption{\footnotesize $t_f \sim \ln(N)$ for $N \geq 20$, \\ $t_c \sim N^{-0.966 \pm 0.044}$ for $N \geq 8$.}
    \end{subfigure}
    \begin{subfigure}[b]{0.45\textwidth}
        \includegraphics[width=\textwidth]{K}
        \centering
        \caption{\footnotesize $t_f \sim \ln(K)$, \\ $t_c \sim ~ \mathrm{e}^{(-0.083 \pm 0.019) K}$ for $K \geq 19$.}
    \end{subfigure}
    \begin{subfigure}[b]{0.45\textwidth}
        \includegraphics[width=\textwidth]{Nm}
        \caption{\footnotesize $t_f \sim ~ \mathrm{e}^{(0.43 \pm 0.13) N_m}$, \\ $t_c$ scaling unknown.}
    \end{subfigure}
    \begin{subfigure}[b]{0.45\textwidth}
        \includegraphics[width=\textwidth]{DLT}
        \caption{\footnotesize $t_f \sim \mathrm{e}^{(2.373 \pm 0.087) \Delta}$ for $\Delta \leq 3.986$, \\ $t_c \sim \mathrm{e}^{(-0.323 \pm 0.019) \Delta}$ for $\Delta \geq 0.5$.}
    \end{subfigure}
    \begin{subfigure}[b]{0.45\textwidth}
        \includegraphics[width=\textwidth]{Cb}
        \caption{\footnotesize $t_f \sim C_b^{-1.177 \pm 0.059}$ for $C_b \geq 0.0057$, \\ $t_c \sim C_b^{-4.41 \pm 0.48}$ for $C_b \geq 0.25207$.}
    \end{subfigure}
    \begin{subfigure}[b]{0.45\textwidth}
        \includegraphics[width=\textwidth]{Cm}
        \caption{\footnotesize $t_f \sim C_m^{-0.3847 \pm 0.0051}$, \\ $t_c \sim C_m^{-1.0734 \pm 0.0044}$.}
    \end{subfigure}
    \caption{\small The individual subfigures show the scaling of $t_f$ and $t_c$ with the individual parameters $N$, $K$, $N_m$, $\Delta$, $C_b$, $C_m$. Note the different parameter ranges for the top and bottom, and left and right axes. The parameters altered are: (a) for only the $t_f$ scan: $p=5 \cdot 10^{-41}$, additionally, $N$ is varied up to $150$ for a better fit; the permitted range is $N \leq 100$. The points corresponding to $N \in [101,150]$ and the respective fit are marked green. (c) for both $t_f$ and $t_c$ scans: $K=8$.}
    \label{fig_indiv}
\end{figure}

Nevertheless, we are able to extract the scaling of $t_f$ and $t_c$ with respect to all other parameters. Our main finding is that
\begin{equation}
t_f \sim \ln(K).
\end{equation}
That is, the time $t_f$ for the initial state of the system to spread over the entire basis $\{ \ket{v_i} \}$ with $|C_i|^2 > p \; \forall i \in [1,\mathcal{N}]$ for a fixed $p$ scales logarithmically with the number of modes $K$. Since the model we are investigating possesses all-to-all couplings between the modes of each sector, we expect the above logarithmic scaling of the wavefunction's ``diffusion'', ``spreading'', ``distribution'' or ``pre-scrambling'' to be the fastest possible for any physical system.
\bigskip
\\
\makebox[\textwidth][c]{%
\begin{minipage}{0.8\textwidth}
We therefore define a system to be a \emph{fast pre-scrambler}, if for an arbitrarily large but fixed range of number of degrees of freedom $K$ there exists a sufficiently small minimum state probability threshold $p<1/\text{dim}(\{ \ket{v_i} \})$, such that the state of the system $\ket{\psi(t)}$, initially in one basis state $\ket{in}$, spreads over the entire Hilbert space $\{ \ket{v_i} \}$ into a superposition of all $\text{dim}(\{ \ket{v_i} \})$ basis states each with a probability $|C_i|^2 = |\braket{\psi(t)|v_i}|^2 > p$, in a time logarithmic in $K$.
\end{minipage}}
\bigskip
\\
There are two reasons why we name this effect ``pre-scrambling''. Firstly, we want to distinguish it from diffusion in~\cite{Susskind_2011} and minimal scrambling in~\cite{Dvali_2016_2}. Secondly, since the probability threshold $p$ can be set arbitrarily small, this type of the state's spreading will precede scrambling where the state is required to be sufficiently distributed with respect to some measure. One could, of course, define a scrambling measure identical to our pre-scrambling measure. However, in that case there is a possibility that while the state does spread over the entire basis to some extent, the most probable state of the system is still its initial state. In other words, this can occur if $t_f < t_c$. In this case it is questionable whether or not we can call the state of the system truly scrambled.

\subsection{Application to black holes} \label{subsect_BH}

Following the black hole's quantum $N$-portrait~\cite{Dvali_2013_1}, the parameters are constrained as mentioned in section~\ref{sect_model}. Unfortunately, an extensive scan with the corresponding scaling of parameters is very time consuming. However, we can still argue what we expect the behavior of $t_f$ and $t_c$ to be in this regime. Firstly, as was discussed previously, we expect $t_f$ to depend logarithmically on $N$, $K$ and $N_m$. Secondly, we performed scans of $t_f$ and $t_c$ for $N=K \leq 15$ with $N_m$ fixed constant and equal to $2$, with precisions $p=5 \cdot 10^{-41}$ and $p=5 \cdot 10^{-21}$, respectively. From these scans we observe that both $t_f$ and $t_c$ each scale as their individual $N$ scans.

Moreover, as is argued in~\cite{Hayden_2007,Sekino_2008}, black holes are fast scramblers, and are expected to scramble in a time logarithmic in the number of degrees of freedom. Based on the discussion above, we expect pre-scrambling to occur not later than scrambling, irrespective of the specific measure used to define the latter. The model in Eq.\ (\ref{eqn_model}) possesses direct all-to-all couplings between the modes of every sector. We therefore expect it to have the fastest pre-scrambling time possible in terms of the number of degrees of freedom. We have argued above that for this model pre-scrambling should occur in a time logarithmic in the number of degrees of freedom when varying $N$, $K$ and $N_m$ individually. From the arguments above it is then natural to expect the model in question to pre-scramble in a time logarithmic in $N$ for the above black hole parameter scaling regime.

\section{Conjectures}

Summarizing the results of the investigation above, we arrive at the following conjectures:
\bigskip
\\
\makebox[\textwidth][c]{%
\begin{minipage}{0.8\textwidth}
\begin{itemize}
\item[(1)] The fastest pre-scramblers spread the initial state of the system over its entire Hilbert space in a time logarithmic in the number of degrees of freedom, given a minimum state probability threshold.
\item[(2)] The enhanced memory capacity model in Eq.\ (\ref{eqn_model}) is a fast pre-scrambler.
\item[(3)] The time of (fast) pre-scrambling is $\leq$ than that of (fast) scrambling, irrespective of the choice of a specific measure of the uniformity of the state distribution for the latter.
\item[(4)] Consequently, fast scramblers are fast pre-scramblers.
\item[(5)] In particular, black holes are fast pre-scramblers.
\end{itemize}
\end{minipage}}
\bigskip
\\
Note, however, that we do \emph{not} claim any of the following: we do not claim that all fast pre-scramblers are fast scramblers, nor do we claim that all fast pre-scramblers have all-to-all couplings between the degrees of freedom. Furthermore, the above are conjectures and should by no means be regarded as concrete proofs.

\section{Conclusion}

In addition to formulating the above conjectures, we would like to conclude by several complementary comments. First, note that the diffusion of the wavefunction as measured by $t_f$ possesses an at least exponential dependence on the energy gap difference $\Delta$, and for $\Delta$ close to $N$ increases even more rapidly than for lower values of $\Delta$. Second, we find it encouraging and fascinating that even the simple prototype model of a system with enhanced memory capacity in Eq.\ (\ref{eqn_model}) exhibits properties that we expect black holes to possess. As an additional remark on the prototype model in Eq.\ (\ref{eqn_model}), note that the all-to-all direct coupling of modes is in full agreement with the description of a black hole as a fast scrambler in~\cite{Susskind_2011}. Third, we allow ourselves to speculate that the Sachdev-Ye-Kitaev model~\cite{Sachdev_1993,Kitaev_2015} is a fast pre-scrambler, since it possesses all-to-all direct couplings and was shown to be a fast scrambler.

We also note that from the numerical data we are able to observe that a large amount of quantum information stored in the system's state slows down the process of the system's pre-scrambling via the memory burden effect.

As an outlook, let us note that the temperature dependent scrambling as discussed in~\cite{Sekino_2008,Shenker_2014,Kitaev_2015,Maldacena_2016} will be investigated for the model in Eq.\ (\ref{eqn_model}) in a future work~\cite{scr_to_appear}. Finally, the scrambling for the model in Eq.\ (\ref{eqn_model}) in relation to the eigenstate thermalization hypothesis~\cite{Deutsch_1991,Srednicki_1994} will also be considered in a future work~\cite{eth_to_appear}.

\section*{Acknowledgments}

We thank Gia Dvali for insightful discussions and valuable comments.

\begin{appendix}

\section*{Appendix}

\begin{table}[H]
\centering
\begin{adjustbox}{width=\textwidth}
\begin{tabular}{l|l|l|l|l|l|l|l|l}
$t_f$ & $\overline{R}^2$ & $\text{RMSE}$ & $a$ & $\sigma_{a}$ & $b$ & $\sigma_{b}$ & $c$ & $\sigma_{c}$ \\
\hline
$a + b\ln(N)$ & $9.9991 \cdot 10^{-1}$ & $6.00 \cdot 10^{-2}$ & $-9.645$ & $4.3 \cdot 10^{-2}$ & $3.6420$ & $9.8 \cdot 10^{-3}$ & $\text{n/a}$ & $\text{n/a}$ \\
$a + b\ln(K)$ & $9.88 \cdot 10^{-1}$ & $4.80 \cdot 10^{-2}$ & $1.97 \cdot 10^{-1}$ & $2.8 \cdot 10^{-2}$ & $7.92 \cdot 10^{-2}$ & $8.8 \cdot 10^{-3}$ & $\text{n/a}$ & $\text{n/a}$ \\
$a + b\exp(cN_m)$ & $9.77 \cdot 10^{-1}$ & $2.86 \cdot 10^{-1}$ & $-1.3 \cdot 10^{-1}$ & $4.3 \cdot 10^{-1}$ & $1.9 \cdot 10^{-1}$ & $1.9 \cdot 10^{-1}$ & $4.3 \cdot 10^{-1}$ & $1.3 \cdot 10^{-1}$ \\
$a + b\exp(c\Delta)$ & $9.9996 \cdot 10^{-1}$ & $2.44 \cdot 10^{-3}$ & $3.183 \cdot 10^{-1}$ & $1.0 \cdot 10^{-3}$ & $7.7 \cdot 10^{-6}$ & $2.7 \cdot 10^{-6}$ & $2.373$ & $8.7 \cdot 10^{-2}$ \\
$a + b C_b^c$ & $9.98 \cdot 10^{-1}$ & $7.81 \cdot 10^{-2}$ & $1.76 \cdot 10^{-1}$ & $3.8 \cdot 10^{-2}$ & $7.2 \cdot 10^{-3}$ & $2.2 \cdot 10^{-3}$ & $-1.177$ & $5.9 \cdot 10^{-2}$ \\
$a + b C_m^c$ & $9.9993 \cdot 10^{-1}$ & $6.10 \cdot 10^{-3}$ & $5.34 \cdot 10^{-2}$ & $7.5 \cdot 10^{-3}$ & $1.086 \cdot 10^{-1}$ & $4.1 \cdot 10^{-3}$ & $-3.847 \cdot 10^{-1}$ & $5.1 \cdot 10^{-3}$

\bigskip
\bigskip \\

$t_c$ & $\overline{R}^2$ & $\text{RMSE}$ & $a$ & $\sigma_{a}$ & $b$ & $\sigma_{b}$ & $c$ & $\sigma_{c}$ \\
\hline
$a + b N^c$ & $9.996 \cdot 10^{-1}$ & $4.52 \cdot 10^{-2}$ & $9.25 \cdot 10^{-1}$ & $5.9 \cdot 10^{-2}$ & $22.9$ & $1.9$ & $-9.66 \cdot 10^{-1}$ & $4.4 \cdot 10^{-2}$ \\
$a + b\exp(c K)$ & $9.94 \cdot 10^{-1}$ & $3.97 \cdot 10^{-1}$ & $3.54$ & $3.2 \cdot 10^{-1}$ & $19.2$ & $6.3$ & $-8.3 \cdot 10^{-2}$ & $1.9 \cdot 10^{-2}$ \\
$a + b\exp(c\Delta)$ & $9.99994 \cdot 10^{-1}$ & $1.69 \cdot 10^{-2}$ & $8.205$ & $7.2 \cdot 10^{-2}$ & $-2.662$ & $5.3 \cdot 10^{-2}$ & $-3.23 \cdot 10^{-1}$ & $1.9 \cdot 10^{-2}$ \\
$a + b C_b^c$ & $9.991 \cdot 10^{-1}$ & $7.01 \cdot 10^{-2}$ & $9.72 \cdot 10^{-1}$ & $7.2 \cdot 10^{-2}$ & $3.9 \cdot 10^{-3}$ & $2.7 \cdot 10^{-3}$ & $-4.41$ & $4.8 \cdot 10^{-1}$ \\
$a + b C_m^c$ & $9.9998 \cdot 10^{-1}$ & $5.55 \cdot 10^{-2}$ & $3.0 \cdot 10^{-2}$ & $4.7 \cdot 10^{-2}$ & $5.288 \cdot 10^{-1}$ & $9.4 \cdot 10^{-3}$ & $-1.0734$ & $4.4 \cdot 10^{-3}$
\end{tabular}
\end{adjustbox}
\caption{Tables of best obtained fit functions for $t_f$ and $t_c$ on the top and bottom, respectively. The corresponding coefficients of determination $\overline{R}^2$ and unbiased root-mean-square errors (RMSE) are both adjusted for the number of fit-model parameters. The values of the free fit-model parameters are given with the corresponding standard errors. No dependence of $t_c$ on $N_m$ could be extracted from the corresponding data points.}
\label{tab}
\end{table}

\end{appendix}


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

While gamma-ray bursts (GRBs) are well known to emit in the energy band keV~$ \leq
\epsilon_\gamma\leq$~GeV, in recent years there are accumulating observations at
even higher energies, in the range of several GeV to TeV. Early detections in this
band date back to
the 1990's. The Energetic Gamma-Ray Experiment Telescope (EGRET) on-board CGRO,
was the first instrument to monitor a photon with an energy of
$\sim18\hspace{0.3em}$GeV in coincidence with GRB 910503 \citep{1992A&A...255L..13S, 1993ApJ...413..281B,
1994Natur.372..652H}, followed by a $\sim1\hspace{0.3em}$TeV photon
observation associated with GRB 970417a in the Milagrito field of
view  \citep{2003ApJ...583..824A}. Successive detections of 14 bursts/year,
on average, have been done by \textit{Fermi}-LAT (Large Area Telescope) in
a broad range of energies up to $\sim 100$~GeV \citep{2018IJMPD..2742003N}.
Very high sub-TeV band detections were recently
reported by the ground-based Cherenkov
observatories, namely the High Energy Stereoscopic System  (HESS) and the Major Atmospheric
Gamma Imaging Cherenkov (MAGIC). HESS is actively detecting GRBs
at energies larger than $100$~GeV, e.g., the 100 to 440 GeV photons reported in  GRB
180720B \citep{2019Natur.575..464A}. Similarly, MAGIC is
also exploring the window of VHE emissions with the detection of
GRB 190114C in the sub-TeV band \citep{MAGIC:2019irs}.

These recent discoveries have  called for much attention
as GRB afterglow spectra were long predicted to have VHE
signatures of this kind \citep{1994ApJ...432..181M, Dermer_2000, 2002ARA&A..40..137M, 2004RvMP...76.1143P, 2015PhR...561....1K}.  Furthermore, such detections of GRB
photons with energy larger than 1 TeV are expected to become ubiquitous in the near
future, once the upcoming new generation high sensitive $\gamma-$detector Cherenkov
Telescope Array (CTA) begins operations with its unprecedented precision \citep{knodlseder2020cherenkov}.

Of particular interest is the recent detection of the very high energy ($0.2-1\hspace{0.3em}$
TeV) emission from GRB 190114C by MAGIC. Thanks to its proximity, this burst turned to be
one of the brightest bursts ever detected. It was observed by diverse space
observatories from near the Fornax constellation 4.5 billion light years away  with
a redshift of $z=0.4245 \pm 0.0005$ \citep{web:nasa:23708, web:nasa:stats} in a dense environment right in the middle of the luminous galaxy.
It was reported with the prompt phase isotropic equivalent energy $E_{\rm
iso} \approx 2.5 \times 10 ^{53}$~erg
for a duration of $T_{\rm 90}\sim 25$~s.
The afterglow of this burst was observed in an unprecedentedly large energy band for
an unprecedented duration. The broadband  afterglow spectrum can be split into a \textit{low-energy} component ranging from optical
to X-rays and a \textit{high-energy} component (MeV to GeV) extended
by a non-thermal VHE tail (at least up to the TeV band)
\citep{1997Natur.389..261F, 1997Natur.386..686V, 2000ApJ...543..722L, zhang_2018}.
The low energy spectral component is
usually observed to be a broken power-law
\citep{costa1997discovery, van1997transient, wijers1997shocked, harrison1999optical,
2000ApJ...529..635B}. It can be interpreted in the framework of synchrotron emission
from shock accelerated relativistic electrons which naturally leads to a broken
power-law spectrum \citep{1993ApJ...418L...5P, 1997ApJ...476..232M, Sari_1998}. However,
the high-energy emission is too bright to be explained within the framework of this
basic synchrotron model \citep{1998ApJ...501..772P, 2005ApJ...633.1018P}. For this reason,
it is required to extend the basic afterglow model based
on electron synchrotron to explain the emission at the highest observed energies.

Various theoretical models were proposed to explain the observed signal at the
high-energy end of the spectrum. A leading suggestion is the synchrotron
self-Compton (hereinafter SSC) mechanism, where the low-energy synchrotron photons
inverse-Compton scatter on the electrons that previously emitted them
\citep{ ghisellini1998quasi, Dermer_1999ApJ, Dermer_2000, 2001ApJ...548..787S, 2009ApJ...703..675N, liu2013interpretation, 10.1093/mnras/stw1175, fraija2019synchrotron, Derishev:2021ivd}.
Alternatively, it was suggested that the high and very-high energy emission might
be produced by (co-)accelerated baryons. Relativistic hadronic-induced emission
processes such as proton-synchrotron,
photo-pion and photo-pair processes and emission from their secondaries were previously suggested
as viable candidates; see \textit{e.g.} \cite{1998ApJ...499L.131B, 2009ApJ...699..953A, 2010OAJ.....3..150R, 2022arXiv220901940G}.

Direct indication of proton acceleration to high energies is given by the
explicit detection of high-energy cosmic-rays \citep{asakimori1998cosmic,
sanuki2000precise, lipari2020shape} and of astrophysical PeV
neutrinos by Icecube \citep{2013PhRvL.111b1103, 2015MNRAS.453..113, 2017ARNPS..67...45}.
Their potential association with blazars \citep{PhysRevD.99.103006, 2012MNRAS.426..462P} indicates that relativistic jets might harbor
high-energy protons, further supporting hybrid and hadronic models for the
high-energy emission; see \textit{e.g.} \cite{2001APh....15..121M, 2002MNRAS.332..215A,
2003APh....18..593M, 2022MNRAS.509.2102G}. Such high energy protons
could also be present in GRB jets, emitting high energy synchrotron photons \citep{PhysRevLett.78.4328, 1998ApJ...totania, 2001ApJ...559..110Z, 2007MNRAS.380...78G, 2010OAJ.....3..150R}. As we show here, the spectral theoretical expectations from the proton-synchrotron model are consistent with the collective data available for GRB 190114C. We therefore suggest that proton synchrotron may be the leading radiative emission that shape the spectra at the very high ($\sim$ TeV) emission range.

This paper is organised as follows. In section \ref{sec:Model} we briefly review the hydrodynamical properties of the gas behind the relativistically propagating shock, within the assumptions set by the fireball evolution framework. We then recall the dependence of the Lorentz factor, $\Gamma$, the radius, $r$ and the magnetic field, $B$ in terms of the independent parameters.
Following that, we describe the characteristics of the synchrotron radiation model of the relativistic electrons and protons in section  \ref{sec:sync}.
Constraints on the model parameters are set from the observations of the afterglow of GRB 190114C in two discrete time intervals in section \ref{sec:Constraints}.
In section \ref{sec:RESULTS}, we present our results and display the corresponding spectral energy distribution functions under different assumptions. In particular, we explore three special cases in subsections \ref{sec:case1},
\ref{sec:case2} and \ref{best fit}, all linked to uncertainty regarding the characteristics of the particle acceleration mechanism. We discuss contribution of various addition radiative processes in section \ref{sec:5}, and show that they are sub-dominant for the parameters we consider. In section \ref{sec:discussion}, we discuss our results in light of previous works. The conclusions follow in section \ref{sec:Conclusion}.

\section{Afterglow model: dynamics and emission} \label{sec:Model}

\subsection{Dynamics}

We adopt the standard fireball model scenario \citep{1986ApJ...308L..43P,  1986ApJ...308L..47G, 1990ApJ...348..485P,Mesazoras1993ApJ...405..278M, Piran:1993b, 1998ApJ...499..301M}, in which the afterglow emission is explained by the interaction between a highly relativistic jet and its surrounding medium \citep{1993ApJ...418L...5P, Sari_1998}.
We hereby focus on the self-similar expansion phase, which is expected after a few to a few tens of seconds, during which a relativistically expanding forward shock collects and heats the material from the surrounding interstellar medium (ISM). This heated material, in turn, radiates thereby produces the observed afterglow signal \citep{10.1093/mnras/263.4.861, 1995ApJ...455L.143S, Sari_1998, 2000ApJ...543...66P}.

When an expanding shell of ultra-relativistic plasma having a Lorentz factor $\Gamma$
interacts with the cold ISM (assumed here to have a  constant density $n$) at rest,
a collisionless shock is formed, and propagates into the ISM.
As the cold ISM material crosses the shock front, its (comoving) density increases to  $4\Gamma n$ and its internal (comoving) energy density becomes $e = 4 \Gamma^2 n m_p c^2$, where $m_p$ is the proton mass and $c$ is the speed of light.
The time evolution of the Lorentz factor and radius of the shock, $\Gamma$ and $r$, respectively, are given by the self-similar solution of a relativistic blast wave expansion  derived by \cite{1976PhFl...19.1130B},
\begin{align}
   \Gamma (E, n; t) &= \left[ \frac{17 E (1+z)}{1024 \pi n m_p c^5 t^{3}} \right]^{1/8} = 9.4 \hspace{0.3em} E_{53}^{1/8} t_{day}^{-3/8} n_0^{-1/8}, \\
    r(E,n; t) &= \left[\frac{17 E t }{4 \pi n m_p c (1+z)}\right]^{1/4} = 6.4 \times 10 ^{17}\hspace{0.3em} E_{53}^{1/4} t_{day}^{1/4} n_0^{-1/4} \quad\text{cm},
\end{align}
where $n = 1 \times n_0 \hspace{0.3em} \text{cm}^{-3}$, $t_{day}$ is the observed
time in days, and the isotropic equivalent kinetic energy is $E = 10^{53} E_{53} \hspace{0.3em} \text{erg}$. Hereinafter any quantity $X_n$ is such that $X = 10^n X_n$ and cgs units are adopted. In addition, here and in all numerical expressions below, we
consider the redshift to be $z=0.4245$, \textit{i.e.} the observed cosmological redshift of GRB 190114C \citep{web:nasa:stats, web:nasa:23708}. The comoving shell-expansion time (the dynamical timescale) is $t_{dyn} \cong \Gamma t$.

We further use the standard assumption that the magnetic field is generated by the shock, and carries an unknown fraction  $\epsilon_B$ of the internal energy density behind the shock front. As a consequence, in the co-moving frame, the magnetic field strength is
\begin{align}
   B &= \sqrt{32\pi \epsilon_B \Gamma^2 n m_p c^2} =  0.4 \hspace{0.3em} E_{53}^{1/8} \epsilon_{B,-2}^{1/2} t_{day}^{-3/8} n_0^{3/8}\text{~G},
\end{align}

where we took $\epsilon_B = 10^{-2}\epsilon_{B, -2}$.

\subsection{Synchrotron Emission Process} \label{sec:sync}

We assume that both electrons and protons are accelerated by the propagating shock wave, and attain a power-law  distribution $N(\gamma) d\gamma \propto \gamma^{-p} d\gamma$ in the range $\gamma_{\min}\leq \gamma \leq \gamma_{\max}$ as they reach the downstream region. Here, $p$ is the post-shock spectral index which we consider is similar for both the proton and electron populations. Typically, $p$ is assumed to be $p>2$, which is in agreement with both theoretical \citep{1996ApJ...473..204S, PhysRevLett.80.3911, 2000ApJ...542..235K, 2011ApJ...726...75S} and observational analysis of GRB afterglows \citep{2006MNRAS.371.1441S}.
Let $\gamma_p$ ($\gamma_e$) be the Lorentz factor of a single proton (electron). As the particles radiate, they cool. The radiative cooling time of a proton having Lorentz factor $\gamma_p$ is
\begin{align}
t_{p, \rm cool}= \frac{\gamma_p m_p c^2}{\frac{4}{3}\sigma_T c \frac{m_e^2}{m_p^2}\gamma_p^2 \frac{B^2}{8\pi}}  =\frac{ 6\pi m_p^3 c}{\sigma_T m_e^2\gamma_p B^2}, \label{eq:proton_cooling_synch}
\end{align}
and that of the electron is $t_{e, \rm cool}  = ({m_e^3}/{m_p^3}) (\gamma_p/\gamma_e) t_{p, \rm cool} $  \citep{george1979radiative}. Here, $\sigma_T$ is the Thomson cross-section.
The energy dependence of the cooling time implies the existence of a characteristic Lorentz factor, denoted by $\gamma_c$, for which the cooling time is equal to the dynamical time.
There are two categories of solution defined by the values of $\gamma_{\min}$ and $\gamma_{c}$: the fast cooling regime, where $\gamma_c < \gamma_{\min}$, and the slow cooling regime, where $\gamma_{\min} < \gamma_c$.

To find the value of $\gamma_{\min}$, we assume that a fraction $\xi_p$ ($\xi_e$) of the proton (electron) population is injected to a power law distribution.
We further assume that the protons (electrons) accelerated to the power-law distribution carry a fraction $\epsilon_p$ ($\epsilon_e$) of the available energy released as internal energy by the shock.
Since the average energy per particle in the downstream region is $\Gamma m_p c^2$, and as the average Lorentz factor of an energetic (power-law distributed) particle is $\bar{\gamma} =  [({p-1})/({p-2})] \gamma_{\min} \equiv f(p)^{-1} \gamma_{\min}$ (assuming $p>2$), one finds
the minimum injection Lorentz factor of the accelerating protons and electrons to be
\begin{align}
   \gamma_{e, \min} &= f(p) \left[\frac{m_p}{m_e}\right]\left[\frac{\epsilon_e}{\xi_e}\right] \Gamma = 2 \times 10^3 \hspace{0.3em} f(p) \xi_e^{-1}E_{53}^{1/8} n_0^{-1/8} t_{day}^{-3/8} \epsilon_{e,-1}, \\
    \gamma_{p, \min} &= f(p)\left[\frac{\epsilon_p}{\xi_p}\right]\Gamma \gtrsim 1 \hspace{0.3em} f(p) \xi_p^{-1}E_{53}^{1/8} n_0^{-1/8} t_{day}^{-3/8}\epsilon_{p,-1}.
\end{align}
Here and below, we consider the number density of the accelerated electrons, $n_e = \xi_e n$ and of the accelerated protons,  $n_p = \xi_p n$.

By balancing the co-moving dynamical time
$t_{dyn}$ and the cooling time by synchrotron emission $t_{cool}$, the cooling
Lorentz factors of the particle species are $ \gamma_{p,c} = ({6\pi m_p^3 c})/({\sigma_Tm_e^2B^2\Gamma t})$
and $\gamma_{e,c} = ({m_e}/{m_p})^3\gamma_{p,c}$.
A third characteristic Lorentz factor is the maximum Lorentz factor achieved by the particles \citep{1996ApJ...457..253D}. It can be estimated as follows. As the particles cool by the synchrotron process, their maximum achievable Lorentz factors are
calculated by equating the energy loss time and the particle
acceleration time, $t_{\rm acc} = 2 \pi \alpha \epsilon/qcB $
\footnote{Note that $\epsilon/ (cqB) = r_L$ is the
Larmor radius.
Here, the acceleration time $t_{\rm acc}$ is for Fermi-type accelerations.}
where $\epsilon = \gamma m c^2$ is the particle's energy, and $\alpha \geq 1$ is a numerical coefficient prescribing the acceleration efficiency. The
maximum Lorentz factors are thus
\begin{align}
   \gamma_{e,\max} &= \left[\frac{6\pi q}{\alpha \sigma_T B} \right]^{1/2} = 2 \times 10^8 \hspace{0.3em} \alpha^{-1/2} E_{53}^{-1/16} n_0^{-3/16} t_{day}^{3/16} \epsilon_{B,-2}^{-1/4} \\
   \text{and}
   \qquad\gamma_{p,\max} &= \left[\frac{6 \pi qm_p^2}{\alpha m_e^2 \sigma_T B} \right]^{1/2} = 3.6 \times 10^{11} \hspace{0.3em} \alpha^{-1/2} E_{53}^{-1/16} n_0^{-3/16} t_{day}^{3/16} \epsilon_{B,-2}^{-1/4}. \label{eq:proton_gamma_max}
\end{align}
In principle the coefficient $\alpha$ for electrons can differ from that of the protons. In our analysis, the only constraining one is that of the protons.
For the parameters of the model presented here, a large magnetic field is required and therefore the electrons are in the fast cooling regime with $\gamma_{\min} > \gamma_c$. The resulting electron distribution is a broken power-law with index 2 between $\gamma_c$ and $\gamma_{\min}$ and index $(p-1)$ between $\gamma_{\min}$ and $\gamma_{\max}$. On the other hand, the protons are in the slow cooling regime. Therefore their distribution function is a single power-law with index $p$ between $\gamma_{p,  \min}$ and $\gamma_{p, \max}$ (for the parameters we find, $\gamma_{p, \max} < \gamma_{p, c}$, see below).

The observed synchrotron spectra from these particle distributions have the broken-power law shape with three characteristics frequencies \citep{1970RvMP...42..237B, 1999ApJ...511..852G,Granot:2000tta, Granot_2002}. The observed characteristics spectral peak frequency of photons emitted by the particles at $\gamma_{\min}$ is
\begin{equation} \label{equ:vmin}
    \nu_{\min} = \frac{3}{4\pi}\frac{q B}{m_x c}\gamma_{\min}^2 \left[\frac{\Gamma}{1+z}\right] = \left\{ \begin{array}{ll}
    3 \times 10^{13} \hspace{0.3em} f(p)^2 E_{53}^{1/2} t_{day}^{-3/2} \epsilon_{B,-2}^{1/2}\epsilon_{e,-1}^2 \xi_e^{-2} \quad \text{Hz} & ({\rm electrons}), \\
    5 \times 10^{3} \hspace{0.3em} f(p)^2 E_{53}^{1/2} t_{day}^{-3/2} \epsilon_{B,-2}^{1/2}  \epsilon_{p,-1}^2 \xi_p^{-2} \quad \text{Hz} & ({\rm protons}), \\ \end{array} \right.
\end{equation}
(here and below, we add the subscript $x$ which is $e$ for electrons and $p$ for protons).

The cooling frequencies of the electrons and of the protons are
\begin{equation} \label{equ:vc}
    \nu_{c} = \frac{3}{4\pi}\frac{q B}{m_x c}\gamma_{c}^2
        \left[\frac{\Gamma}{1+z}\right] = \left\{ \begin{array}{ll}
     5.3 \times 10^{14} \hspace{0.3em}E_{53}^{-1/2} t_{day}^{-1/2} \epsilon_{B,-2}^{-3/2} n_0^{-1} \quad \text{Hz} & ({\rm electrons}), \\
     1.1 \times 10^{31} \hspace{0.3em}E_{53}^{-1/2} t_{day}^{-1/2} \epsilon_{B,-2}^{-3/2} n_0^{-1}\quad \text{Hz} & ({\rm protons}). \\ \end{array} \right.
\end{equation}
Note that the proton cooling frequency $\nu_{p,c}$ {is} $ \sim 10^{16}$ times larger than the electron cooling frequency $\nu_{e,c}$.

The maximum synchrotron frequencies for electrons and protons are
\begin{equation} \label{equ:vmax}
    \nu_{\max} =\frac{3}{4\pi}\frac{q B}{m_x c}\gamma_{\max}^2
     \left[\frac{\Gamma}{1+z}\right] = \left\{ \begin{array}{ll}
     4 \times 10^{23} \alpha^{-1} E_{53}^{1/8}n_{0}^{-1/8}t_{day}^{-3/8} \quad \text{Hz} & ({\rm electrons}), \\
     7 \times 10^{26} \alpha^{-1} E_{53}^{1/8} n_{0}^{-1/8}t_{day}^{-3/8}\quad \text{Hz} & ({\rm protons}). \\ \end{array} \right.
\end{equation}
Consequently, for fiducial values of the parameters, the highest energy of photons produced by the electron-synchrotron process is limited to few tens of MeV. Therefore, photons having energies much higher than few tens of MeV must have a different origin.
From the previous equation, the cooling and injection frequencies are similar. However, a model in which proton synchrotron explains the TeV  observations requires a large magnetic field, which translates into a fast cooling regime for the electron $\nu_{e,c} \ll \nu_{e,\min}$. However, for protons, $\nu_{p,c} \gg \nu_{p,\min}$ implying that the protons are in the slow cooling regime, and lose their energy inefficiently.

The maximum observed radiative power from a single particle (proton or electron) at observed frequency $\nu_{obs} = \nu_{peak}\Gamma/(1+z)$
is given by
\begin{equation}
P_{\nu_{p,\max}} = (1+z)\frac{2}{9} \frac{m_e^2 c^2 \sigma_T}{q m_p} B\Gamma
\end{equation}
and $P_{\nu_{e,\max}} = ({m_p}/{m_e})P_{\nu_{p,\max}}$ \citep{george1979radiative, 1996ApJ...473..204S, Sari_1998}.
\begin{equation}\label{equ:Fv}
     F_{\nu_{x,\rm peak}} =\frac{N_x P_{\nu_{x,\max}}}{4\pi d_L^2} =\left\{ \begin{array}{ll}
     1 \times 10^{-24}  E_{53} \epsilon_{B,-2}^{1/2} \xi_e n_0^{1/2} d_{L, 28}^{-2} \quad \text{ergs~cm$^{-2}$~s$^{-1}$~Hz$^{-1}$} & ({\rm electrons}),\\
     6 \times 10^{-28} E_{53} \epsilon_{B,-2}^{1/2} \xi_p n_0^{1/2} d_{L, 28}^{-2}\quad \text{ergs~cm$^{-2}$~s$^{-1}$~Hz$^{-1}$} & ({\rm protons}). \\ \end{array} \right.
\end{equation}
where $N_x = (4\pi/3) n_x r^3 $ is the number of particles swept by the blast-wave which are actively radiating, where $n_x = \xi_x n$. Here, $d_{L}$ is the luminosity distance.

For particles accelerated to a power law, $N(\gamma) d\gamma \propto \gamma^{-p}$ above $\gamma_{\min}$ and below $\gamma_{\max}$, the expected photon spectrum thus is a broken power-law shape, with $F_\nu \propto \left\{ \nu^{1/3}, \nu^{-(p-1)/2}, \nu^{-p/2} \right\}$ for $\left\{ \nu < \nu_{\min}, \nu_{\min} < \nu < \nu_c, {\nu_c < \nu < \nu_{\max}} \right\}$ in the slow cooling regime and $F_\nu \propto \left\{ \nu^{1/3}, \nu^{-1/2}, \nu^{-p/2} \right\}$ for $\left\{ \nu < \nu_{c}, \nu_{c} < \nu < \nu_{\min}, {\nu_{\min} < \nu < \nu_{\max}} \right\}$ in the fast cooling regime \citep{george1979radiative, Sari_1998}.

\section{Model constraints derived from the available data of  GRB 190114C} \label{sec:Constraints}

We proceed to interpret the available data of GRB 190114C
within the framework of a hybrid model,
for which the low energy component is explained by synchrotron
radiation from electrons while the high energy TeV component
is required to be proton synchrotron. GRB 190114C was a long-GRB with a prompt energy released $E_{\rm iso} \simeq 2.5 \times 10^{53}$~erg
\citep{2020ApJ...890....9A}. It is seen at the
redshift $z=0.4245$; identified by the
Nordic Optical Telescope \citep{web:nasa:stats} and further established by Gran
Telescopio Canarias \citep{web:nasa:23708}.
Its prompt phase was recorded over an energy band of $8$~keV - $100$~GeV
by the \textit{Swift}-Burst Alert Telescope (BAT) \citep{web:nasa:swift}, the Gamma-ray Burst Monitor
(GBM) \citep{web:nasa:Fermi}, and the Large Area Telescope (LAT) \citep{web:nasa:FermiLAT}. The reported duration is $T_{90} = 25$~s \citep{web:nasa:swift, MAGIC:2019irs}, although the GBM collaboration reported a duration of $T_{90} = 116$~s \citep{web:nasa:Fermi}. This longer duration is explained by the observation of a weak second emission episode after the initial signal. This episode was interpreted as emission from the afterglow \citep{2020ApJ...890....9A}, and in this paper we make the same assumption.

The afterglow  follow-up observations were carried out by many instruments. As a
result a good temporal and multi-wavelength data set
exists during the early afterglow phase when the burst was bright enough in the TeV band to allow the MAGIC instrument to
measure the spectrum in five time bins; 68 - 110 s, 110 - 180 s, 180 - 360
s, 360 - 625 s and 625 - 2400 s \citep{MAGIC:2019irs}. The first two time intervals
have a spectral data in the keV to  GeV band (XRT, GBM, LAT) along with
data in the TeV band from MAGIC.
Inspection of the available data reveals that the afterglow spectra during the first
two time bins centered at $90$~s
and $120$~s are characterized by two peaks, with the
lowest peak at
energy around 10 keV which we attribute to synchrotron emission
from electrons, and a second
peak
between the GeV and the TeV bands which we interpret as synchrotron emission from
protons. In this section, we derive the relations between the model parameters to satisfy these two assumptions.

According to the XRT online repository \citep{2009MNRAS.397.1177E}, the
X-ray spectrum of GRB 190114C is fitted by an absorbed power-law, characterised by a photon index of $1.7^{\pm
0.04}$. With the XRT online tool \citep{2009MNRAS.397.1177E}, we checked that this holds in the first and second time bins independently. We find the photon index to be $1.68^{+0.12}_{-0.11}$ between 68 and 110s, and $1.58^{\pm 0.1}$ between 110 and 180s. This power-law can naturally be produced in the electron-synchrotron model. Since for
the fiducial values of the model parameters discussed herein (see discussion above) electrons are expected to be in the fast cooling regime,
the spectral slope can be obtained in two cases. The first one
corresponds to (i) $\nu_c < \nu_{\rm XRT} < \nu_{\min}$, producing a
spectral index of $-1.5$ comparable to the spectral index found from
the XRT repository. In this case the electron index cannot be
determined and we can take $p_e \sim 2.2$. We discuss the impact of
this assumption below. In the other case, $\nu_c < \nu_{\min} < \nu_{\rm XRT} < \nu_{\rm max}$, the electron index is deduced to be $p_e = 1.2 $ as the observed photon spectral index is 1.6. It also
requires  $\nu_{\rm max}$ to be close to 10~keV to produce the first
hump. However, inspection of Equation \eqref{equ:vmax} reveals that
this case is difficult to achieve as it requires the observed time, $t_{day}$ to be long
(which is in contrast with the observation at hundreds seconds) and a large $\alpha$, meaning
a low particle acceleration efficiency.

The available multi-wavelength data are now used to constrain the free model
parameters. For this purpose, we use the available data  at the two time bins
centered at $90$~s and $120$~s.
\citet{2020ApJ...890....9A}  reported a break at 4.72 keV (at 68 - 110s) and  5.6 keV (at 110 - 180s), and therefore
we assume that the injection frequency
$\nu_{e,\min}$ is equal to 5.5 keV consistent with those findings.  The large magnetic field requirement in our model implies that the cooling frequency $\nu_{e,c}$
to be below the XRT band, \textit{i.e.} smaller than 0.3 keV. Using Equation \eqref{equ:vmin} for the value of the injection frequency yields

\begin{equation}
E_{53}^{1/2}\epsilon_{B,-2}^{1/2} \epsilon_{e,-1}^2 \xi_e^{-2} = \left \{
\begin{aligned} \label{equ: v_min limit}
 &  53.65   & ~~~~~~ &{\rm at~}  t=90\hspace{0.3em}{\rm s},\\
 &  82.60   & & {\rm at~} t=120\hspace{0.3em}{\rm s},
\end{aligned}  \right.
\end{equation}
while using Equation \eqref{equ:vc} for the cooling frequency gives
\begin{equation} \label{equ: v_c limit}
E_{53}^{-1/2}\epsilon_{B,-2}^{-3/2} n_0^{-1} \leq \left \{
\begin{aligned}
& 4.42 & ~~~~~~ &{\rm at~} t=90\hspace{0.3em}{\rm s},\\
& 5.10  & ~~~~~~ &{\rm at~} t=120\hspace{0.3em}{\rm s}.
\end{aligned}  \right.
\end{equation}
Similarly, the maximum observed photon energy, $\sim 0.23$~TeV, can be used to set an estimation of the maximum frequency of photons radiated from the protons, see Equation \eqref{equ:vmax}. The spectrum suggests that $\nu_{\rm max}$ is in fact smaller than this value as the $\nu F_\nu$ spectrum is reported by MAGIC  to decrease. We will assume for the sake of estimating the parameter relations that $\nu_{p,\rm max} \equiv 0.23$ TeV. We get :
\begin{equation} \label{equ: v_maxp limit}
\frac{1}{\alpha} \left ( \frac{E_{53}} {n_0} \right ) ^{1/8} = \left \{
\begin{aligned}
& 6.25 \times 10^{-3} & ~~~~~~ &{\rm at~}  t=90\hspace{0.3em}{\rm s},\\
& 6.96  \times 10^{-3} & ~~~~~~ &{\rm at~}  t=120\hspace{0.3em}{\rm s}.
\end{aligned}    \right.
\end{equation}
Note that when considering the maximum observed photon energy, we took similar values at both 90~s and 120~s.

We now use the specific fluxes at $\nu_{\min}$ and $\nu_{p, \max}$ to constrain the parameters. We first assume a power-law index $p=2.2$ for both electrons and protons,
the fluxes corresponding to the observed energies at the time period of 90 s ($F_{\nu_{obs} =5.5
~\rm keV} \approx 2.56 \times
10^{-26}$ ergs cm$^{-2}$ s$^{-1}$ Hz$^{-1}$ and $F_{\nu_{obs}=0.23~\rm TeV} \approx 9.5 \times
10^{-34}$ ergs cm$^{-2}$ s$^{-1}$ Hz$^{-1}$) and 120 s ($F_{\nu_{obs} = 5.5~\rm keV} \approx 1.41 \times 10^{-26}$ and $F_{\nu_{obs} = 0.23~\rm TeV} \approx4.15 \times 10^{-34}$ ergs cm$^{-2}$ s$^{-1}$ Hz$^{-1}$),
can be explained using the synchrotron
fluxes ($F_{\nu_e}$ and $F_{\nu_p}$) calculated from Equations
\eqref{equ:vmin}, \eqref{equ:vc} and \eqref{equ:Fv}. Further using Equation \eqref{equ: v_min limit} provides
restrictions on the parameter space as follows,
\begin{equation} \label{equ: Fve}
E_{53}^{0.75} \epsilon_{B,-2}^{-0.25} \xi_e \approx \left \{
\begin{aligned}
&0.21 & ~~~~~~ &{\rm for~}  t=90\hspace{0.3em}{\rm s},\\
&0.12 & ~~~~~~ &{\rm for~}   t=120\hspace{0.3em}{\rm s}.
\end{aligned} \right.
\end{equation}
and
\begin{equation} \label{equ: Fvp}
 E_{53}^{1.3} \epsilon_{B,-2}^{0.8} n_0^{0.5} \epsilon_{p,-1}^{1.2}\xi_p^{-0.2} \approx \left \{
\begin{aligned}
& 2.4 \times 10^{5} & ~~~~~~ &{\rm for~}  t=90\hspace{0.3em}{\rm s},\\
& 1.3 \times 10^{5} & ~~~~~~ &{\rm for~}   t=120\hspace{0.3em}{\rm s}.
\end{aligned} \right.
\end{equation}
Equation \eqref{equ: Fve} is derived from the relation $F_{\nu_e} = F_{\nu_e,peak} (\nu_{e,\min}/\nu_{e,c})^{-1/2} (\nu/\nu_{e,\min})^{-p/2}$ for the synchrotron emission of the fast-cooling electrons and Equation \eqref{equ: Fvp} is obtained from $F_{\nu_p} = F_{\nu_p,peak} (\nu/\nu_{p,\min})^{-(p-1)/2}$ relevant for the synchrotron emission from the protons.

Using those two last equations and assuming the value $p = 2.2,$ one can express $\epsilon_B$ and $\xi_e$ as,

\begin{equation} \label{eq:epsilonB}
 \epsilon_{B,-2}
  = \left \{
\begin{aligned}
& 5.3 \times 10^6 ~ \xi_p^{\frac{1}{4}}E_{53}^{-\frac{13}{8}} \epsilon_{p, -1}^{-\frac{3}{2}} n_0^{-\frac{5}{8}}& ~~~~~~ &{\rm for~}  t=90\hspace{0.3em}{\rm s},  \\
& 2.5 \times 10^6 ~ \xi_p^{\frac{1}{4}}E_{53}^{-\frac{13}{8}} \epsilon_{p, -1}^{-\frac{3}{2}} n_0^{-\frac{5}{8}}& ~~~~~~ &{\rm for~}  t=120\hspace{0.3em}{\rm s},
\end{aligned} \right.
\end{equation}
and for $\xi_e$:
\begin{equation} \label{eq:xi_e}
\xi_e =
  \left \{
\begin{aligned}
  &10.08 ~ \xi_p^{\frac{1}{16}} E_{53}^{-\frac{37}{32}} \epsilon_{p, -1}^{-\frac{3}{8}} n_0^{-\frac{5}{32}}& ~~~~~~ &{\rm for~}  t=90\hspace{0.3em}{\rm s},\\
  &4.76 ~ \xi_p^{\frac{1}{16}} E_{53}^{-\frac{37}{32}} \epsilon_{p,-1}^{-\frac{3}{8}} n_0^{-\frac{5}{32}}& ~~~~~~ &{\rm for~}  t=120\hspace{0.3em}{\rm s}.
\end{aligned} \right.
\end{equation}
We further use Equations \eqref{eq:epsilonB} and \eqref{eq:xi_e} with Equation \eqref{equ: v_min limit} to obtain the value of $\epsilon_e$ :
\begin{equation}\label{eq:epsilone}
    \epsilon_{e, -1} =
    \left \{
\begin{aligned}
& 1.54 ~ E_{53}^{-1}
& ~~~~~~ &{\rm for~}  t=90\hspace{0.3em}{\rm s},  \\
& 1.09 ~ E_{53}^{-1}
& ~~~~~~ &{\rm for~}  t=120\hspace{0.3em}{\rm s},
    \end{aligned} \right.
\end{equation}

In our analysis below, we also use the physical condition that the electrons, protons and magnetic field energy are all obtained from the post-shock thermal energy, namely $\epsilon_B+\epsilon_e+\epsilon_p\leq 1$. Inspections of those equations reveals that the time bin centered at 90 s gives the adequate parameter magnitudes. Therefore, hereinafter we will use
the constraints obtain from this time bin only and we will show that the time evolution can be well reproduced by the afterglow dynamics alone.

It is clear that satisfying the condition  $\epsilon_B < 1$ requires a combination of (i) a small $\xi_p$  \textit{i.e.} a small fraction of protons being accelerated into a power-law, (ii) a large total kinetic energy $E_{53}$, (iii) a large circumburst medium density $n_0$ and (iv) a large (but smaller than the unity) $\epsilon_p$, \textit{i.e.} a large fraction of energy given to accelerated protons. We note that this combination of parameters also leads to a small value of $\xi_e$, namely, only a small fraction of the electron population is accelerated to a power law.

We show the constraints on the parameters in Figure \ref{fig:epsilon_B plot}, which displays the value of $\epsilon_B$ as a function of $n$, $E_{53}$ and $\alpha$. Inspection of Equation \eqref{equ: v_min limit} reveals that the dependence of $\epsilon_e$ on $E_{53}, n_0$~and $\epsilon_p$, implies that the condition $\epsilon_e < 1$ is automatically satisfied, when satisfying the limitation given by Equation \eqref{eq:epsilonB}. Finally, inspecting Equation \eqref{equ: v_maxp limit}, it comes that the dependence on $E_{53}$ and $n_0$ requires $\alpha$ to not be too large, with $\alpha \sim 100$ owning for the weak dependence on $E_{53}$ and $n_0$.

\begin{figure}[t!]
     \centering
     \begin{tabular}{ccc}
        \includegraphics[width=0.31\textwidth]{Figures/epsilon_B_for_xi_p_=_0.01.pdf} &
        \includegraphics[width=0.31\textwidth]{Figures/epsilon_B_for_xi_p_=_0.01.pdf} &
        \includegraphics[width=0.31\textwidth]{Figures/epsilon_B_for_xi_p_=_1.pdf} \\
        $\xi_p = 0.01$ & $\xi_p = 0.1$ & $\xi_p = 1$
     \end{tabular}
        \caption{ Colors give the value of the magnetization parameter $\epsilon_B$ (in logarithmic scale) as a function of the ambient density, $n_0$ and total energy, $E_{53}$ from Equation \eqref{eq:epsilonB} for the time bin centered at 90~s. In making the plot, we assumed $\epsilon_p = 0.8$ and display the results for three different values of $\xi_p$, namely, $\xi_p = 10^{-2},10^{-1},1$, respectively presented from left to right. The thin black lines represent combination of parameters that result in $\epsilon_B = 10^{-1}, 10^{-2}, 10^{-3}$.}
        \label{fig:epsilon_B plot}
\end{figure}

\section{Results : Explaining the TeV observation of GRB 190114C with proton synchrotron}\label{sec:RESULTS}

In this section, we present several solutions for a proton synchrotron model
to explain the TeV emission of GRB 190114C. In particular we consider three assumptions
based on the uncertainty of the particle acceleration process: (i) similar particle
spectral indices, $p_e = p_p$ and similar fraction of accelerated particles, $\xi_e
= \xi_p$ \textit{i.e} the acceleration process is similar for protons and electrons
in terms of number of accelerated particles and the obtained spectral shape. (ii)
Similar particle spectral indices, $p_e = p_p$ but $\xi_e \neq \xi_p$, namely the
acceleration process accelerates different proportions of electrons and protons, but
produces a similar spectral shape. And (iii) different spectral indices, $p_e \neq
p_p$ with $\xi_e = \xi_p$, meaning that the acceleration process produces a different spectral shape. We find that the effect of $\xi_p$ on the resulting
spectra is not significant, therefore we did not have to assume an extra degree of freedom (see below)\footnote{Since we normalize the flux to the observed flux,
reducing $\xi_p$ enforces an increase in the values of other free model parameters,
such as the magnetization.}.

The results are presented in Figures \ref{fig:xie = xip graphs}, \ref{fig:xie neq xip graphs} and \ref{fig:two graphs} respectively. The data are extracted from \cite{MAGIC:2019irs}
and are presented here for convenience. In particular, we did not aim at producing
the statistical best fit to the data, but only to demonstrate the ability of our
model to reproduce the fluxes and characteristic breaks in the XRT and MAGIC bands.
This allows us to constrain the values of the free parameters of our model.

\subsection{Case (i): \texorpdfstring{$\xi_e = \xi_p$}{TEXT}} \label{sec:case1}
We first examine the  assumption of similarity between the injection fraction of
electrons and protons in the acceleration process.
We impose $\xi_e =
\xi_p$ in Equation \eqref{eq:xi_e} for $t =90$~s which results in
\begin{align}  \label{eq:xie=xip}
    \xi_e \equiv 11.76 \quad E_{53}^{-\frac{37}{30}}\epsilon_{p,-1}^{-\frac{6}{15}} n_0^{-\frac{1}{6}}.
\end{align}
Also in order for $\epsilon_B$ to be smaller than unity, the kinetic energy is
required to be greater than the observed released prompt energy
$E_{\rm iso}$ as seen explicitly from Equation \eqref{eq:epsilonB}.
Here we choose $E = 4 \times 10^{54}$~ergs, $\epsilon_p = 0.8$, and ambient density $n = 80$~cm$^{-3}$, then Equation \eqref{eq:xie=xip} yields $\xi_e =\xi_p = 0.026$, for which $\epsilon_B = 0.152$ and $\epsilon_e = 3.85 \times 10^{-3}$ and so the relation $\epsilon_B+\epsilon_e+\epsilon_p\leq 1$
is satisfied. For this set of parameters, only $2.6\%$ of the electrons and
protons injected attain the power-law distributions behind the shock front and the rest assume a thermal distribution with temperature lower than $\gamma_{\min}$. We do
not attempt to model the radiation from these thermal particles, as their contribution is below the observed band.
Equation \eqref{equ: v_maxp limit} constrain the numerical coefficient that determines the acceleration efficiency to be $\alpha = 147$.
Above $\gamma_{p, \max}$ we assume an exponential cutoff representing the inability of the acceleration process to accelerate protons to energies above $\gamma_{p, \max}$ (see dashed line in Figure \ref{fig:xie = xip graphs}).
For this set of parameters, the spectra obtained in the time bins
centered at 90 s and 120 s are displayed in Figure \ref{fig:xie = xip graphs}.
This set of parameters results in a spectra and flux consistent with the observed data in both time intervals. The large value of $E$ implies that the
efficiency $\eta=E_{\rm iso}/(E+E_{\rm iso})$ of the prompt emission $\eta \simeq 5.9\%$ is low, but not extremely low. Similarly, the required energy is large, but acceptable. We comment on the low efficiency in the discussion below in section \ref{sec:discussion}.

\begin{figure}[t!]
     \centering
     \begin{tabular}{cc}
        \includegraphics[width=0.45\textwidth]{Figures/SED_at_90s_xie_=_xip.pdf} &
        \includegraphics[width=0.45\textwidth]{Figures/SED_at_120s_xie_=_xip.pdf} \\
        a) $t=90$~s &  b) $t=120$~s
     \end{tabular}
     \caption{Resulting spectra for parameters in our case (i) : $\xi_e = \xi_p $ plotted on top of the data of GRB 190114C at 90~s and 120~s. The SED obtained from
     our model for the time bins (a) [$\rm T_0 + 68- T_0 + 110$] and (b) [$\rm T_0 +
     110- T_0 + 180$]. Continuous line: electron synchrotron, dashed line: proton
     synchrotron and dotted-dashed line: inverse Compton. The data are from \cite{MAGIC:2019irs}
     and are displayed to guide the eye. The parameters are $E_{53} = 40$, $n_0 = 80$, $\epsilon_{p, -1} = 8$, $p = 2.2$ and $\alpha = 147$, resulting in $\epsilon_{B,-2} = 15.2$, $\xi_e = \xi_p = 0.026$ and $\epsilon_{e,-1} = 0.0385$. The cut-off at $\gamma_{p,\max}$ is attained by $\exp (-\nu/(3\nu_{\max}))$.}
     \label{fig:xie = xip graphs}
\end{figure}

\subsection{\texorpdfstring{Case (ii): $\xi_e \neq \xi_p$}{TEXT}} \label{sec:case2}
We next examine a model in which the electrons and protons are
accelerated in different numbers, $\xi_e \neq \xi_p$. The parameters we consider, $E_{53} = 50$,
$n_0 = 130$, $\epsilon_p = 0.8$, $\alpha = 142$, and $\xi_p = 0.2$ provides
adequate spectra that are consistent with the observed data at both 90~s and 120~s. Using Equation \eqref{eq:epsilonB} we obtain
$\epsilon_B = 0.13$, and from Equation \eqref{eq:xi_e} one finds $\xi_e = 0.02$ and $\epsilon_e =  3.08 \times 10^{-3}$.

Our results are presented in Figure
\ref{fig:xie neq xip graphs} for the time intervals centered at 90 s and 120 s.
The parameters of this solution are nearly the same as the ones obtained
for the case $\xi_e = \xi_p$. This is because the protons dominate the
energetic requirements. In particular, it means that the total kinetic energy
of the blast-wave has to be large, and correspondingly, the efficiency is low. In this case it is $\eta \sim 4.8\%$. Owning for the large energy, the interstellar
medium density is also large. Again we point out that all the values we obtain are in range with estimates and uncertainties of GRB energetics and ambient densities.

\begin{figure}[t!]
     \centering
     \begin{tabular}{cc}
         \includegraphics[width=0.45\textwidth]{Figures/SED_at_90s_xie_neq_xip.pdf} &
         \includegraphics[width=0.45\textwidth]{Figures/SED_at_120s_xie_neq_xip.pdf} \\
        a) $t=90$~s & b) $t=120$~s \\
     \end{tabular}
     \caption{Case (ii). Same as Figure \ref{fig:xie = xip graphs}. The parameters are $E_{53} = 50$, $n_0 = 130$, $\epsilon_{p, -1} = 8$, $p = 2.2$, $\alpha = 142$ and $\xi_p = 0.2$, $\epsilon_{B,-2} = 13$, $\epsilon_{e,-1} = 0.0308$, $\xi_e = 0.02$ and so $\xi_e/\xi_p \simeq 0.1$. The cut-off at $\gamma_{p,max}$ is attained by $\exp (-\nu/(3\nu_{\max}))$.}
     \label{fig:xie neq xip graphs}
\end{figure}

\subsection{Case (iii) \texorpdfstring{$\xi_e \equiv \xi_p$}{TEXT} and \texorpdfstring{$p_e \neq p_p$}{TEXT}} \label{best fit}

The critical parameter which determines the  energetic budget is
the proton injection index $p_p$. We have assumed the electron index to be 2.2 in agreement with theory \citep{Sari_1998} and observations \citep{2009MNRAS.397.1177E}. The proton index was then assumed to be equal to the electron index. Here we relax
this assumption. Since the TeV flux is modeled by emission from protons at $\gamma_{p, \max}$, a high value of $p_p$ would require higher energy budget, and lower efficiency. We therefore consider $p_p = 2.1 < p_e$ in order to decrease the energy
budget and increase the efficiency. The solution we present here is given
by the parameters $n = 15$~cm$^{-3}$, $E = 3 \times 10^{54}$~ergs, $\epsilon_p = 0.8$. This gives $\alpha = 35$,
$\epsilon_e = 4.7 \times 10^{-3}$, $\epsilon_B = 0.12$, and $\xi_e = \xi_p = 0.023$. The resulting
spectrum in the two time bins centered at 90 s and 120 s is shown in
Figure \ref{fig:two graphs}.

A consequence of the assumptions used here is that the prompt phase efficiency is increased to $7.7\%$, which is higher than in the two previous cases. Furthermore, the density of the interstellar medium is substantially lower than the value obtained in our previous scenarios.

\begin{figure}[t!]
     \centering
     \begin{tabular}{cc}
        \includegraphics[width=0.45\textwidth]{Figures/SED_at_90s.pdf} &
        \includegraphics[width=0.45\textwidth]{Figures/SED_at_120s.pdf} \\
        a) $t=90$~s & b) $t=120$~s
     \end{tabular}
     \caption{Case (iii) : $\xi_e = \xi_p = 0.023$, and $p_e = 2.2$ while $p_p = 2.1$. Same as Figure \ref{fig:xie = xip graphs}. The parameters we used here are $E_{53} = 30$, $n_0 = 15$, $\epsilon_{p,-1} = 8$, and $ \alpha = 35 $, $\epsilon_{B,-2} = 12$ and $\epsilon_{e,-1} = 0.047$. The cut-off at $\gamma_{p,max}$ is attained by $\exp (-\nu/\nu_{\max})$.}
     \label{fig:two graphs}
\end{figure}

In addition to the data in the two time bins [$\rm T_0 + 68- T_0 + 110$] and [$\rm T_0 +110- T_0 + 180$], MAGIC data exists for three other time bins, namely
$180-360$~s, $360-625$~s and $625-2400$~s \cite{MAGIC:2019irs}. As a consistency
check of our model, we examined the evolution of the proton-synchrotron component
with time, using the same physical parameters, and the self-similar solution to
obtain the dynamical evolution of the Lorentz factor, energy density and magnetic
field, from which the evolution of the flux and characteristic frequencies are
readily obtained. The resulting synchrotron emission from the proton component at
later times is presented in Figure \ref{fig:three graphs} alongside the MAGIC data.
The similarity between the model results and the data is another independent indication for the ability of our hybrid model to provide acceptable explanations to the data at different times.

\begin{figure}[t!]
     \centering
     \begin{tabular}{ccc}
        \includegraphics[width=0.31\textwidth]{Figures/fig_3.pdf} &
        \includegraphics[width=0.31\textwidth]{Figures/fig_4.pdf} &
        \includegraphics[width=0.31\textwidth]{Figures/fig_5.pdf} \\
        $t=260$~s &    $t=520$~s & $t=2200$~s
     \end{tabular}
        \caption{The proton synchrotron component for our model in case (iii), together with the data from the three additional MAGIC observations \citep{MAGIC:2019irs}. The depletion of the observed flux is in agreement with our model.}
        \label{fig:three graphs}
\end{figure}

\section{Additional radiative processes}
\label{sec:5}
\subsection{Inverse Compton Scattering}

The high energy emission of several GRBs  have been interpreted as originating from the synchrotron
self-Compton process \citep{Dermer_2000, 2001ApJ...548..787S,
2009ApJ...703..675N, liu2013interpretation, 10.1093/mnras/stw1175,
fraija2019synchrotron, Derishev:2021ivd}. This process considers the
inverse Compton (hereinafter IC) interaction of synchrotron photons with the electrons that emitted them. Indeed, it is possible to interpret the observed TeV data using IC rather than proton synchrotron, if one assumes different values of the free model parameters, as follows.

The break frequencies for the IC component can be obtained within the Thompson
regime, since Klein-Nishina corrections are important only at early times, typically $t^{obs}\lesssim 60$~s, \textit{i.e.} earlier than the time of the first MAGIC observation, see \citet{2001ApJ...559..110Z}. The $\nu F_{\nu}$ spectral
peak is achieved at frequency $\nu_{IC,\min}= 2 \gamma_{e,\min}^2 \nu_{e,\min} $ or
$\nu_{IC,c}= 2 \gamma_{e,c}^2 \nu_{e,c}$ in the slow and fast cooling
regimes respectively. For the asserted model specifications (\ref{best fit}), these frequencies are defined as
\begin{align}
    \nu_{IC,\min} &= 1.75 \times 10^{20}\hspace{0.3em}f(p)^2 E_{53}^{3/4} n_{0}^{-1/4} t_{day}^{-9/4} \epsilon_{B,-2}^{1/2} \xi_e^{-4} \epsilon_{e,-1}^4 \text{~Hz},
\end{align}
and
\begin{align}
    \nu_{IC,c} &= 5.52 \times 10^{22}\hspace{0.3em}E_{53}^{-5/4} n_{0}^{-9/4} t_{day}^{-1/4} \epsilon_{B,-2}^{-7/2} \text{~Hz} .
\end{align}
The peak-flux of the IC spectrum is estimated as \citep{2001ApJ...548..787S}
\begin{align}
F_{IC,\nu_{\rm peak}} = \frac{1}{3}\sigma_T n_e r F_{e, \nu_{\rm peak}} \approx 1.5 \times 10^{-31}\hspace{0.3em} \xi_e^2 E_{53}^{5/4} n_0^{5/4} t_{day}^{1/4} D_{28}^{-2} \epsilon_{B,-2}^{1/2} \text{~erg~cm$^{-2}$~s$^{-1}$~Hz$^{-1}$}.
\end{align}
The specific flux of the IC spectrum at the time interval $t = 90$~s when compared to the observed flux at $0.23$~TeV (MAGIC observed band) gives the constraints for the parameters of this model as
\begin{align}
\xi_e^2 E_{53}^{\frac{5}{8}} n_0^{\frac{1}{8}} \epsilon_{B,-2}^{-\frac{5}{4}} = 0.236.
\end{align}
From the above relation we can get
\begin{align}
\epsilon_{B,-2} = 3.17~E_{53}^{\frac{1}{2}} n_0^{\frac{1}{10}} \xi_e^{\frac{8}{5}}.
\end{align}
Thus, for $\xi_e \sim 10^{-2}$, the data can be explained if $\epsilon_B \approx 10^{-3}$.

However, within the framework of the proton-synchrotron model, for the parameters we considered to explain the data of GRB 190114C in the section \ref{best fit}, the energy peak
flux of the IC spectrum is around $\sim 40$~TeV (90 s) and $\sim 21$~TeV (120 s)  with
the specific flux of $\sim 8.3 \times 10^{-10}$~ergs~cm$^{-2}$~s (90~s) and  $\sim6.2\times 10^{-10}$~ergs~cm$^{-2}$~s (120~s), respectively. The expected contribution from the inverse Compton scattering process is shown
in Figures \ref{fig:xie = xip graphs}, \ref{fig:xie neq xip graphs} and \ref{fig:two graphs}, corresponding to the three models we presented. These figures
show that the IC component is always subdominant with respect to the proton
synchrotron component in the TeV band. This was expected since our model requires $\epsilon_{B} \gg \epsilon_e$, resulting in a strong suppression of the SSC component. Thus it states that our model presents an alternative
to the IC component models, as suggested by \cite{Derishev:2021ivd}.

\subsection{Photo-pion Production}

Another mechanism to produce photons
with energy around 1 TeV is photon-pion interactions and their subsequent
electro-magnetic pair cascade. Pions ($\pi^0\hspace{0.3em}\&
\hspace{0.3em}\pi^+$) originated from the collisions of the high energy
protons and the seed synchrotron photon-field emitted from the electrons,
produce $\gamma-$ray photons accompanied by neutrinos,
$p^+ + \gamma \rightarrow \pi^0 + p^+$ and $p^+ + \gamma \rightarrow
\pi^+ + n$. Neutral pions quickly decay to photons, $\pi^0 \rightarrow
\gamma + \gamma$, and charged pions decay to positrons
and neutrinos in the decay chain ($p^+ \rightarrow \mu^+ + \nu_\mu$, $\mu^+
\rightarrow e^+ +\nu_e + \Bar{\nu}_\mu$). High energy electrons and
positrons ($e^+$) generated from the muon ($\pi^\pm$) decay
as well as muons and charged pions can emit a
substantial amount of synchrotron radiation. We seek here to compute the
relevance of this process.

In appendix \ref{sec:appendix_ppi}, we derive the expression for the cooling time
by photohadronic interaction, given by Equation \eqref{eq:t_ppi}. This time is
compared to the cooling time by synchrotron given by Equation
\eqref{eq:proton_cooling_synch}:
\begin{align}
    \left. \frac{t_{synch}}{t_{p \pi}} \right |_{\gamma_{p,\max}} = 2.9 \times 10^{-8} \left( \frac{p-2}{p-1} \right ) E_{53}^{\frac{45}{32}} \epsilon_p^{\frac{15}{8}} n_0^{\frac{37}{32}} t_{\rm day}^{\frac{9}{8}} \xi_p^{-\frac{5}{16}}
\end{align}
were we also use the constraints on $\epsilon_B$, $\epsilon_e$ and $\xi_e$ given by Equations \eqref{eq:epsilonB}, \eqref{eq:xi_e} and \eqref{eq:epsilone}.
Therefore, for the protons, the synchrotron cooling rate is orders of magnitude faster than the photohadronic interaction, and most
of the energy radiated by the protons is done via synchrotron radiation.  In other words,
the contribution from the hadronic cascade in our scenario can be neglected and no modification
of the TeV component is expected. This also means that, for our model, the neutrinos fluence at PeV energy is expected to be small, challenging observations and constraints by IceCube and future instruments.

\section{Discussion} \label{sec:discussion}
\subsection{Radiative efficiency}

As we show here, our proton synchrotron model is able to explain the VHE emission of GRB 190114C. Within the framework of this model, we find that the prompt phase signal has an energy conversion efficiency of a few percent. Indeed,
our most favorable model with $\xi_e = \xi_p$ and $p_e \neq p_p$ requires
a blast-wave with energy $3\times 10^{54}$ ~ergs, resulting in a prompt efficiency
around eight percents. The other models we considered all have a few percent efficiency as well.
Interestingly, this efficiency is comparable to the radiative efficiency predicted in the
internal shock scenario used by many authors to explain the prompt phase; see e.g. \cite{KPS97, PSM99, GSW01}.

At first, our model requirements on the efficiency seem inconsistent with
previous findings, specifically for the burst observed by Fermi LAT for which a high
(around 50\%) prompt radiative efficiency is usually determined, \textit{e.g.}
\cite{CFH11}. This trend seems to also be retrieved for bursts observed by the
Neils Gherel Observatory \citep{CFH10}, as well as older bursts \citep{YHS03}. However, we note that:
\begin{itemize}
    \item [(i)] the analysis presented in those papers heavily relies on the
numerical coefficients chosen for the emission  process. Using updated coefficients,
\cite{FP06} re-evaluated the efficiency of several bursts, finding it to be lower,
in the order of few percents, so of similar magnitude than the requirements from our model.
    \item [(ii)] their analysis also further relies on assuming that all the electrons participate in the radiative process, namely $\xi_e =1$, resulting in an
under-evaluation of the kinetic energy of the blast-wave compared to the case $\xi_e < 1$ as seen from rearranging Equation \eqref{eq:xi_e},
\begin{equation}
    E_{53} \simeq 7.38 \quad \xi_p^{\frac{2}{37}} \epsilon_{p,-1}^{-\frac{12}{37}} n_0^{-\frac{5}{37}} \xi_e^{-\frac{32}{37}},
\end{equation}
for the parameters that satisfy the relation \eqref{eq:epsilonB} at $t = 90$~s.
Note that our proton synchrotron model forbids $\xi_e = 1$, as it would necessarily results in $\epsilon_B > 1$. Instead, our model requires $\xi_e$ to be in the order of $10^{-2}$.
Recent studies analysed the spectral effect of a thermal population resulting in
an incomplete acceleration of particles, see \textit{e.g.} \citet{WBI18, WDB22}.
Yet, it remains to understand how the modification to the spectral energy distribution impacts the recovery of the blast-wave parameters. Thus, one can conclude that currently there is still no reliable measurement of a high efficiency during GRB prompt emission.
\end{itemize}

\subsection{Constraints on the electron injection fraction.}

The strong constraint on the magnetic field equipartition parameter $\epsilon_B$, requires that the number of electrons participating in the radiation process is smaller than unity, of the order of $10^{-2}$, weakly sensitive to all parameters but the kinetic energy of the blast-wave, $E_{53}$. The 5.5 keV break and the low-energy slope below it, suggest that the thermal component made of the bulk of the electrons, should have a temperature much smaller than $\gamma_{\min}$. If this was not the case, the low-energy slope below $\nu_{\min}$ would be different and entailed to the exact injection function.

In other words, we require an acceleration scenario in which the thermal component and accelerated particles are two clearly separate entities. Such an injection
function was originally proposed by \citet{2005ApJ...627..861E}. In that paper, they
also studied the effect of this assumption on the emission properties of GRB
afterglow. In particular they found that such model can be constrained by early
afterglow emission, underlying the necessity of multi-wavelength observations.
Multi-wavelength observations of GRB 190114C afterglow, which include the temporal evolution in the radio and optical bands, was reported by \citet{2021MNRAS.504.5685M}. And indeed, the constraint upon the parameter $\epsilon_e$ at time interval of 65~s reported in that work, resulted with an electron injection fraction of 2$\%$, which is comparable to the fraction $\xi_e$ obtained in our model. These consistent results therefore serve as an independent support to our model.

\section{Conclusion}

\label{sec:Conclusion}

The broken power-law spectrum of the electron-synchrotron emission  model is a
standard prediction of GRB afterglow theories, which was successful in explaining the spectra at energy lower than a few MeV. Yet, the source of the VHE afterglow component is still being
investigated. GRB 190114C is one such a burst with a high energy ($0.2 \sim 1$~TeV) peaking
afterglow spectral component as reported by MAGIC within the epoch $60$~s to $2400$~s. We argue here that the source of the VHE segment
is of proton synchrotron origin, while the sub-MeV component is explained by electron synchrotron radiation within the framework of the classical fireball evolution scenario. In our model, the protons are co-accelerated with the electrons in the blast wave such that they can radiate in the (sub-)TeV band.

We provided the constraints that this model must satisfy under different
conditions for the parameters describing the uncertainty of the particle acceleration process. We presented three models
with different injection parameters characterizing the particle acceleration and injection : (i)
our first model has similar fractions of electrons and protons accelerated to a high-energy power law, $\xi_e = \xi_p = 2.6 \times 10^{-2}$; (ii) the second one has $\xi_e \neq \xi_p $ and (iii) the last one has $\xi_e = \xi_p = 2.3\times 10^{-2}$, but (slightly) different power law indices, $p_e = 2.2$ and $p_p = 2.1$, that is to say electrons and protons are described by a
different injection index. All the scenarios in our model exhibit external medium density $n \sim 10 - 100$~cm$^{-3}$, which is somewhat higher than the fiducial value often assumed in explaining many GRB afterglows, $n \approx 1$~cm$^{-3}$, but not unreasonably high. In fact a similar value was inferred from fitting the much later time radio and optical afterglow data  \citep{2021MNRAS.504.5685M}.  Similarly, a reasonably high blast-wave
kinetic energy ($E \approx 3\times10^{54}$~erg) which, given the inferred energy from the observed photons of $\approx 2.5 \times 10^{53}$~erg during the prompt phase, leads to a considerable prompt
radiative efficiency of $\approx 8\%$.  We point out that this value of efficiency is similar to the efficiency expected in kinetic energy conversion by internal shocks. Finally, we note the degeneracy in determining the values of the explosion energy and the ambient density: for a given observed flux, a low ambient density can be compensated by a higher explosion energy for a fixed magnetization, $\epsilon_B$, thereby reducing the radiative efficiency of the prompt phase.

\citet{2021MNRAS.504.5685M} analysed and interpreted the radio, optical and X-ray data
of GRB 190114C at both early and later times. Very interestingly, they concluded that the kinetic energy $E$
was at least one order of magnitude higher than the observed isotropic
equivalent energy. This result is in excellent agreement with
the constraints on the energy we obtained here. Further examination of our model requirements conveyed low
acceleration
efficiency with $\alpha > 10$. In addition, we find that in order to reproduced the TeV band data via the proton synchrotron process, the fractions of both electrons and protons accelerated to a high energy power law, $\xi_e$ and $\xi_p$ are at the order of few percent. It also requires a large magnetic field, hence a large kinetic energy, from which, using the observed data at the {\it Swift-XRT} band, a constraint on the fraction of electrons accelerated, $\xi_e$ is imposed as well.

 The injection fractions we derive, of a few percents, are similar to the ones derived using the radio, optical and X-ray observational data in \citet{2021MNRAS.504.5685M}. Furthermore, these values are consistent with
the theoretical values found in particle-in-cell simulations of particle acceleration in relativistic shock waves \citep{2008ApJ...682L...5S}. Using the same values, we showed in Figure \ref{fig:three graphs} that the proton-synchrotron mechanism is in agreement with the MAGIC data at later times as well, until the last existing MAGIC observation at 2400~s.  Importantly, we showed that for the constraints set by the data, the proton synchrotron emission is dominant at the TeV band over other emission processes such as SSC and photo-pion production mechanism. This makes this model a viable alternative to leptonic models.

To conclude, we presented a proton synchrotron model to
explain the VHE afterglow spectrum of GRB 190114C and found that it can accommodate the
available observed MAGIC data-set from the epoch of 68 s up until 2400~s. The parameters we found are consistent with both independent measurements, as well as theoretical predictions of particle acceleration. Furthermore, the uncertain values of the explosion energy and ambient density, are within one-two orders of magnitude from 'fiducial' values often assumed in the literature. These results, therefore, point to the need of a more thorough investigation of the role of proton-synchrotron emission during the afterglow phase in GRBs.

The VHE band limit of the GRB afterglows is anticipated to
reach higher magnitudes with CTA in the near future. Therefore, highlighting the role
played by highly relativistic protons in the GRB afterglow theory is of crucial
importance. Conversely, the acceleration mechanism and radiative efficiency of protons are expected to
be further constrained by the observation of very-high energy emission ($\sim 1$~TeV)
in the afterglow spectrum.