\section{Introduction}
The spectrum of hadrons provides a unique connection between QCD and experiment. Many experiments in hadron spectroscopy have utilized elastic and charge-exchange $\pi N$ scattering,
but certain resonances are difficult to discern from this production mechanism, which leads to the ``missing resonance problem'' in the spectrum of baryons \cite{Koniuk:1979vw}. Production processes other than $\pi N$ scattering allow us to better probe the spectrum of hadrons and increase our understanding of the underlying physics \cite{Ronchen:2015vfa}. One possible alternative production process that is being utilized at modern-day facilities like JLab, ELSA, and MAMI is photoproduction. There, a photon is absorbed on a stable hadron, like a proton, which is then excited into a resonance that decays to various lighter hadrons. Recent theoretical developments \cite{Briceno:2014uqa,Briceno:2015csa} now allow lattice QCD calculations of the photoproduction of two-body resonances. However, as the baryon sector is very challenging, it is helpful to first investigate the photoproduction of the simplest known resonance, the $\rho$ meson. Here, we give an overview of our lattice QCD calculation of this process \cite{Alexandrou:2018jbt}.

\section{Description of the photoproduction}
Because the $\rho$ is a resonance and thus not stable under the strong interaction, the photoproduction process $\pi\gamma\to \rho$ is obtained by analytical continuation of the more general process $\pi\gamma\to\pi\pi$ to the $\rho$ pole. Here the $\pi\pi$ state is in $P$-wave, has isospin $I=1$ and $J^{PC}=1^{--}$. The process $\pi\gamma\to\pi\pi$ is described by the infinite-volume matrix element $\langle \pi \pi |J_\mu |\pi\rangle$, which has the Lorentz decomposition
\begin{align}
\label{eq:LLdecomp}
\langle \pi\pi | J^\mu(0) | \pi \rangle = \frac{2\mathrm{i} {\cal V}_{\pi\gamma \to \pi\pi}(q^2,s)}{m_\pi} \epsilon^{\nu\mu\alpha\beta} \epsilon_{\nu}(P,m) (p_\pi)_\alpha P_\beta.
\end{align}
Above, ${\cal V}_{\pi\gamma \to \pi\pi}$ is the transition amplitude, which depends on the momentum transfer $q^2 = (p_\pi - P)^2$, $P$ is the $\pi\pi$-state four-momentum vector, $p_\pi$ is the $\pi$-state four-momentum vector and $\epsilon_{\nu}(P,m)$ is the $\pi\pi$-state polarization vector. Further details on the states, their normalization and specific implementations can be found in Ref.~\cite{Alexandrou:2018jbt}. Near the $\rho$ pole, $s_P=m_R^2 + \mathrm{i} m_R \Gamma_R$, the transition amplitude ${\cal V}_{\pi\gamma \to \pi\pi}$ takes the form
\begin{align}
\label{eq:ampV}
{\cal V}_{\pi\gamma \to \pi\pi}(0, s) \sim \frac{G_{\rho\pi\pi} \:G_{\rho\pi\gamma}}{s_P - s},
\end{align}
where $G_{\rho\pi\gamma}$ is the $\rho$ resonance photocoupling and $G_{\rho\pi\pi}$ is the coupling of the $\rho$ to the $\pi\pi$ channel. Making use of Watson's theorem, the transition amplitude can be reformulated as
\begin{align}
\label{eq:V_BW_d}
{\cal V}_{\pi\gamma \to \pi\pi}(q^2,s) & =  \sqrt{\frac{16 \pi}{k\Gamma(s)}} \frac{F(q^2,s) }{\cot{\delta(s) - \mathrm{i}}},
\end{align}
where $\delta$ is the phase shift describing the strong $\pi\pi$ scattering. We investigate the two Breit-Wigner forms described in Ref.~\cite{Alexandrou:2017mpi}, where BWI is the classic $P$-wave Breit-Wigner and BWII is modified by a Blatt-Weisskopf barrier factor. For complete generality, the form factor $F(q^2,s)$ still carries an $s$ dependence, but it no longer has a pole in the $s$ plane within the energy region of interest, below the $K\bar{K}$ threshold.

\section{Finite-volume matrix elements with lattice QCD}
We use a single gauge-field ensemble with $N_f=2+1$ clover-Wilson fermions on a $N_s^3 \times N_t = 32^3\times96$ lattice. The lattice spacing determined from the $\Upsilon(1S)$-$\Upsilon(2S)$ splitting with NRQCD is $a=0.11403(77)$ fm. Our light-quark masses correspond to a pion mass of $am_\pi=0.18295(36)$. Our calculation of the $\rho$ resonance parameters using the L\"uscher method in several moving frames and irreducible representations can be found in Ref.~\cite{Alexandrou:2017mpi}.
To determine the finite-volume matrix elements $\langle \pi,\vec{p}_{\pi} | J_\mu(0,\vec{q})|n,\vec{P},\,\Lambda,\,r\rangle_{FV}$, we calculate the three-point functions
\begin{align}
\label{eq:C3_decomp}
& C_{3,\;\mu,i}^{\vec{p}_{\pi},\,\vec{P},\,\Lambda,\,r}(t_{\pi},t_{J},t_{\pi\pi}) = \langle  O_{\pi}^{\vec{p}_{\pi}}(t_{\pi})  J_\mu(t_{J},\vec{q}) O_{i}^{n,\vec{P},\,\Lambda,\,r}(t_{\pi\pi}) \rangle \cr
&= \sum_{n} Z_\pi^{\vec{p}_{\pi}}\: Z_i^{n,\,\vec{P},\,\Lambda\,\dag}\: \langle \pi,\vec{p}_{\pi} | J_\mu(0,\vec{q})|n,\vec{P},\,\Lambda,\,r\rangle_{FV}
 \times \frac{e^{-E_\pi^{\vec{p}_{\pi}} (t_{\pi} - t_{J})} e^{-E_n^{\vec{P}, \Lambda} (t_{J} - t_{\pi\pi})} }{2E_n^{\vec{P}, \Lambda } 2E_\pi^{\vec{p}_{\pi}}},
\end{align}
where $i$ indexes the single- and two-hadron interpolating operators used to construct the $\pi\pi$ state with $J^{PC}=1^{--}$ and $I=1$ \cite{Alexandrou:2018jbt}. The initial two-pion state has momentum $\vec{P}$, is in the $r$th row of the finite-volume irreducible representation $\Lambda$, and has energy $E_{n}^{\vec{P},\Lambda}$ and overlap factors $Z_i^{n,\vec{P},\Lambda}$. The final state has momentum $\vec{p}_{\pi}$ and energy $E_\pi^{\vec{p}_{\pi}}$ with the overlap factor $Z_\pi^{\vec{p}_{\pi}}$. The electromagnetic current insertion $J_\mu = Z_V (\frac{1}{3}\bar{u}\gamma_\mu u -\frac{2}{3} \bar{d}\gamma_\mu d)$ is renormalized as described in Ref.~\cite{Green:2017keo}. The multi-hadron three-point functions are a linear combination of the Wick contractions presented in Fig.~\ref{fig:3pt} and are calculated with forward, sequential and stochastic propagators. We have checked that the the current-disconnected diagrams are statistically compatible with $0$ and omit them from our construction.
\begin{figure}[htb!]
  \centering
  \includegraphics[width=0.29\textwidth]{wick_3pt_qq_disc.pdf}%
  \includegraphics[width=0.29\textwidth]{wick_3pt_MM_disc.pdf
  \includegraphics[width=0.29\textwidth]{wick_3pt_qq.pdf}
  \includegraphics[width=0.29\textwidth]{wick_3pt_T1.pdf}%
  \includegraphics[width=0.29\textwidth]{wick_3pt_T2.pdf}%
  \includegraphics[width=0.29\textwidth]{wick_3pt_T3.pdf}
  \caption{
  \label{fig:3pt} The different topologies of Wick contractions that make up the three-point function $C_{3,\;\mu,i}^{\vec{p}_{\pi},\,\vec{P},\,\Lambda,\,r}(t_{\pi},t_{J},t_{\pi\pi})$. The current-disconnected diagrams a) and b) are not included in the construction of the three-point functions.
  }
\end{figure}
The three-point functions as written in Eq.~(\ref{eq:C3_decomp}) are still a sum over all states in the $\pi\pi$ channel and need to be projected to the finite-volume states which have well-defined invariant mass $s$ and momentum transfer $q^2$. We achieve this by making use of the generalized eigenvectors $v_i^{n\,\vec{P}, \Lambda}(t_0)$ from the variational analysis in Ref.~\cite{Alexandrou:2017mpi}, and construct the optimized three-point functions \cite{Dudek:2009kk,Becirevic:2014rda,Shultz:2015pfa}:
\begin{align}
\label{eq:omega}
&\Omega_{3,\;\mu,\, n}^{\vec{p}_{\pi},\,\vec{P},\,\Lambda,\,r}(t_{\pi},t_{J},t_{\pi\pi}, t_0) =  v_i^{n\,\vec{P}, \Lambda}(t_0)\:C_{3,\;\mu,i}^{\vec{p}_{\pi},\,\vec{P},\,\Lambda,\,r}(t_{\pi},t_{J},t_{\pi\pi}) \cr
&= \sqrt{2E_n^{\vec{P}, \Lambda}}\, e^{E_n^{\vec{P}, \Lambda} t_0 / 2} \, Z_\pi^{\vec{p}_{\pi}} \, \langle \pi,\vec{p}_{\pi} |J_\mu(0,\vec{q}) | n, \vec{P}, \Lambda,r \rangle_{FV}\times\frac{ e^{-E_\pi^{\vec{p}_{\pi}} (t_{\pi} - t_{J})} e^{-E_n^{\vec{P}, \Lambda} (t_{J} - t_{\pi\pi})}}{2E_n^{\vec{P}, \Lambda} 2E_{\pi}^{\vec{p}_{\pi}}}.
\end{align}
Here, the sum over all states $n$ in the irreducible representation is reduced to a single term which enables us to use all individual excited states within the energy region of interest to determine the matrix element at as many points as allowed by the quark masses and volume of the gauge ensemble.

The matrix elements $\langle \pi,\vec{p}_{\pi} | J_\mu(0,\vec{q})|n,\vec{P},\,\Lambda,\,r\rangle_{FV}$ are determined using the ratio
\begin{align}
\label{eq:ratio}
& R_{\mu,\, n}^{\vec{p}_{\pi},\,\vec{P},\,\Lambda,\,r}(t_{\pi},t_{J},t_{\pi\pi}) = \frac{|\langle \pi,\vec{p}_{\pi} | J_\mu(0,\vec{q})|n,\vec{P},\,\Lambda,\,r\rangle_{FV}|^2}{4 E_n^{\vec{P}, \Lambda} E_\pi^{\vec{p}_{\pi}}} \\
& = \frac{\Omega_{3,\;\mu,\, n}^{\vec{p}_{\pi},\,\vec{P},\,\Lambda,\,r}(t_{\pi},t_{J},t_{\pi\pi}, t_0)\: \Omega_{3,\;\mu,\, n}^{\vec{p}_{\pi},\,\vec{P},\,\Lambda,\,r\,\dag}(t_{\pi}, t^\prime,t_{\pi\pi}, t_0)}{C_{\pi}^{\vec{p}_{\pi}}(\Delta t)\:\lambda_n^{\vec{P}, \Lambda}(\Delta t,t_0)},
\end{align}
which depends on the optimized three-point functions as well as the two point functions.

\section{Mapping to infinite volume}
To obtain the infinite-volume matrix elements we follow the Brice\~no-Hansen-Walker-Loud approach \cite{Briceno:2014uqa,Briceno:2015csa}, which is a generalization of the seminal work of Lellouch and L\"uscher \cite{Lellouch:2000pv} and allows us to utilize any gauge ensemble to study general transitions between single-hadron and two-hadron states. The mapping in the case of $\pi\gamma\to\pi\pi$ looks like
\begin{align}
\label{eq:LLmap}
& \frac{|\langle \pi, \vec{p}_\pi |J_\mu(0) |  s,q^2 ; \vec{P},\,\Lambda,\,r \rangle_{IV}|^2}{|\langle \pi,\vec{p}_{\pi} | J_\mu(0,\vec{q})|n,\vec{P},\,\Lambda,\,r\rangle_{FV}|^2}  = \frac{1 }{2E_n^{\vec{P},\Lambda}}    \frac{16 \pi \sqrt{s_n^{\vec{P},\Lambda}}}{k^{\vec{P},\Lambda}_n} \left( \frac{\partial \delta}{\partial E} + \frac{\partial \phi^{\vec{P},\Lambda}}{\partial E} \right) \bigg|_{E=E_n^{\vec{P},\Lambda}}, \qquad
\end{align}
where $\phi^{\vec{P},\Lambda}$ are the combinations of the finite-volume quantization condition (cf.~Eq.~49 of Ref.~\cite{Alexandrou:2017mpi}) and $\delta$ is the scattering phase shift as determined in Ref.~\cite{Alexandrou:2017mpi}. We make use of two different parametrizations of the scattering phase shift and perform the infinite-volume mapping for both. The systematic uncertainties associated with the different choices in fitting the energies and matrix elements are propagated to the final results.

\section{The $\pi\gamma \to \pi\pi$ amplitude}
We build the transition amplitude ${\cal V}_{\pi\gamma \to \pi\pi}$ (see Eq.~\ref{eq:V_BW_d}) from phase shifts BW I and BW II and parametrize the transition form factor $F$ with a two-dimensional Taylor series:
\begin{align}
F(q^2,s) = \frac{1}{1 - \frac{q^2}{m_P^2}} \sum_{n,m} A_{nm} z^n \mathcal{S}^m,
\label{eq:Fseries}
\end{align}
where ${\cal S} = \frac{s - m_R^2}{m_R^2}$, $z = \frac{\sqrt{4m_\pi^2 - q^2} - \sqrt{4m_\pi^2 }}{\sqrt{4m_\pi^2 - q^2} + \sqrt{4m_\pi^2 }}$, and $m_R$ corresponds to the $\rho$ resonance mass. Because we neglect the current-disconnected diagrams, we also set $m_P=m_R$. We investigate three different systematic truncations: {\bf F1)} combined order $K$, $\sum_{n+m \leq K}$,  {\bf F2)} order $N$ in $z$ and combined order $K$, $\sum_{n=0}^{N} \sum_{m=0}^{K-n}$, and {\bf F3)} order $N$ in $z$ and order $M$ in $\mathcal{S}$, $\sum_{n=0}^{N} \sum_{m=0}^{M}$. We present a three-dimensional projection of the transition amplitude ${\cal V}_{\pi\gamma \to \pi\pi}$ for the chosen parametrization ``BW II F1 K2'' in Fig.~\ref{fig:V_3d}.
\begin{center}
\begin{figure}[htb]
\centering
  \includegraphics[width=0.67\textwidth]{V_TII_F1_K2.pdf}
  \caption{\label{fig:V_3d} A three-dimensional plot of ${\cal V}_{\pi\gamma \to \pi\pi}$. The surface shows the central value of the nominal fit function (``BWII F1 K2'').  The lattice results are shown as the vertical bars, with the widths and depths corresponding to the uncertainties. The magenta sections in the vertical direction cover the range from ${\cal V}_{\pi\gamma \to \pi\pi}-\sigma_{{\cal V}_{\pi\gamma \to \pi\pi}}$ to ${\cal V}_{\pi\gamma \to \pi\pi}+\sigma_{{\cal V}_{\pi\gamma \to \pi\pi}}$. Data points with larger uncertainties have reduced opacity.}
\end{figure}
\end{center}

\section{Resonant form factor and photocoupling}

The resonant form factor is determined by analytically continuing the transition form factor $F$ to the $\rho$ pole:
\begin{equation}
 F_{\pi\gamma\to\rho}(q^2) = F(q^2,\, m_R^2 + i m_R \Gamma_R). \label{eq:Fresonant}
\end{equation}
At zero momentum transfer, $q^2=0$, the resonant form factor becomes equal to the photocoupling $G_{\rho\pi\gamma}$, which enters the formula for the $\rho$ radiative decay width $\Gamma(\rho\to\pi\gamma)$:
\begin{align}
&\Gamma(\rho \to \pi \gamma) =\frac{2}{3}\alpha \left(\frac{m_R^2 - m_\pi^2}{2m_R}\right)^3  \frac{|G_{\rho\pi\gamma}|^2}{m_\pi^2}. \label{eq:radwidth}
\end{align}
The left panel of Fig.~\ref{fig:PC} shows the resonant transition form factor as a function of momentum transfer. The inner shaded region corresponds to the statistical and systematical errors and the outer shaded region represents the model uncertainty, evaluated point-by-point as the root-mean-square deviation of the other parametrizations from the nominal parametrization.
\begin{figure}[!t]
\includegraphics[width=0.45\textwidth]{FF_VIIV_model_TII_F1_K2.pdf}%
\includegraphics[width=0.45\textwidth]{photocoupling_model_dep.pdf}
\caption{\label{fig:PC} Left panel: Resonant transition form factor $F_{\pi\gamma\to\rho}$ as a function of momentum transfer. Right panel: Absolute value of the photocoupling for the various parametrizations under consideration.}
\end{figure}
The right panel of Fig.~\ref{fig:PC} shows the photocoupling determined from those parametrizations that lead to a good description of the data without redundant parameters. We find the $\rho$ photocoupling at $m_\pi\approx320$ MeV to be: $|G_{\rho\pi\gamma}| = 0.0802(32)(20)$.
Because our lattice QCD calculation was performed at an unphysical pion mass, our thresholds are much closer to the resonances than in nature and thus our phase space is unphysical. However, if we assume the photocoupling to be pion-mass independent and use the physical values of hadron masses to have the proper thresholds, we find the radiative decay width $\Gamma(\rho\to\pi\gamma) \,=\, 84.2(6.7)(4.3)\, {\rm keV}$. The number in the first bracket represents the combined statistical and systematical uncertainty and the second uncertainty is associated with the parametrization.

\section{Conclusions}

In this talk we presented a $(2+1)$-flavor lattice QCD calculation of the $\rho$ photoproduction process at $m_\pi\approx 320$ MeV. Because at this light-quark mass the $\rho$ is a resonance, we made use of the the Brice\~no-Hansen-Walker-Loud formalism to calculate the general $\pi\gamma\to\pi\pi$ transition amplitude at many discrete points in the kinematic space of $q^2$ and $\sqrt{s}$. Fitting a multitude of parametrizations to the data and analytically continuing them to the $\rho$ pole enabled us to determine the $\rho$ photocoupling and radiative decay width.  Additional computations are needed to perform the continuum extrapolation and the chiral extrapolation to the physical point, where direct comparisons with experimental data are possible.
The methods used in this work are also applicable to other electroweak processes with two hadrons in the final state, such as $B \to \rho (\to \pi \pi) \ell \bar{\nu}$ and $B \to K^* (\to K \pi) \ell^+ \ell^-$.

\section*{Acknowledgments}
We are grateful to Kostas Orginos for providing the gauge field ensemble, which was generated using resources provided by XSEDE (supported by National Science  Foundation Grant No.~ACI-1053575).
LL was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under contract DE-AC05-06OR23177.
SM and GR were supported in part by National Science Foundation Grant No.~PHY-1520996; SM and GR were also supported in part by the U.S.~Department of Energy Office of High Energy Physics under Grant No.~DE-{}SC0009913.
SM and SS further acknowledge support by the RHIC Physics Fellow Program of the RIKEN BNL Research Center.
JN and AP were supported in part by the U.S. Department of Energy Office of Nuclear Physics under Grant Nos.~DE{}-SC-0{}011090 and DE-{}FC02-06ER41444.
We acknowledge funding from the  European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 642069. SP is a Marie Sklodowska-Curie  fellow supported by the HPC-LEAP joint doctorate program. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-{}AC02-05CH11231.

\providecommand{\href}[2]{#2}\begingroup\raggedright


\section{Introduction}\label{intro}
First a brief review of `classic' parametrizations of Bose-Einstein Correlations (BEC) is given and contrasted with the
parametrization of the $\tau$-model~\cite{Tamas;Zimanji:1990,ourTauModel}.
which has been found~\cite{L3_levy:2011} to describe well
Bose-Einstein correlations in hadronic \PZ\ decay.

The data used in this paper are taken from  Ref.~\citen{L3_levy:2011}. They comprise both two-jet and three-jet events,
as determined using the Durham jet algorithm~\cite{durham,durham2,durham3} with resolution parameter \ycut=0.006,
from \Pee\ annihilation at the \PZ-pole.

\subsection{`Classic' Parametrizations}    \label{classic}
The    Bose-Einstein correlation function, $R_2$,
is measured by $R_2(Q)=\rho(Q)/\rho_0(Q)$,
where $\rho(Q)$ is the density of identical boson pairs with invariant four-momentum difference
$Q=\sqrt{-(p_1-p_2)^2}$
and $\rho_0(Q)$ is the similar density in an artificially constructed reference sample,
which should differ from the data only in that it does not contain the effects of Bose symmetrization of identical bosons.
It is often parametrized as
\begin{equation} \label{eq:R2parametriz}
  R_2 =    \gamma \left[ 1+ \lambda G \right] (1+ \epsilon Q) \;,
\end{equation}
with
\begin{equation} \label{eq:gauss_param}
  G   =                             \exp \left(-\left(rQ\right)^{2} \right) \;.
\end{equation}
The corresponding distribution of boson emission points in space-time is a spherically symmetric Gaussian with
standard deviation $r$.

The factor $(1+\epsilon Q)$ is included to account for non-BEC which are not removed by $\rho_0$, \ie,
to make up for slight inadequacies in $\rho_0$, and $\gamma$ is a normalization parameter.
The parameter $\lambda$ is introduced to account for effects reducing the
amount of BEC, \eg, some of the bosons coming from resonance decays, or the identical bosons being partially coherent.

However, this description was found~\cite{L3_levy:2011} not to describe the \PZ-decay data, a better description being provided by
the Edgeworth expansion about the Gaussian~\cite{L3_3D:1999}:
\begin{equation} \label{eq:edgew_param}
  G   =                             \exp \left(-\left(rQ\right)^{2} \right)
           \cdot\left[1+\frac{\kappa}{3!} H_{3}(rQ)\right]  \;,
\end{equation}
where $H_3$ is the third-order Hermite polynomial.

A different way to depart from the Gaussian is the generalization to a symmetric L\'evy stable distribution. Then
\begin{equation} \label{eq:symlevy_param}
  G   =                             \exp \left(-\left(rQ\right)^{\alpha} \right) \:,
\end{equation}
where $0<\alpha\le2$ is the so-called index of stability, which was introduced to BEC in Ref.~\cite{Tamas:Levy2004}.

A fit of the Edgeworth parametrization to the two-jet data of Ref.~\citen{L3_levy:2011} finds
$\kappa=0.71\pm0.06$, while a fit of the symmetric L\'evy parametrization yields $\alpha=1.34\pm0.04$.
Both values are far from the corresponding Gaussian values of $\kappa=0$ and $\alpha=2$, respectively.
Although the \chisq\ of these fits are a great improvement over that of the Gaussian fit, they are still very large.
The corresponding confidence levels are $10^{-15}$ for the Gaussian and $10^{-5}$ and $10^{-8}$ for the Edgeworth and L\'evy
fits, respectively.   The symmetric L\'evy fit is shown in \Fig{fig:symlevy2jet}a.

\begin{figure}[h]
  \mbox{\includegraphics[width=.45\textwidth,clip,viewport=56 87  518 680]{run_fitQ1SL4-ISMDproc}}
  \mbox{\includegraphics[width=.45\textwidth,clip,viewport=56 87  518 680]{run_fitQR4-ISMDproc}}
        \caption{(a) Symmetric L\'evy fit for two-jet events;
                 (b) Simplified $\tau$-model\ fit for two-jet events.
        }
        \label{fig:symlevy2jet}
\end{figure}

The reason for the failure of the `classic' parametrizations is readily apparent from
\Fig{fig:symlevy2jet}a. There is a region of anti-correlation ($R_2<1$) extending from about $Q=0.5$ to $1.5\,\GeV$.
The `classic' parametrizations, which are of the form $1+P$, where $P$ is a positive-definite quantity, are unable
to accomadate the anti-correlation.
This was not realized for a long time because experiments only plotted the correlation function up to $Q=2\;\GeV$
or less.  In Ref.~\citen{L3_levy:2011} $Q$ was plotted to 4\,\GeV, and the
anti-correlation became apparent. This anti-correlation,
which one might term Bose-Einstein Anti-Correlations (BEAC),
as well as the Bose-Einstein correlations (BEC) are both well described by the $\tau$-model.

\subsection{The $\tau$-model}    \label{taumodel}
The $\tau$-model~\cite{Tamas;Zimanji:1990,ourTauModel}
rests on several assumptions:
\begin{itemize}
   \item The average production point is proportional to the momentum of the emitted boson.
     Dimensionally the momentum must be multiplied by a time divided by a mass to yield a spatial dimension.
     The description of a two-jet event is invariant to Lorentz boosts along the direction of the colour field.
     The relevant boost-invariant quantities are the ``longitudinal'' proper time,
     $\tau=\sqrt{\overline{t}^2-\overline{r}^2_z}$
     and the ``transverse'' mass,
     $\mt=\sqrt{E^2-p_z^2}$, resulting in
     \begin{equation} \label{eq:xpcorr}
        \overline{x}^\mu (p^\mu)  =  a \tau\ p^\mu \;, \qquad a=1/\mt \;.
     \end{equation}
   \item The spatial distribution of production points about their mean is very narrow,
     although the distribution of proper time may be broad.
   \item The distribution of $\tau$ is a one-sided L\'evy distribution,
     one-sided because no particles are emitted before the \Pee\ collision.
\end{itemize}
Then $R_2$ turns out to depend on three variables, $Q$ and the transverse mass of each of the particles making up the pair:
\begin{eqnarray}
   R_2 (Q, a_1, a_2) &=&
       \gamma
       \left\{ 1+\lambda
      \cos\left[\frac{\tau_0 Q^2 (a_1+a_2)}{2} +
           \tan\left(\frac{\alpha\pi}{2}\right)
           \left(\frac{\Delta\tau {Q^2}}{2}\right)^{\!\alpha}
           \frac{a_{1}^{\alpha}+a_{2}^\alpha}{2} \right]
  \right. \nonumber  \\
     & & \hfill
  \left.     \cdot
      \exp\left[-\left(\frac{\Delta\tau {Q^2}}{2}\right)^{\!\alpha}
          \frac{ a_{1}^{\alpha}+a_{2}^{\alpha}}{2} \right]
      \right\}
        \cdot\left(1+\epsilon Q\right) \;.
  \label{eq:R2taumodel}
\end{eqnarray}

Fits in three dimensions are problematic with the available statistics. Hence we simplify this expression by
introducing an effective radius, $R$, defined by
     \begin{equation} \label{eq:effR}
         R = \left(\frac{\Delta\tau}{2}\right)^{\!\alpha}
               \frac{a_{1}^{\alpha}+a_{2}^{\alpha}}{2}      \;.
     \end{equation}
Further, we assume that particle production begins immediately, \ie, $\tau_0=0$.
Then
\begin{subequations} \label{eq:asymlevR2}
\begin{align}
    R_2(Q) &= \gamma \left[ 1+ \lambda \cos \left(\left(\ensuremath{R_\mathrm{a}} Q\right)^{2\alpha} \right)
             \exp \left(-\left(RQ\right)^{2\alpha} \right) \right] (1+ \epsilon Q) \;, \label{eq:asymlevR2_}   \\
    R_\mathrm{a}^{2\alpha} &= \tan\left(\frac{\alpha\pi}{2}\right) R^{2\alpha} \;.     \label{eq:asymlevRaR}
\end{align}
\end{subequations}
The fit of the simplified $\tau$-model\ to the two-jet data is shown in \Fig{fig:symlevy2jet}b.
Unlike the fits of the classic parametrizations, the \chisq\ is acceptable, and the residuals lack structure.

Note that the difference between the parametrizations of \Eqs{eq:symlevy_param} and (\ref{eq:asymlevR2})
is the presence of the cosine term, which provides the description of the BEAC dip.
The parameter $R$ describes the BEC peak, and \ensuremath{R_\mathrm{a}}\  describes the anti-correlation region.
While one might have had the insight to add, \textit{ad hoc}, a $\cos$ term to \Eq{eq:symlevy_param}, it is the $\tau$-model\
which provides a physical reason for it and
which predicts a relationship,     \Eq{eq:asymlevRaR}, between $R$ and \ensuremath{R_\mathrm{a}}, \ie, between
the correlation and the anti-correlation.

\section{Expansions} \label{sect:expansion}
Recall that the Edgeworth expansion of the Gaussian parametrization provided evidence (in addition to the poor \chisq)
that the Gaussian was inadequate.
In this section we look at expansions of the Symmetric L\'evy and the $\tau$-model\ parametrizations.

\subsection{Symmetric L\'evy}
The symmetric L\'evy distribution can be expanded in terms of L\'evy polynomials~\cite{DeKock:WPCF2011,Csorgo:Odderon}, $l_i$,
which are orthonormal.
The resulting expression for $R_2$ is
\begin{equation} \label{eq:R2levypol}
         R_2(Q) = \gamma
                   \left[ 1+\lambda \exp\left(-\abs{rQ}^\alpha\right)
                             (1 + \sum c_i l_i )
                   \right]
                   \cdot\left(1+\epsilon Q\right) \;.
\end{equation}
Fits to the two-jet data are shown in \Fig{fig:R2expan}a for orders 0 through 3 of the L\'evy polynomials.
The order-0 fit (also shown in \Fig{fig:symlevy2jet}a) has a very poor \chisq, but the order-1 fit has a good \chisq.
Higher orders show only marginal further improvement.

\begin{figure}[h]
   \mbox{\includegraphics[width=.45\textwidth,clip,viewport=56 87  518 680]{run_fitQ1SL0-3npoly4}}
   \mbox{\includegraphics[width=.45\textwidth,clip,viewport=56 87  518 680]{run_fitQ1tau-SL1npoly4}}
   \raisebox{75mm}{\hspace{-.8\textwidth}(a)\hspace{.43\textwidth}(b)}
        \caption{(a) Symmetric L\'evy fit for two-jet events up to orders 0 through 3 L\'evy polynomials;
                 (b) First-order L\'evy polynomial fit (black) for two-jet events compared to the
                 $\tau$-model\ fit (red).
        }
   \label{fig:R2expan}
\end{figure}

The first-order symmetric L\'evy polynomial fit is shown together with the simplified $\tau$-model\ fit in
\Fig{fig:R2expan}b.  The \chisq\ of the L\'evy polynomial fit is slightly better than the $\tau$-model\ fit,
but the fit curves are nearly identical, what difference there is being mainly for $Q>1.5\,\GeV$.

Comparing \Eqs{eq:asymlevR2} and (\ref{eq:R2levypol}), we see that in the symmetric L\'evy parametrization
the cosine of the $\tau$-model\ parametrization is replaced by the L\'evy polynomial expansion.
Also, in the exponential $2\alpha$ becomes simply $\alpha$.
\Fig{fig:cospolycompare} compares the cosine and the L\'evy polynomials.  We see a rather similar behaviour:
Both decrease more or less linearly with $Q$, which explains why both fit the data approximately equally well.

\begin{figure}[h]
   \mbox{\includegraphics[width=.99\textwidth]{taumodel-levypol-diff-proc.pdf}}
   \caption{Comparison of the cosine term of the $\tau$-model\ parametrization with the L\'evy polynomial term of the
            symmetric L\'evy parametrization.}
   \label{fig:cospolycompare}
\end{figure}

\subsection{$\tau$-model\ (asymmetric L\'evy)}
Lacking an orthogonal polynomial expansion for the asymmetric L\'evy distribution $H(\tau)$ of the $\tau$-model,
we use, motivated by the results of Ref.~\citen{DeKock:WPCF2011}, a derivative expansion:
\begin{eqnarray}
     R_2(Q) &=& {\gamma}
        \left[ 1+ \lambda \left\{
        \cos \left(\left(\ensuremath{R_\mathrm{a}} Q\right)^{2\alpha}\right)
             \exp\left(-\left(RQ\right)^{2\alpha}\right)
          \right. \right. \qquad\qquad\qquad  \nonumber \\
           & &  \left. \left
             + \sum   c_n \frac{{\textrm d}^n\phantom{Q}}{{\textrm d}Q^n}
        \cos \left(\left(\ensuremath{R_\mathrm{a}} Q\right)^{2\alpha}\right)
             \exp\left(-\left(RQ\right)^{2\alpha}\right)
             \right\}
             \right]
             \cdot\left(1+ \epsilon Q\right)       \;.
\end{eqnarray}
We also consider letting \ensuremath{R_\mathrm{a}}\ be a free parameter rather than as defined in \Eq{eq:asymlevRaR}.
The results of these fits are shown in \Tab{tab:fits} and \Fig{fig:R2tauexpan}a.
The \chisq\ of the order-1 fit is, of course, smaller than that of the order-0 fit, as is the \ensuremath{R_\mathrm{a}}-free fit,
and the confidence levels are somewhat better.
To test the significance of the improvement in \chisq, we make use of the \chisq-difference.
For the order-0 and order-1 fits this is $\chisq_\text{diff}=94.7-90.9=3.8$, and the difference in the number of degrees of freedom
is 1.  The confidence level for a \chisq\ of 3.8 with 1 degree of freedom is 5.1\%.
A small value of this confidence level, say less than 5\%, would be grounds for rejecting the order-0 parametrization. (This
corresponds to the 95\% commonly used in making decisions.)
Thus we conclude that for two-jet events the order-0 fit is adequate.

\begin{table}[b]
\centering
\caption{Fit results of $\tau$-model\ parametrizations for two-jet events.}
\label{tab:fits}
   $ \begin{array}{lccc}
      \hline
              &\text{order 0}&\text{order 1}& \text{order 0, \ensuremath{R_\mathrm{a}}\ free} \\
      \hline
       \alpha & 0.44\pm0.01  & 0.43\pm0.01     & 0.41\pm0.02      \\
 R\ \text{(fm)}& 0.78\pm0.04  & 0.84\pm0.05     & 0.79\pm0.04     \\
 \ensuremath{R_\mathrm{a}}\ \text{(fm)}& -            & -             & 0.69\pm0.04     \\
       \lambda& 0.61\pm0.03  & 0.67\pm0.05     & 0.63\pm0.03      \\
       \gamma & 0.979\pm0.002& 0.979\pm0.002   & 0.988\pm0.005    \\
       \epsilon&0.005\pm0.001& 0.005\pm0.001   & 0.001\pm0.002    \\
        c_1   & -            & 0.0008\pm0.0005 & -                \\
      \hline
       \chisq/{\text{DoF}}\hspace{-2mm} &94.7/95 &90.9/94 &91.0/94\\
       \text{CL}& 49\%  & 57\% & 57\%  \\
     \hline
    \end{array} $
\end{table}

The order-0, \ensuremath{R_\mathrm{a}}-free fit also provides an adequate description, having a \chisq\ nearly identical to that of the order-1,
\ensuremath{R_\mathrm{a}}-constrained fit.
Note that the physical parameters ($\alpha$, $R$, $\lambda$) differ at most by about 1 standard deviation in going from
order-0 to order-1 or to \ensuremath{R_\mathrm{a}}\ free.  Thus conclusions based on these values, such as the reconstructions of the space-time picture
in Ref.~\citen{L3_levy:2011}, remain valid.

\begin{figure}[h]
   \mbox{\includegraphics[width=.45\textwidth,clip,viewport=56 87  518 680]{run_fitQR4-QR4c1-QRR4}}
   \mbox{\includegraphics[width=.45\textwidth,clip,viewport=56 87  518 680]{run_fitQR4-QR4c1-QRR4-QRR4c1_3j}}
   \raisebox{75mm}{\hspace{-.8\textwidth}(a)\hspace{.43\textwidth}(b)}
        \caption{(a) Fits of the $\tau$-model\ (black), the order-1 expansion of the $\tau$-model\ (red), and the $\tau$-model\ with \ensuremath{R_\mathrm{a}}\
                     free (green), for two-jet events;
                 (b) Fits of the $\tau$-model\ (black), the order-1 expansion of the $\tau$-model\ (red), the $\tau$-model\ with \ensuremath{R_\mathrm{a}}\ free
                     (green), and the order-1 expansion of the $\tau$-model\ with \ensuremath{R_\mathrm{a}}\ free (blue), for three-jet events.
        }
   \label{fig:R2tauexpan}
\end{figure}

\section{Conclusions for two-jet events} \label{sect:2jetConclusions}
We have used expansions about the hypothesized form to test whether it provides an adequate description
of the data or is only a (poor) approximation.  In the latter case the shape of the expansion terms provide an indication of how
to modify the original parametrization.  For the symmetric L\'evy parametrization this showed that an approximately linearly
decreasing function of $Q$ is necessary, which in fact is what is provided by the $\tau$-model.

An expansion in the case of the $\tau$-model\ is found not to be significant, \ie, the one-sided L\'evy distribution of the $\tau$-model\
is adequate.

\section{Three-jet events} \label{sect:3jet}
For two-jet events hadronization occurs basically in 1+1 dimensions, which lead to the dependence of $R_2$ on
$\tau$, the longitudinal proper time and \mt, the transverse mass.
For three-jet events, the \Pq\Paq\Pg\ system no longer forms a linear system (in the overall centre of mass), but a planar one.
There is no event axis by which the transverse mass and longitudinal proper time are defined.
Therefore we might expect the $\tau$-model, as formulated for a two-jet system, to work less well.

The results of fits of the $\tau$-model\ and its first-order expansion, without and with \ensuremath{R_\mathrm{a}}\ a free parameter, are shown in
\Fig{fig:R2tauexpan}b and \Tab{tab:fits3j}.  Applying the \chisq-difference test to the order-0 and order-1 fits yields a
confidence level of 37\% for the case that \ensuremath{R_\mathrm{a}}\ is constrained and 58\% when \ensuremath{R_\mathrm{a}}\ is free.

However, the \chisq-difference between the order-0 \ensuremath{R_\mathrm{a}}\ constrained and free cases yields a confidence level of $6\cdot10^{-8}$.
Thus regarding \ensuremath{R_\mathrm{a}}\ as a free parameter does give significant improvement.

But it must be pointed out that the parameters \ensuremath{R_\mathrm{a}}\ and
$\tau_0$ are expected to be highly correlated. While for two-jet events it was found in Ref.~\citen{L3_levy:2011} that $\tau_0$ is
consistent with zero, such studies have not yet been performed for three-jet events.
Further investigation is ongoing.
Also, note that the value of $\alpha$ is significantly less for the fits with \ensuremath{R_\mathrm{a}}\ free than for those with \ensuremath{R_\mathrm{a}}\ constrained.
This was not the case for two-jet events.  This too requires additional investigation.

\begin{table}[t]
\centering
\caption{Fit results of $\tau$-model\ parametrizations for three-jet events.}
\label{tab:fits3j}
   $ \begin{array}{lcccc}
      \hline
              &\text{order 0}&\text{order 1}& \text{order 0, \ensuremath{R_\mathrm{a}}\ free} & \text{order 1, \ensuremath{R_\mathrm{a}}\ free} \\
      \hline
       \alpha & 0.42\pm0.01  & 0.42\pm0.01     & 0.35\pm0.01       & 0.35\pm0.01 \\
 R\ \text{(fm)}& 0.98\pm0.04  & 0.96\pm0.05     & 1.06\pm0.05      & 1.00\pm0.11 \\
 \ensuremath{R_\mathrm{a}}\ \text{(fm)}& -            & -             & 0.87\pm0.04      & 0.96\pm0.30 \\
       \lambda& 0.84\pm0.04  & 0.81\pm0.05     & 0.92\pm0.05       & 0.65\pm0.49 \\
       \gamma & 0.977\pm0.001& 0.977\pm0.001   & 0.997\pm0.005     & 0.994\pm0.007 \\
       \epsilon&0.008\pm0.001& 0.008\pm0.001   & 0.0003\pm0.0017    & 0.001\pm0.002 \\
        c_1   & -            &-0.0005\pm0.0005 & -                 &-0.04 \pm0.08  \\
      \hline
       \chisq/{\text{DoF}}\hspace{-2mm} &113.2/95 &112.4/94 &83.7/94  & 83.4/93 \\
       \text{CL}& 10\%  &  9\% & 77\%    & 75\% \\
     \hline
    \end{array} $
\end{table}

\section{Conclusions for three-jet events} \label{sect:3jetConclusions}
As for the two-jet case, expansion of the $\tau$-model\ expression does not lead to significant improvement in the fits.
This validates the use of an asymmetric L\'evy distribution for the longitudinal proper time.

However, significant improvement of the fit is obtained by letting \ensuremath{R_\mathrm{a}}\ be a free parameter. \ie,
by lessening the connection of the simplified $\tau$-model\ between the BEC peak and the antisymmetric dip.
Whether letting $\tau_0$ also be a free parameter would also give significant improvement is the subject of ongoing investigation,
as is the question whether $\alpha$ decreases as the number of jets increases.

\begin{acknowledgement}
S.~L\"ok\"os greatfully acknowledges the support from the ERASMUS project of the EU and thanks W. Metzger for
his kind hospitality at the Radboud University Nijmegen, The Netherlands.
This research was partially supported by the  NKTIH FK 123842 and FK 123 959 grants (Hungary),
by the Circles of Knowledge Club (Hungary) and by  the  EFOP 3.6.1-16-2016-00001 project (Hungary).
\end{acknowledgement}

\newpage
\input p.bbl

\end{document}


\section{Grzegorczyk representation of natural numbers and Grzegorczyk sequences}
Let us represent natural numbers in terms of Grzegorczyk functions $F_{n}$.
The representation is implicit in the technique of H. Friedman\cite{Simpson85}.

\begin{definition}\label{df:Grzegorczyk}
{\rm
\begin{enumerate}
\item
$f^{(x)}$ denotes the $x$-th iterate of a unary function $f$ on $\mathbb{N}$, defined by
$f^{(0)}(y)=y$ and $f^{(x+1)}(y)=f(f^{(x)}(y))$.
\item
The \textit{$n$-th Grzegorczyk function}
$F_{n}$ denote the function on integers defined recursively on $n$ as follows.
$F_{0}(x)=x+1$ and $F_{n+1}(x)=F_{n}^{(x)}(x)$.
\end{enumerate}
}
\end{definition}

\begin{proposition}\label{prp:Grzegorczyk}
\begin{enumerate}
\item\label{prp:Grzegorczyk.1}
$F_{n}(x)>x$ for $x>0$.
\item\label{prp:Grzegorczyk.2}
$F_{n}^{(y)}(x)\geq x$.
\item\label{prp:Grzegorczyk.3}
$(x,y)\mapsto F_{n}^{(y)}(x)$ is strictly increasing both in $x>0$ and $y$.
\item\label{prp:Grzegorczyk.4}
$F_{n+1}(x)> F_{n}(x)$ for $x\geq 2$.
\end{enumerate}
\end{proposition}

\begin{lemma}\label{lem:Grzegorczyk}
For each primitive recursive function $f(\vec{x})$, there exists an $n$ for which
$f(\vec{x})\leq F_{n}(\max\{\vec{x},2\})$ holds for any $\vec{x}=x_{1},\ldots,x_{k}$.
\end{lemma}

\begin{definition}\label{df:Frep}
{\rm
Let $k_{1}\geq 2$.
$F_{n}$ denotes the $n$-th Grzegorczyk function.
We define \textit{$F$-representation}
of natural numbers $x$ \textit{with base $k_{1}$}
by induction on $x\dot{-}k_{1}$
as follows.

First if $x<k_{1}$, then $x$ itself is the $F$-representation of $x$ with base $k_{1}$.
Assume $k_{1}\leq x$.
Let $x_{1}$ be the least natural number such that $k_{1}\leq x<F_{x_{1}+1}(k_{1})$, and
$i_{1}$ the maximal number such that
$F_{x_{1}}^{(i_{1})}(k_{1})\leq x<F_{x_{1}}^{(k_{1})}(k_{1})=F_{x_{1}+1}(k_{1})$.
If $x_{1}=0$, then $[(0,i_{1})]_{k_{1}}$ with $0=x_{1}$ is the $F$-representation of $x$.
Otherwise let $x_{2}<x_{1}$ be the least natural number such that
$k_{2}\leq x<F_{x_{2}+1}(k_{2})$ for $k_{2}=F_{x_{1}}^{(i_{1})}(k_{1})\leq x<F_{x_{1}}(k_{2})$.
If $x=k_{2}$, then $[(x_{1},i_{1})]_{k_{1}}$ is the $F$-representation of $x$ with base $k_{1}$.
Otherwise
$[(x_{1},i_{1}),(x_{2},i_{2}),\ldots,(x_{\ell},i_{\ell})]_{k_{1}}$ is the
$F$-representation of $x$ with base $k_{1}$,
where $[(x_{2},i_{2}),\ldots,(x_{\ell},i_{\ell})]_{k_{2}}\neq[(0,0)]_{k_{2}}$ is the $F$-representation of $x$ with base $k_{2}$.
}
\end{definition}

We define \textit{Grzegorczyk sequences} $\{z_{k}\}_{k}$ of natural numbers $z$ as follows.
$z_{k+1}=z_{k}\dot{-}1$ if $z_{k}<2+k$.
Let $z_{k}\geq 2+k$.
Let $x=[(x_{1},i_{1}),(x_{2},i_{2}),\ldots,(x_{\ell},i_{\ell})]_{k}$ be the $F$-representation of
$x\geq k$ with base $k$.
Then the shift of base $k$ to $m$ in the representation as well as the representations of $x_{i}$
hereditarily in the left is defined recursively by
$x[k:=m]=
[(x_{1}[k:=m],i_{1}),(x_{2}[k:=m],i_{2}),\ldots,(x_{\ell}[k:=m],i_{\ell})]_{m}$.
Then
$z_{k+1}=z_{k}[2+k:=3+k]-1$.

Since the function $(n,x)\mapsto F_{n}(x)$ is provably computable in {\sf PA}
(but not primitive recursive), so is the function $(x,k,m)\mapsto x[k:=m]$.

\begin{theorem}\label{th:GrzegorczykseqPA}
${\sf PA}$ does not prove the true $\Pi_{2}$-statement $\forall z \exists k[z_{k}=0]$
for the Grzegorczyk sequence $\{z_{k}\}_{k}$ of $z$.
\end{theorem}

\section{Slowly well-foundedness}
\begin{definition}\label{df:normepsilon}
{\rm
We define integers $C(\alpha)$ for ordinal terms $\alpha<\varepsilon_{0}$.
$C(0)=0$.
In the following let $\alpha=\omega^{\alpha_{1}}n_{1}+\cdots+\omega^{\alpha_{k}}n_{k}$
with $\alpha_{1}>\cdots>\alpha_{k}$ and $0<n_{1},\ldots,n_{k}<\omega$ for $k>0$.
$C(\alpha)=\max\{C(\alpha_{\ell}), n_{\ell}:1\leq\ell\leq k\}$.
$C(\alpha)$ is the maximal coefficient $n_{\ell}$ in the ordinal term $\alpha$.
}
\end{definition}

In the $F$-representation of Definition \ref{df:Frep},
for the case $k_{1}\leq x$,
note that $k_{1}<F_{n}(k_{1})<F_{n+1}(k_{1})$ holds by $k_{1}\geq 2$ and Proposition \ref{prp:Grzegorczyk}.\ref{prp:Grzegorczyk.1}.
When $x_{1}>0$, $i_{1}>0$ holds by the minimality of $x_{1}$.
Hence $k_{1}<F_{x_{1}}^{(i_{1})}(k_{1})=k_{2}$.

It is clear that for $x\geq k_{1}$, $[(0,0)]_{k_{1}}$ is the $F$-representation of $x$ iff $x=k_{1}$.
Let $[(x_{1},i_{1}),(x_{2},i_{2}),\ldots,(x_{\ell},i_{\ell})]_{k_{1}}\neq[(0,0)]_{k_{1}}$
 be the $F$-representation of
$x$ with base $k_{1}$ for $2\leq k_{1}< x$.
Proposition \ref{prp:Grzegorczyk}.\ref{prp:Grzegorczyk.4} yields $x_{1}<F_{x_{1}}(k_{1})\leq x$.
Hence $x>x_{1}>x_{2}>\cdots>x_{\ell}\geq 0$,
$\forall p\leq\ell(i_{p}>0)$.

Let $k_{\ell}\, (1\leq p\leq\ell+1)$ denote the integers recursively defined as follows.
$k_{1}$ is the given number, and
$k_{p+1}=F_{x_{p}}^{(i_{p})}(k_{p})$.
Then
$i_{p}<k_{p}$, and $k_{1}<k_{2}<\cdots<k_{\ell}<k_{\ell+1}=x$, i.e.,
$x=F_{x_{\ell}}^{(i_{\ell})}(\cdots(F_{x_{1}}^{(i_{1})}(k_{1}))\cdots)$.
Moreover $\forall p\leq \ell(i_{p}<k_{p}<x)$.
Hence we obtain in the $F$-representation $[(x_{1},i_{1}),(x_{2},i_{2}),\ldots,(x_{\ell},i_{\ell})]_{k_{1}}$
of $x\geq k_{1}\geq 2$
\begin{equation}\label{eq:Frepless}
x>x_{1}>x_{2}>\cdots>x_{\ell}\geq 0 \,\&\, \forall p\leq \ell(i_{p}<x)
\end{equation}
Note that the function assigning $x\geq k$ to its $F$-representation
\\
$[(x_{1},i_{1}),(x_{2},i_{2}),\ldots,(x_{\ell},i_{\ell})]_{k_{1}}$
is elementary recursive since the ternary relation $R=\{(n,x,y)\in\mathbb{N}^{3}:F_{n}(x)=y\}$
is elementary recursive.
This is seen as follows.
By Proposition \ref{prp:Grzegorczyk}.\ref{prp:Grzegorczyk.1} we see that
$(n,x,y)\in R$ iff there exists a matrix $A=(a_{ij})$ such that its size is at most $n\times y$,
and its entry $a_{ij}\leq\max\{x,y\}$ is a number $F_{m}(y_{ij})$ appearing in the computation of
$F_{n-1}^{(p)}(x)\,(p\leq x)$.

Thus $x$ is represented by
\[
[(x_{1},i_{1}),(x_{2},i_{2}),\ldots,(x_{\ell},i_{\ell})]_{k_{1}}:=
F_{x_{\ell}}^{(i_{\ell})}(\cdots(F_{x_{1}}^{(i_{1})}(k_{1}))\cdots)=x
.\]

Let
$[(x_{1},i_{1}),(x_{2},i_{2}),\ldots,(x_{\ell},i_{\ell})]_{k_{1}}$ and
$[(y_{1},j_{1}),(y_{2},j_{2}),\ldots,(y_{m},j_{m})]_{k_{1}}$ be the $F$-representations of
$x\geq k_{1}$ and $y\geq k_{1}$, resp.
Then we see the following.
\begin{equation}\label{eq:Frep}
x<y \Leftrightarrow
\left((x_{1},i_{1}),(x_{2},i_{2}),\ldots,(x_{\ell},i_{\ell})\right)
\prec
\left((y_{1},j_{1}),(y_{2},j_{2}),\ldots,(y_{m},j_{m})\right)
\end{equation}
where for the lexicographic order $<_{lx}$ on pairs,
\[
\left((x_{1},i_{1}),(x_{2},i_{2}),\ldots,(x_{\ell},i_{\ell})\right)
\prec
\left((y_{1},j_{1}),(y_{2},j_{2}),\ldots,(y_{m},j_{m})\right)
\]
iff either
$\exists p\leq\min\{\ell,m\}[\forall q<p\{(x_{q},i_{q})=(y_{q},j_{q})\}\land (x_{p},i_{p})<_{lx}(y_{p},j_{p})]$
or $\ell<m\land \forall q\leq\ell\{(x_{q},i_{q})=(y_{q},j_{q})\}$.

Let $k_{1}\leq x<F_{n}(k_{1})$, and
$[(x_{1},i_{1}),(x_{2},i_{2}),\ldots,(x_{\ell},i_{\ell})]_{k_{1}}$ the $F$-representation of $x$
with base $k_{1}$, where $n> x_{1}>x_{2}>\cdots>x_{\ell}\geq 0$.
$s(x;n,k_{1}):=(j_{1},j_{2},\ldots,j_{n})$ denotes the sequence of natural numbers $j_{q}$,
which is obtained from the sequence $(i_{1},i_{2},\ldots,i_{\ell})$ by filling the gaps with zeros,
i.e., $j_{x_{p}}=i_{p}$ and $j_{q}=0$ else.
Moreover let $(m_{1},m_{2},\ldots,m_{n})$ be the sequence defined recursively by
$m_{1}=k_{1}$ and $m_{q+1}=F_{n-q}^{(j_{q})}(m_{q})$.

Let $y_{1},y_{2}$ be integers such that $k_{1}\leq y_{1},y_{2}<F_{n}(k_{1})$ for $n>0$.
For simplicity let us write $s(x)$ for $s(x;n,k_{1})$.
We see from Proposition \ref{prp:Grzegorczyk}.\ref{prp:Grzegorczyk.1} that for $k_{1}\geq 2$
\begin{equation}\label{eq:slows}
y_{1}<y_{2} \Leftrightarrow s(y_{1})<_{lx}s(y_{2})
\end{equation}
where $<_{lx}$ denotes the lexicographic order on $n$-tuples.

Now let us flip over the integers $i_{\ell}$.
Namely let ${}^{t}s(x)=(m_{1}-j_{1},m_{2}-j_{2},\ldots,m_{n}-j_{n})$ for
$s(x)=(j_{1},j_{2},\ldots,j_{n})$.
${}^{t}s(x)$ is an $n$-tuple of positive integers $m_{q}-j_{q}\leq x$ such that
\[
y_{1}<y_{2}\Leftrightarrow {}^{t}s(y_{2})<_{lx} {}^{t}s(y_{1})
\]
We obtain a descending chain $\{g_{n}(k_{1},x)\}_{x<F_{n}(k_{1})}$ of ordinals
$g_{n}(k_{1},x)<\omega^{n+1}$ for $n>0$ as follows.
When $k_{1}\leq x<F_{n}(k_{1})$, let
\begin{equation}\label{eq:gn}
g_{n}(k_{1},x)=\omega^{n-1}(m_{1}-j_{1})+\omega^{n-2}(m_{2}-j_{2})+\cdots+\omega^{0}(m_{n}-j_{n})
\end{equation}
for $s(x;n,k_{1})=(j_{1},j_{2},\ldots,j_{n})$.
When $x<k_{1}$, let
$g_{n}(k_{1},x)=\omega^{n}(k_{1}-x)$.

Let us extend $g_{n}$ by letting $g_{n}(k_{1},x)=0$ for $x\geq F_{n}(k_{1})$, and
$g_{0}(k_{1},x)=(k_{1}+1)\dot{-}x$.
Then it is easy to see that for each $n$ there exists a constant $c_{n}$ for which
$C(g_{n}(k_{1},x))\leq\max\{n,k_{1}+1,x\}$ holds
for any $k_{1},x$.

\begin{lemma}\label{lem:slowdown}{\rm (H. Friedman\cite{Simpson85})}\\
For each unary primitive recursive function $f$, there exist $n, c$ and a primitive recursive function
$g:\mathbb{N}^{2}\to\omega^{n}$ such that
$x<f(k) \Rarw g(k,x+1)<g(k,x)$, and
$C(g(k,x))\leq\max\{c,k+1,x\}$ for any $k,x$.
\end{lemma}
{\bf Proof}.\hspace{2mm}
Pick an $n$ such that $f(k)\leq F_{n}(\max\{2,k\})$ by Lemma \ref{lem:Grzegorczyk}.
\hspace*{\fill} $\Box$

\begin{theorem}\label{th:slowdown}
Let $\varepsilon_{0}>\alpha_{0}>\alpha_{1}>\cdots$ be a primitive recursive descending chain of ordinal terms $\alpha_{k}$.
We can find a primitive recursive descending chain $\{\gamma_{i}\}_{i}$ of ordinal terms such that
$C(\gamma_{i})\leq i+1$.
\end{theorem}
{\bf Proof}.\hspace{2mm}
Consider the primitive recursive functions
$f_{C}(k)=C(\alpha_{k+1})$.
By Lemma \ref{lem:slowdown},
let $h:\mathbb{N}^{2}\to\omega^{n}$ be a primitive recursive function and $c$ a constant such that
$x<C(\alpha_{k+1}) \Rarw h(k,x+1)<h(k,x)$, $C(h(k,x))\leq\max\{c,k+1,x\}$.

Llet $\ell=\max\{c,C(\alpha_{0})\}$.
For $i\geq \ell$, let $k,x$ be numbers such that
$i=C(\alpha_{0})+\cdots+C(\alpha_{k})+x$ and $x<C(\alpha_{k+1})$.
Then set
\[
\gamma_{i}  =  \omega^{\omega}\alpha_{k}+h(k,x)
\]
Note here that $C(\alpha_{k})>0$ for $\alpha_{k}>0$.

Let $N$ be a number such that $\omega_{N-\ell}>\omega^{\omega}\alpha_{0}$.
For $i<\ell$ let
$\gamma_{i}=\omega_{N-i}$.
Then the primitive recursive sequence $\{\gamma_{i}\}_{i}$ is descending, and $C(\gamma_{i})\leq i+1$.
\hspace*{\fill} $\Box$

\begin{cor}\label{cor:slowdown}
Over {\sf PRA}, the 1-consistency $\mbox{{\rm RFN}}_{\Sigma_{1}}({\sf PA})$ of {\sf PA}
is equivalent to the fact that
there is no primitive recursive and
infinite descending chain $\{\gamma_{i}\}_{i}$ of ordinals $\gamma_{i}<\varepsilon_{0}$
such that $C(\gamma_{i})\leq i+1$.
\end{cor}
{\bf Proof}.\hspace{2mm}
This is seen from the fact that $\mbox{{\rm RFN}}_{\Sigma_{1}}({\sf PA})$ is equivalent to the fact that
there is no primitive recursive and
infinite descending chain of ordinals$<\varepsilon_{0}$.
\hspace*{\fill} $\Box$

\section{Proof of Theorem \ref{th:GrzegorczykseqPA}}

For $x\geq k\geq 2$,
let $x=[(x_{1},i_{1}),(x_{2},i_{2}),\ldots,(x_{\ell},i_{\ell})]_{k}$
be the $F$-representation of $x$ with base $k$.
Then an ordinal $o_{k}(x)<\varepsilon_{0}$ is associated recursively by
\[
o_{k}(x)=\omega^{o_{k}(x_{1})}i_{1}+\omega^{o_{k}(x_{2})}i_{2}+\cdots+\omega^{o_{k}(x_{\ell})}i_{\ell}
\]
From (\ref{eq:Frep}) we see that for $x,y\geq k\geq 2$,
$o_{k}(x)$ is in Cantor normal form with base $\omega$ and
\begin{equation}\label{eq:Disoleft}
x<y\Leftrightarrow o_{k}(x)<o_{k}(y)
\end{equation}
Let
\[
D_{k}:=\{o_{k}(x):x\geq k\}
\]
and $L_{k}(\alpha)$ be the inverse of $o_{k}$ for $\alpha\in D_{k}$, $L_{k}(o_{k}(x))=x$.

Let us prove Theorem \ref{th:GrzegorczykseqPA}.
First we show the termination of Grzegorczyk sequences $\{z_{k}\}_{k}$ of $z\geq 2$.
Suppose $\forall k(z_{k}>0)$. We have $\forall k(z_{k}\geq 2+k)$, for
if $z_{k}=m<2+k$, then $z_{k+m}=0$.
Hence we obtain $\forall k(0<o_{2+k}(z_{k})\in D_{2+k})$.

From (\ref{eq:Disoleft}) we see that $\omega$ is the order type of the set $D_{k}$ of ordinals$<\varepsilon_{0}$ when $k\geq 2$.
For $0<\alpha\in D_{k}$, let $Q_{k}\alpha=\max\{\beta<\alpha: \beta\in D_{k}\}$.
Then it is easy to see that
$z_{k+1}=L_{3+k}Q_{3+k}o_{2+k}(z_{k})$ from $o_{2+k}(z_{k})\in D_{2+k}\subset D_{3+k}$.
Therefore $\{o_{2+k}(z_{k})\}_{k}$ would be an infinite descending chain of ordinals$<\varepsilon_{0}$.
Thus we show the fact that each Grzegorczyk sequences $\{z_{k}\}_{k}$ of $z$ eventually terminates.

Next let us show the independence following M. Rathjen\cite{Rathjen15}.
For the independence of the fact from ${\sf PA}$, it suffices to show that
the fact implies the $1$-consistency of ${\sf PA}$ over {\sf PA}.
Argue in {\sf PA}.
Suppose that ${\sf PA}$ is $1$-inconsistent.
By Corollary \ref{cor:slowdown}
let $\{\gamma_{k}\}_{k}$ be a primitive recursive and infinite descending chain of
ordinals $\gamma_{i}<\varepsilon_{0}$ such that $C(\gamma_{k})\leq 1+k$.
Then $\gamma_{k}\in D_{2+k}$ since $k+2<F_{x}^{(i)}(k+2)$ for $i>0$.
Hence $\{\gamma_{k}\}_{k}$ would be a primitive recursive and infinite descending chain of
ordinals$<\varepsilon_{0}$ such that $\gamma_{k}\in D_{2+k}$.
Let $v_{k}=L_{2+k}\gamma_{k}$, and $\{z_{k}\}_{k}$ be the Grzegorczyk sequence of $z_{0}=v_{0}$.
We show that $\forall k(v_{k}\leq z_{k})$ by induction on $k$.
Suppose $v_{k}\leq z_{k}$. Then $\gamma_{k}=o_{2+k}(v_{k})\leq o_{2+k}(z_{k})$.
From $\gamma_{k},\gamma_{k+1}\in D_{3+k}$ and $\gamma_{k+1}<\gamma_{k}\leq o_{2+k}(z_{k})$ we obtain
$\gamma_{k+1}\leq Q_{3+k}o_{2+k}(z_{k})$.
Hence $v_{k+1}=L_{3+k}\gamma_{k+1}\leq L_{3+k}Q_{3+k}o_{2+k}(z_{k})=z_{k+1}$.
Thus $\forall k(v_{k}\leq z_{k})$ is shown.
Since $\gamma_{k}> 0$, we have $v_{k}>0$.
Thus $\{z_{k}\}_{k}$ would be a non-terminating Grzegorczyk sequence of $z_{0}$.
\hspace*{\fill} $\Box$
\\

\noindent
The \textit{total representation}  $T(x)_{k_{1}}$ of $x$ \textit{with base} $k_{1}\geq 2$
represents $x$  hereditarily in the both components $x_{p},i_{p}$.
These are obtained by modifying Definition \ref{df:Frep} of $F$-representations as follows.

Let $k\geq 2$.
First if $x<k$, then $T(x)_{k}=x$ itself is
the hereditary $F$-representations of $x$ with base $k$.
Assume $k\leq x$, and let $[(x_{1},i_{1}),(x_{2},i_{2}),\ldots,(x_{\ell},i_{\ell})]_{k}$
be the $F$-representation of $x$ with base $k$.
Define hereditary representations recursively as follows, cf.\,(\ref{eq:Frepless}).
\\
$T(x)_{k}:=[(T(x_{1})_{k},T(i_{1})_{k}),(T(x_{2})_{k},T(i_{2})_{k}),\ldots,(T(x_{\ell})_{k},T(i_{\ell})_{k})]_{k}$,
\\
where $T(x_{p})_{k}, T(i_{p})_{k}$ are the total representations of $x_{p},i_{p}$ with base $k$.

We define \textit{hereditarily Grzegorczyk sequences} $\{w_{k}\}_{k}$ of natural numbers $w$ as follows.
$w_{k+1}=w_{k}\dot{-}1$ if $w_{k}<2+k$.
Let $w_{k}\geq 2+k$.

Let $x=[(x_{1},i_{1}),(x_{2},i_{2}),\ldots,(x_{\ell},i_{\ell})]_{k}$ be the hereditary $F$-representation of
$x\geq k$ with base $k$.
Then the shift of base $k$ to $m$ in the \textit{both} representations of $x_{p}, i_{p}$
hereditarily is defined recursively by
$x[k:=m]=
[(x_{1}[k:=m],i_{1}[k:=m]),(x_{2}[k:=m],i_{2}[k:=m]),\ldots,(x_{\ell}[k:=m],i_{\ell}[k:=m])]_{m}$.
Then
$w_{k+1}=w_{k}[2+k:=3+k]-1$.
\\

\noindent
{\bf Open problem}.
How strong the true statement $\forall w\exists k[w_{k}=0]$ for the hereditary Grzegorczyk sequence $\{w_{k}\}$ of $w$?
Specifically
does ${\sf ID}_{1}$ prove the $\Pi_{2}$-statement?