text:
\section{Introduction}
\label{sec:intro}
Supersonic turbulence is a key ingredient in various astrophysical
contexts, from gamma ray bursts~\citep{lazar-et-al:09,
  narayan-kumar:09} or stellar accretion~\citep{walder-folini:08,
  hobbs-et-al:11} to molecular clouds and star
formation~\citep{chabrier-hennebelle:11, federrath-klessen:12,
  padoan-et-al:12, kritsuk-et-al2:13}. A key question is whether this
turbulence, like incompressible turbulence, is characterized by
universal statistics. Results from 3D periodic box simulations of
driven, isothermal, supersonic turbulence~\citep{kritsuk-et-al:07,
  schmidt-et-al:08, pan-et-al:09} are indeed consistent with the
highly compressible variant~\citep{boldyrev:02} of the hierarchical
structure model that was put forward by~\citet{she-leveque:94} for
incompressible turbulence and that was further scrutinized
by~\citet{dubrulle:94} and~\citet{she-waymire:95}. This model is
correspondingly popular in astrophysics. It is employed, for example,
in the interpretation of molecular cloud
observations~\citep{gustafsson-et-al:06, hily-blant-et-al:08} or to
derive a theoretical expression for the density distribution in
supersonic turbulence~\citep{boldyrev-et-al1:02}, which enters
theories of the stellar initial mass
function~\citep{hennebelle-chabrier:08}.

Nevertheless, some doubts remain whether the model really applies to
simulation data of supersonic turbulence and, if so, with what
parameter values. The best-fit model parameters that we are aware
of~\citep{pan-et-al:09} still come with a 5-10\% uncertainty range
that is only marginally compatible with theoretically predicted
parameter values (see below). Here we argue that today's computer
resources should allow for 1-2\% accurate parameter fits in the highly
compressible regime, thereby likely settling the issue. Our claim is
based on a Monte Carlo approach to mimic actual simulation data.

The hierarchical structure model predicts the ratios $Z_{p}$ of
(observable) structure function scaling exponents $\zeta_{p}$,
$p=1,2,3...$ etc., of a 3D velocity field $\mathbf{u}$ as
\begin{equation}
Z_{p} = \frac{\zeta_{p}}{\zeta_{3}} =
(1 - \Delta)\frac{p}{3} + \frac{\Delta}{1-\beta}\left( 1 - \beta^{p/3} \right).
\label{eq:dubrulle}
\end{equation}
Here, $D=3-C$ is the dimension of the most intermittent structure,
$C=\Delta/(1-\beta)$ the associated co-dimension, $\beta \in
[0,1]$ measures the intermittency of the energy cascade, and
$\Delta \in [0,1]$ measures the divergent scale dependence of
the most intermittent structures. The $\zeta_{p}$ are defined in the inertial range by
\begin{equation}
S_{p}(r) \equiv < | \mathbf{u}(\mathbf{x} + \mathbf{r}) - \mathbf{u}(\mathbf{x})  |^{p} >
         \,\,\, \propto r^{\zeta_{p}},
\label{eq:struc_func_def}
\end{equation}
where $<\dots>$ denotes the average over all positions $\mathbf{x}$
within the sample and over all directed distances $\mathbf{r}$. The
$Z_{p}$ should be well defined over a larger range because of extended
self-similarity~\citep{benzi-et-al:93} and Eq.~\ref{eq:dubrulle}
should remain formally valid for generalized structure
functions $\tilde{S}_{p}(r)$, computed from mass-weighted velocities
$\mathbf{v} \equiv \rho^{1/3}\mathbf{u}$~\citep{kritsuk-et-al:07,
  kritsuk-et-al2:07}.

Several special cases of the model that differ in their parameter
values exist in the literature~\citep[see e.g. the review
by][]{she-zhang:09}. The original model by~\citet{she-leveque:94}
applies most successfully to incompressible turbulence with 1D vortex
filaments as most dissipative structures ($D=1$) and parameter values
$\beta = \Delta = 2/3$. For highly compressible turbulence, parameter
values remain debated.  \citet{boldyrev:02} argues that the most
dissipative structures are 2D shocks, thus $D=2$, and chose to keep
$\Delta=2/3$ and set $\beta=1/3$. By contrast,
\citet{schmidt-et-al:08} argue that $\Delta=1$ (implying $\beta=0$)
to be consistent with Burgers turbulence. A few studies used 3D
simulation data, derived sets of $Z_{p}$, and attempted simultaneous
fits of $\beta$ and $\Delta$~\citep{kritsuk-et-al2:07,
  schmidt-et-al:08,schmidt-et-al:09, folini-et-al:14}. The results
  are inconclusive in that fits of similar quality are obtained for
widely different $\beta$-$\Delta$-pairs. Also using 3D simulation data
(1024$^{3}$, Mach 6) but working with density-weighted moments of the
dissipation rate, \citet{pan-et-al:09} simultaneously fitted $\Delta$
and $D$ to their data. They find $\Delta \in [0.67,0.78]$ and $D \in
[2.04,1.60]$, with a best estimate of $\Delta=0.71$ and $D=1.9$, thus
$\beta=0.35$. The range for $\Delta$ is not compatible with the
suggested $\Delta=1$ (see above), and also $\Delta=2/3$ lies only at
the lower-most bound of the inferred range. Both $\Delta$ and
$\beta$ may thus deviate from their incompressible values ($\beta =
\Delta = 2/3$) as the Mach number increases, making simultaneous
determination of $\beta$ and $\Delta$ a must.

The present study is motivated by this still inconclusive situation.
We want to better understand what factors (accuracy / order of
$Z_{p}$; mildly versus highly compressible turbulence) are most
relevant for the fit quality and why widely different
$\beta$-$\Delta$-pairs yield fits of similar quality. We use this
insight to formulate quantitative estimates of what is needed to
obtain 1\% accurate estimates of $\beta$ and $\Delta$. We present
results in Sect.~\ref{sec:results}, discuss them in
Sect.~\ref{sec:discussion}, and conclude in Sect.~\ref{sec:conc}.
\section{Results}
\label{sec:results}
We first show that $\beta$ and $\Delta$ can be uniquely determined from an
associated (i.e. computed via Eq.~\ref{eq:dubrulle}) pair
$Z_{p_{\mathrm{1}}}$ and $Z_{p_{\mathrm{2}}}$. We then illustrate how
uncertainties in $Z_{p}$ map onto the $\beta$-$\Delta$-plane. Finally,
we give estimates on how accurate the $Z_{p}$ have to be to achieve a
desired accuracy of $\beta$, $\Delta$, and $C$.
\subsection{$\beta$ and $\Delta$  from exact $Z_{p}$}
\label{sec:exact_zp_2p}
Consider two values $Z_{p_{\mathrm{1}}}$ and $Z_{p_{\mathrm{2}}}$ that
both fulfill Eq.~\ref{eq:dubrulle} for the same values $(\beta,
\Delta)$.  In the following, we show that $(\beta, \Delta)$ can
unambiguously (uniquely) be recovered from $Z_{p_{\mathrm{1}}}$ and
$Z_{p_{\mathrm{2}}}$.

We start by rewriting Eq.~\ref{eq:dubrulle}, factoring out $\Delta$:
\begin{equation}
0 = \Delta \left ( \frac{3(1-\beta^{p/3}) - p(1-\beta)}{3(1-\beta)} \right ) +
           \left ( \frac{p - 3 Z_{p_{j}}}{3}\right ), \quad j=1,2.
\label{eq:dubrulle_recast}
\end{equation}
Using $Z_{p_{\mathrm{1}}}$, we can obtain an expression for $\Delta$,
\begin{equation}
\Delta = \left ( \frac{3 Z_{p_{\mathrm{1}}} - p_{\mathrm{1}}}{3} \right )
         \left ( \frac{3(1-\beta)}{3(1-\beta^{p_{\mathrm{1}}/3}) - p_{\mathrm{1}}(1-\beta)} \right ).
\label{eq:dubrulle_delta}
\end{equation}
By now writing Eq.~\ref{eq:dubrulle_recast} with
$Z_{p}=Z_{p_{\mathrm{2}}}$, using Eq.~\ref{eq:dubrulle_delta} to
replace $\Delta$, do some re-ordering of terms, and abbreviating
$p_{\mathrm{1}}/3 \equiv a$ and $p_{\mathrm{2}}/3 \equiv b$, we end up
with the following equation for $\beta$:
\begin{equation}
0 = \left ( \frac{-\beta^{b} + b \beta + 1 - b}{-\beta^{a} + a \beta + 1 - a}  \right ) -
    \left ( \frac{Z_{p_{\mathrm{2}}} - b}{Z_{p_{\mathrm{1}}} - a} \right )
  \equiv \frac{P_{\beta,b}}{P_{\beta,a}} - R_{a,b},
\label{eq:dubrulle_beta}
\end{equation}
or
\begin{equation}
R_{a,b} P_{\beta,a} = P_{\beta,b}.
\label{eq:dubrulle_beta_short}
\end{equation}
The polynomial $P_{\beta,x} \equiv -\beta^{x} + x \beta + 1 - x$, with
$x>0$ and $\beta \in (0,1)$, is a monotonically decreasing
(increasing) function for $x<1$ ($x>1$), as can be seen by taking the
derivative of $P_{\beta,x}$ with respect to $\beta$ and as illustrated
in Fig.~\ref{fig:poly_ratio}.  Consequently,
Eq.~\ref{eq:dubrulle_beta_short} has a unique solution, $\beta$, from
which $\Delta$ can be recovered via Eq.~\ref{eq:dubrulle_delta}. Thus
Eq.~\ref{eq:dubrulle} defines an exact one-to-one correspondence
between pairs $(Z_{p_{\mathrm{1}}}, Z_{p_{\mathrm{2}}})$ and $(\beta,
\Delta)$.
\begin{figure}[tp]
\centerline{
\includegraphics[width=8.5cm,height=5.5cm]{folini_walder_2016_fig_01.eps}
}
\caption{Ratio of polynomials $P_{\beta,b} / P_{\beta,a}$,
  Eq.~\ref{eq:dubrulle_beta}, (y-axis, shown as logarithm) for
  selected exponents $a$ and $b$ as function of $\beta$ (x-axis).
  Colors indicate $b/a$ = 2 (red), 4 (green), 5 (blue), 6 (magenta), 7
  (cyan), all with $a=p_{\mathrm{1}}/3=1/3$, as well as $b/a$ = 7/6 (black)
  with $a=p_{\mathrm{1}}/3=6/3$.}
\label{fig:poly_ratio}
\end{figure}

Two more points deserve to be highlighted, with the help of
Fig.~\ref{fig:poly_ratio}. The ratio $P_{\beta,b} /
P_{\beta,a} = R_{a,b}$ is shown as a function of $\beta$ for different $a$ and
$b$ or, equivalently, $p_{\mathrm{1}}$ and $p_{\mathrm{2}}$. From the
figure it can be taken that, first, largely different values of
$p_{\mathrm{1}}$ and $p_{\mathrm{2}}$ are advantageous since they result
in stronger stratification of $\beta$ with respect to $R_{a,b} =
(Z_{p_{\mathrm{2}}} - b) / (Z_{p_{\mathrm{1}}} - a)$. The cyan curve
in Fig.~\ref{fig:poly_ratio}, which represents $p_{\mathrm{1}} /
p_{\mathrm{2}} = 1/7$, covers a wider range of values on the y-axis
than the black curve ($p_{\mathrm{1}} / p_{\mathrm{2}} = 6/7$).
Secondly, the stratification is stronger for small $\beta$. Somewhat
anticipating Sect.~\ref{sec:perturbed_zp_2p}, we thus expect
uncertainties in the $Z_{p}$ to be less important if $Z_{p}$ are
available for largely different $p$ and if they are associated with
(yet to be determined) small values of $\beta$.
\subsection{Uncertainty of  $Z_{p}$ in the $\beta$-$\Delta$-plane}
\label{sec:perturbed_zp_2p}
\subsubsection{Single  $Z_{p}$}
\label{sec:zp_in_beta_delta}
\begin{figure}[tbp]
\centerline{
\includegraphics[width=5.1cm]{folini_walder_2016_fig_02a.eps}
\hspace{-0.7cm}
\includegraphics[width=5.1cm]{folini_walder_2016_fig_02b.eps}
\hspace{-0.3cm}
}
\centerline{
\includegraphics[width=5.1cm]{folini_walder_2016_fig_02c.eps}
\hspace{-0.7cm}
\includegraphics[width=5.1cm]{folini_walder_2016_fig_02d.eps}
\hspace{-0.3cm}
}
\caption{$Z_{p}$ in the $\beta$-$\Delta$-plane, role of $p$. Shown is
  $Z_{p}$ (color coded) for $p=1$ (top left), $p=2$ (top right), $p=4$ (bottom left)
  and $p=10$ (bottom right). For $\beta=1/3$ and $\Delta=2/3$ (red dot),
  the curve of constant $Z_{p}$ (black) is shown, as well as curves
  of  $\pm 5$\% different $Z_{p}$ (white).}
\label{fig:zp_of_beta_delta}
\end{figure}
From Eq.~\ref{eq:dubrulle_delta} it is clear that each $Z_{p}$ defines
a curve in the $\beta$-$\Delta$-plane. If $Z_{p}$ is derived from
model data or observations, it will typically come with an uncertainty
estimate, e.g. $\delta Z_{p}/Z_{p} \le 5$\%, with $\delta$ indicating
the uncertainty. In the $\beta$-$\Delta$-plane, this uncertainty range
translates into an area around the $Z_{p}$ curve. An illustration is
given in Fig.~\ref{fig:zp_of_beta_delta}. The following points may be
made.

One value of $Z_{p}$ (a line of constant $Z_{p}$ in the
$\beta$-$\Delta$-plane) is compatible with a (large) range of $\beta$
and / or $\Delta$ that always includes $\Delta=1$ and $\beta=0$.  The
range tends to be smaller for $Z_{p}$ associated with small $\beta$
and large $\Delta$ (i.e. the lower right corner of
$\beta$-$\Delta$-plane). Uncertainties associated with $Z_{p}$ (5\% in
Fig.~\ref{fig:zp_of_beta_delta}, white curves) augment the range,
especially for $p=2$ and $p=4$, as well as for small $\Delta$ and large
$\beta$ (top left corner of the plane).  Also apparent from
Fig.~\ref{fig:zp_of_beta_delta} (or from taking the derivative with
respect to $\Delta$ of Eq.~\ref{eq:dubrulle}): for fixed $\beta$ and
$p<3$ ($p>3$), $Z_{p}$ is a monotonically increasing (decreasing)
function of $\Delta$. A similar statement holds for $Z_{p}$ as a
function of $\beta$ for fixed $\Delta$.

In summary, we expect uncertainties in the $Z_{p}$ to be more of
an issue if only low orders of $p$ (up to about 4) are available
and / or if the (yet to be determined) $\beta$ is large.
\subsubsection{Multiple $Z_{p}$}
\label{sec:multiple_zp}
We now turn to multiple $Z_{p}$ and their associated uncertainty
ranges $\delta Z_{p}$, and ask what area they define in the
$\beta$-$\Delta$-plane. An illustration is given in
Fig.~\ref{fig:dmap}. Starting from one specific pair of $\beta=1/3$
and $\Delta=2/3$ and computing $Z_{p}$ for $p=1,.\,.\,.\,,5$ (left) or $p=1,.\,.\,.\,,7$
(right), we show pairs of 5\% perturbed $Z_{p}$ curves,
i.e. $1.05 \cdot Z_{p}$ and $0.95 \cdot Z_{p}$.

As can be seen, only a small fraction of the $\beta$-$\Delta$-plane
lies between all pairs of perturbed curves. Yet this area comprises a
wide range of $(\beta,\Delta)$ values or co-dimensions.  The 5\%
uncertainty in the $Z_{p}$ translates into a much larger uncertainty
(in per cent) for $\beta$ and $\Delta$. Closer inspection reveals that
the area is actually defined by only two sets of curves: those for
$p=1$ and $p=5$ (left panel) or $p=7$ (right panel). The latter area
is smaller, which indicates that higher order structure functions
constrain the problem of finding $\beta$ and $\Delta$ from a set of
$Z_{p}$ more strongly. Also apparent from Fig.~\ref{fig:dmap} is the
dominant role of the $p=1$ curve for narrowing down the composite area
between all curves.  All this is in line with the expectation (see
Sects.~\ref{sec:exact_zp_2p}) that $Z_{p}$ for largely different $p$
are advantageous for the determination of $\beta$ and $\Delta$.
\begin{figure}[tbp]
\centerline{
\includegraphics[width=5.2cm]{folini_walder_2016_fig_03a.eps}
\hspace{-0.7cm}
\includegraphics[width=5.2cm]{folini_walder_2016_fig_03b.eps}
}
\caption{Part of $\beta$-$\Delta$-plane (color coded in co-dimension
  $C=\Delta / (1-\beta)$, $0<C<3$) within joint reach of (at most) 5\%
  perturbed $Z_{p}$, starting from $Z_{p}$ for $\Delta=2/3$ and
  $\beta=1/3$ (red dot). Left panel: Curves of 5\% perturbed
  $Z_{p}$ values for $p=1$ (white), $p=2$ (blue), $p=4$ (red), and
  $p=5$ (green) and part of the $\beta$-$\Delta$-plane enclosed by all
  of them. Right panel: Same as left panel but also including
  $p=6$ (purple) and $p=7$ (black).}
\label{fig:dmap}
\end{figure}

The relevance of the overall location in the $\beta$-$\Delta$-plane is
illustrated in Fig.~\ref{fig:dmap_var}. Again, the area shown is
contained within 5\% perturbed $Z_{p}$ curves for two additional
$(\beta, \Delta)$ pairs.  As can be seen, smaller values of $\beta$
(lower panels) result in smaller areas, independent of $\Delta$. The
crucial role of the $p=1$ (white) and $p=5$ (green) curves for
confining the area persists.  Table~\ref{tab:delta_c} gives a
quantitative idea of the relevance of $\beta$, $\Delta$, $\delta
Z_{p}$, and $p_{\mathrm{max}}$ for the uncertainty range $\pm \delta
C$ of the co-dimensions $C$. A small $\delta C$ basically requires a
small $\delta Z_{p}$, a large $p_{\mathrm{max}}$, a small $\beta$, and
a large $\Delta$. The concrete numbers highlight the difficulty (or
ill-posedness) of the problem. The situation is
worse for larger $\beta$ (bottom rows in Table~\ref{tab:delta_c}) and
better for smaller values of $\beta$ (not shown).

We emphasize that the above considerations serve only as illustration. We
looked at the area confined by a set of $Z_{p} \pm \delta Z_{p}$
curves. We have not yet considered the problem of estimating best-fit
$\beta_{\mathrm{f}}$, $\Delta_{\mathrm{f}}$, and thus $C_{\mathrm{f}}$
for a set of given $Z_{p}$. Such a best-fit solution may lie outside
the area considered here.

In summary, very accurate $Z_{p}$ are needed to derive reliable best
estimates for $\beta$, $\Delta$, and $C,$ and smaller values of $\beta$
help.
\subsection{Best-fit $\beta_{\mathrm{f}}$ and $\Delta_{\mathrm{f}}$
  from uncertain $Z_{p}$}
\label{sec:best_fit_beta_delta}
We now turn to our actual problem of interest: given a set of
perturbed (uncertain) $\tilde{Z}_{p} = Z_{p} + \delta Z_{p}$, what are
associated best-fit estimates for $\beta_{\mathrm{f}}$ and
$\Delta_{\mathrm{f}}$? Different techniques exist to cope with this kind of
question~\citep[e.g.][]{najm:09, lemaitre-knio:10}. We use a simple
Monte Carlo approach.

We start with a pair $(\beta, \Delta)$ and a maximum order
$p_{\mathrm{max}}$, then use Eq.~\ref{eq:dubrulle} to obtain a set of
$Z_{p} = Z_{p} (\beta, \Delta)$ for $p=1$ to $p=p_{\mathrm{max}}$.
Each of these $Z_{p}$ we perturb randomly (uniformly distributed
random numbers) by, at most, $\alpha$\%, which gives us a perturbed set
of $\tilde{Z}_{p}$. For this set of $\tilde{Z}_{p}$ we then seek to
find best-fit $\beta_{\mathrm{f}}$ and $\Delta_{\mathrm{f}}$. In the
following, we do not consider one set of $\tilde{Z}_{p}$, as would be
the case in a real application (unless multiple time slices are
available, see Sect.~\ref{sec:discussion}). Instead, we take a
statistical view for the problem by looking at a large number (1000 to
100\,000, see below) of randomly generated sets of perturbed
$\tilde{Z}_{p}$. This enables us, in a statistical sense, to relate
the accuracy of the $\tilde{Z}_{p}$ with the accuracy of the fitted
parameters. Our approach leaves us with two free parameters, the
uncertainty $\alpha$ and the maximum order $p_{\mathrm{max}}$.

\begin{table}[ht]
  \caption{Illustration of range $\delta C$ of co-dimension $C$ for given order $p$ and
    uncertainty $\delta Z_{p}$ of structure functions for two $\beta$-$\Delta$ pairs.}
\begin{center}
\begin{tabular}{cc}
\hline
\multicolumn{2}{c}{$\beta = 1/3$, $\Delta = 2/3$, $C=1$} \\
\rule[-2mm]{0pt}{6mm}      $\delta Z_{p}/Z_{p} = 1$\%    &     $p = 1,.\,.\,.\,,9 $ \\
\hline
$p = 1,.\,.\,.\,,5$:\hspace{0.3cm} $0.79-1.46$ \hspace{0.2cm} & \hspace{0.2cm} $\delta Z_{p}/Z_{p} = 1$\%:\hspace{0.3cm} $0.93-1.08$\\
$p = 1,.\,.\,.\,,6$:\hspace{0.3cm} $0.83-1.23$ \hspace{0.2cm} & \hspace{0.2cm} $\delta Z_{p}/Z_{p} = 2$\%:\hspace{0.3cm} $0.84-1.27$\\
$p = 1,.\,.\,.\,,7$:\hspace{0.3cm} $0.87-1.18$ \hspace{0.2cm} & \hspace{0.2cm} $\delta Z_{p}/Z_{p} = 3$\%:\hspace{0.3cm} $0.74-1.52$\\
$p = 1,.\,.\,.\,,8$:\hspace{0.3cm} $0.90-1.13$ \hspace{0.2cm} & \hspace{0.2cm} $\delta Z_{p}/Z_{p} = 4$\%:\hspace{0.3cm} $0.70-1.92$\\
$p = 1,.\,.\,.\,,9$:\hspace{0.3cm} $0.93-1.08$ \hspace{0.2cm} & \hspace{0.2cm} $\delta Z_{p}/Z_{p} = 5$\%:\hspace{0.3cm} $0.65-2.25$\\
\hline
\multicolumn{2}{c}{$\beta = 2/3$, $\Delta = 2/3$, $C=2$} \\
\rule[-2mm]{0pt}{6mm}      $\delta Z_{p}/Z_{p} = 1$\%    &       $p = 1,.\,.\,.\,,9 $ \\
\hline
$p = 1,.\,.\,.\,,5$:\hspace{0.3cm} $0.51-3.00$ \hspace{0.2cm} & \hspace{0.2cm} $\delta Z_{p}/Z_{p} = 1$\%:\hspace{0.3cm} $1.11-3.00$\\
$p = 1,.\,.\,.\,,6$:\hspace{0.3cm} $0.70-3.00$ \hspace{0.2cm} & \hspace{0.2cm} $\delta Z_{p}/Z_{p} = 2$\%:\hspace{0.3cm} $0.71-3.00$\\
$p = 1,.\,.\,.\,,7$:\hspace{0.3cm} $0.86-3.00$ \hspace{0.2cm} & \hspace{0.2cm} $\delta Z_{p}/Z_{p} = 3$\%:\hspace{0.3cm} $0.55-3.00$\\
$p = 1,.\,.\,.\,,8$:\hspace{0.3cm} $1.00-3.00$ \hspace{0.2cm} & \hspace{0.2cm} $\delta Z_{p}/Z_{p} = 4$\%:\hspace{0.3cm} $0.44-3.00$\\
$p = 1,.\,.\,.\,,9$:\hspace{0.3cm} $1.11-3.00$ \hspace{0.2cm} & \hspace{0.2cm} $\delta Z_{p}/Z_{p} = 5$\%:\hspace{0.3cm} $0.36-3.00$\\
\hline
\end{tabular}
\end{center}
\label{tab:delta_c}
\end{table}

\begin{figure}[tbp]
\centerline{
\includegraphics[width=5.2cm]{folini_walder_2016_fig_04a.eps}
\hspace{-0.7cm}
\includegraphics[width=5.2cm]{folini_walder_2016_fig_04b.eps}
}
\centerline{
\includegraphics[width=5.2cm]{folini_walder_2016_fig_04c.eps}
\hspace{-0.7cm}
\includegraphics[width=5.2cm]{folini_walder_2016_fig_04d.eps}
}
\caption{Same as Fig.~\ref{fig:dmap}, left, but for $\beta=\Delta=2/3$
  (top left), $\beta=\Delta=1/3$ (top right),
  $\beta=0.048$ and $\Delta=1/3$ (bottom left), as well as and
  $\beta=0.048$ and $\Delta=2/3$ (bottom right).}
\label{fig:dmap_var}
\end{figure}
\subsubsection{Minimization of least square error in $Z\_{p}$ }
\label{sec:min_least_square}
\begin{figure}[tbp]
\centerline{
\includegraphics[width=4.6cm,height=4.2cm]{folini_walder_2016_fig_05a.eps}
\includegraphics[width=4.3cm,height=4.2cm]{folini_walder_2016_fig_05b.eps}
}
\caption{Best-fit $\beta_{\mathrm{f}}$ and $\Delta_{\mathrm{f}}$ from
  5\% perturbed $Z_{p}$.  Left: 2D histogram (contours, log10,
  spacing 0.5, spanning three orders of magnitude) of best-fit
  $\beta$-$\Delta$-values from 100\,000 perturbed data sets.  We note
  that the 2D histogram shows two peak values (indicated by cyan
  dots), none of them co-located with the unperturbed $(\beta=1/3,
  \Delta=2/3)$ pair (red dot). Right: Underestimation of $Z_{5}$
  favors $\Delta_{\mathrm{f}}=1$. Shown is, for a subset of 1000
  perturbed data sets, $\Delta_{\mathrm{f}}$ as function of
  $\tilde{Z}_{5}/Z_{5}$, again for an unperturbed pair $(\beta=1/3,
  \Delta=2/3)$. Colors indicate $\tilde{Z}_{1}/Z_{1}$.}
\label{fig:stat_map}
\end{figure}
A straightforward way to determine best-fit $\beta_{\mathrm{f}}$ and
$\Delta_{\mathrm{f}}$ for any given set of $\tilde{Z}_{p}$,
$p=p_{\mathrm{min}}.p_{\mathrm{max}}$ is to minimize
\begin{equation}
  \sum_{p=1}^{p_{\mathrm{max}}}
       [\tilde{Z}_{p} - Z_{p}(\beta_{\mathrm{f}},\Delta_{\mathrm{f}})]^{2}
\label{eq:to_minimize}
\end{equation}
over the $\beta$-$\Delta$-plane. To find the minimum, we compare the
$\tilde{Z}_{p}$ with pre-computed values $Z_{p}(\beta , \Delta)$ on a
fine $\beta$-$\Delta$-grid ($\beta , \Delta \in (0,1)$; grid-spacing
0.002). The associated co-dimension is given by $C_{\mathrm{f}} =
\Delta_{\mathrm{f}} / (1 - \beta_{\mathrm{f}})$.

To capture the range of potential outcomes for a range of similarly
perturbed data sets $\tilde{Z}_{p}$, we produced 100\,000 perturbed
data sets, for each of which we determined $\beta_{\mathrm{f}}$ and
$\Delta_{\mathrm{f}}$. For initial values $(\beta, \Delta) = (1/3,
2/3)$ and (at most) 5\% perturbed $Z_{p}$ for $p=1,.\,.\,.\,,5$, the result is
summarized in Fig.~\ref{fig:stat_map}.

Shown in the left panel of Fig.~\ref{fig:stat_map} is a 2D histogram
(contours) of our 100\,000 best-fit $(\beta_{\mathrm{f}},
\Delta_{\mathrm{f}})$ pairs. Two points are noteworthy. First, the
overall area defined by the histogram is similar to the area in
Fig.~\ref{fig:dmap}, left panel. This is remarkable since the area in
Fig.~\ref{fig:dmap} is strictly defined by the 5\% uncertainty of the
$Z_{p}$, whereas the area in Fig.~\ref{fig:stat_map} is defined
through a minimization problem. Second, the 2D histogram has an interior
structure with two peaks, around $(\beta, \Delta)=(0.1, 0.4)$ or
$C=0.4$ and $(\beta, \Delta)=(0.45, 1.0)$ or $C=1.8$ (cyan dots). None
of them is co-located with the initial, unperturbed $(\beta, \Delta) =
(1/3, 2/3)$ pair (red dot, $C=1$).

Three questions come to mind. Where do the two peaks in the 2D
histogram come from? Do other $(\beta, \Delta)$ pairs result in a
qualitatively different picture? Can the minimization procedure be
improved to better recover the initial, unperturbed $(\beta, \Delta)$
pair? We address the first two questions in the following while
postponing the third question for Sect.~\ref{sec:alt_best_fits}.

The existence and location of the two peaks can be understood, at
least qualitatively, from two observations.  First, minimization via
Eq.~\ref{eq:to_minimize} gives more weight to larger $p$, as they are
associated with larger values of $Z_{p}$. Roughly speaking, the
best-fit $(\beta_{\mathrm{f}}, \Delta_{\mathrm{f}})$ pair tends to lie
on or close to the curve defined by $\tilde{Z}_{5}$. Moving away from
that curve results in a large penalty in the form of a large
contribution to the sum in Eq.~\ref{eq:to_minimize}. Second, this
translates the minimization problem into the question of where the
curves for $p<5$ come closest to the curve defined by $\tilde{Z}_{5}$.
For illustration, we consider two extreme values of $\tilde{Z}_{5}$.
To stay on the lower green curve in the left panel of
Fig.~\ref{fig:dmap}, ($\tilde{Z}_{5} = 0.95 Z_{5}$) and, at the same
time, be as close as possible to any of the white curves ($p=1$)
results in a $(\beta_{\mathrm{f}}, \Delta_{\mathrm{f}})$ pair to the
right, at $\Delta_{\mathrm{f}}=1$.  By contrast, the upper green curve
(105\% of the exact $Z_{5}$ curve) only comes closest to (intersects)
any white curve between 95\% $Z_{1}$ and 105\% $Z_{1}$ in a region
further to the left. Clearly, the full problem is more intricate, with
also curves for $Z_{2}$ and $Z_{4},$ and the $Z_{5}$ curve not
necessarily adopting one of its two extreme values. Nevertheless,
Fig.~\ref{fig:stat_map}, right panel, suggests the full data to be in
line with the above reasoning. For 1000 randomly picked data
sets from the left panel, we show $\Delta_{\mathrm{f}}$ as a function of
$\tilde{Z}_{5} / Z_{5}$, with $Z_{5}$ the exact value. Colors indicate
$\tilde{Z}_{1} / Z_{1}$.  As can be seen, $\Delta_{\mathrm{f}}=1$
indeed tends to be associated with small $\tilde{Z}_{5}$ and small
$\tilde{Z}_{1}$ (lower green and upper white curve in
Fig.~\ref{fig:dmap}, left panel). Particularly low values of
$\Delta_{\mathrm{f}}$ (e.g. $\Delta_{\mathrm{f}} \approx 0.4$) tend to
occur for large $\tilde{Z}_{5}$ and any $\tilde{Z}_{1}$ (upper green
curve and any white curve in Fig.~\ref{fig:dmap}, left panel).

Concerning other initial values (other exact $(\beta, \Delta)$ pairs),
a similar situation arises in the sense that double peaked histograms
emerge. Details depend, however, on the concrete values of $\beta$ and
$\Delta$, on the assumed uncertainty (5\% or more / less), and on
$p_{\mathrm{max}}$. An illustration is given in
Fig.~\ref{fig:accu_vs_codim}, by means of 1D histograms of
$C_{\mathrm{f}} = \Delta_{\mathrm{f}} / (1 - \beta_{\mathrm{f}})$.
These 1D histograms are less intricate than the 2D histogram in
Fig.~\ref{fig:stat_map}, left, yet still capture the essentials.
We show histograms for $\Delta=2/3$ and different values of $\beta
\in [0.0476, 2/3]$ (corresponding to $C \in [0.7, 2]$) and accuracies
between 0.1\% and 20\%. Five points may be made.  First, the
double-peaked structure that is apparent in the $\beta$-$\Delta$ plane in
Fig.~\ref{fig:stat_map} re-appears as a double peak in the 1D
$C_{\mathrm{f}}$-histograms of Fig.~\ref{fig:accu_vs_codim} (panel in
row three, column three). Second, the double-peak vanishes as $\beta$
and the uncertainty both become small (lower left corner of the
figure). For the same uncertainty, the double-peak exists for large
$\beta$ but not for small $\beta$ (third row in
Fig.~\ref{fig:accu_vs_codim}). For the same $\beta$, the double-peak
exists for large uncertainties but not for small ones (second column
of Fig.~\ref{fig:accu_vs_codim}). Third, going to really small values
of $\beta$ and the uncertainty, the histogram becomes symmetric with
one central peak.  Fourth, only for these really small values is the
co-dimension of the initially prescribed $(\beta, \Delta)$ pair (show
in red) co-located with the peak of the histogram. Fifth, for
$\beta=\Delta=2/3$ (right column) the histogram peaks at
$C_{\mathrm{f}}=3$ instead of $C=2$, unless the accuracy is really
high (0.1\%, last row). This is understandable from the arguments
presented above, with regards to the origin of the two peaks in the 2D
histogram in Fig.~\ref{fig:stat_map}, and from looking at the green
($p=5$) and white ($p=1$) curves in Fig.~\ref{fig:dmap_var}, upper
left panel.

In summary, unless both, $\delta Z_{p}$ and $\beta$ are small, best-fit
values will preferentially reside in either one of the two peaks of
the histograms in Figs.~\ref{fig:stat_map} or~\ref{fig:accu_vs_codim}
instead of merely scattering around the correct solution, as in
Fig.~\ref{fig:accu_vs_codim}, lower left panel.
\begin{figure}[tp]
\centerline{
\includegraphics[width=8.8cm,height=10cm]{folini_walder_2016_fig_06.eps}
}
\caption{Role of $\beta$ (columns) and accuracy (rows) of associated
  $Z_{p}$, for fixed $\Delta = 2/3$. Shown are PDFs (y-axis) of
  $C_{\mathrm{f}}$ (x-axis), 1000 random data sets, powers $p=1,.\,.\,.\,,5$
  for $\beta$=0.0476 (first column), $\beta$=0.17 (second
    column), $\beta$=1/3 (third column), and $\beta=2/3$
    (fourth column). Corresponding exact co-dimensions (red lines)
  are, from left to right: $C=0.7$, $C=0.8$, $C=1$, and $C=2$.
  Individual rows from top to bottom contain accuracies of 20\%,
  10\%, 5\%, 2\%, 0.5\%, and 0.1\%. As can be seen, the larger $\beta,$
  the more severe are the consequences of inaccuracies in the $Z_{p}$.
  Histograms in the upper right (large $\beta$, low accuracy of the
  $Z_{p}$) look worst.  We note that axis ranges differ among panels,
  to best capture the shape of each histogram.}
\label{fig:accu_vs_codim}
\end{figure}
\subsubsection{Alternative ways to obtain best-fit  $\beta_{\mathrm{f}}$ and $\Delta_{\mathrm{f}}$}
\label{sec:alt_best_fits}
A number of ideas come to mind on how one may improve the best-fit
approach detailed in Sect.~\ref{sec:min_least_square}.

Recalling the findings in Sect.~\ref{sec:zp_in_beta_delta}, including
higher values of $p$ in the best-fit estimate should improve the
situation. From Fig.~\ref{fig:d_b_of_p} it can be taken that this is
indeed the case, at least for the example shown ($\beta=1/3$,
$\Delta=2/3$, accuracy of 5\%). However, the improvement may be
regarded as rather modest. Going from $p=5$ to $p=10$, as is
illustrated in the figure, has about the same effect as staying with
$p=5$ but going from an accuracy of 5\% to an accuracy of 2\%.  From a
practical point of view it also seems questionable whether high order
structure functions can meet the accuracy requirements.  In numerical
simulations, higher order structure functions are probably more prone
to the bottleneck effect~\citep{dobler-et-al:03, kritsuk-et-al:07}.

Another way to improve the situation could be to go to weighted root
mean squares instead of the unweighted sum in
Eq.~\ref{eq:to_minimize}. Hopefully this breaks the dominant
role of the highest order $Z_{p}$ available (see
Sect.~\ref{sec:min_least_square}) and, ultimately, leads to more
accurate best-fit $\beta_{\mathrm{f}}$ and $\Delta_{\mathrm{f}}$.  Two
weightings come to mind.  On the one hand, weights proportional to the
inverse of the $Z_{p}$ with the goal of giving equal weight to each
term in the sum, thus reducing the "overweight" of larger $p$ in the
sum.  On the other hand, we could try to give more weight to $p$ terms
with a higher accuracy (smaller $\delta Z_{p}$).  Corresponding
information may be available, e.g. from the numerical determination of
the $Z_{p}$. We tried both ideas but neither choice of weights
decidedly improved the best-fit values. Weighting tends to change the
relative height of the two peaks in the double peaked histograms of
Fig.~\ref{fig:accu_vs_codim}, but it does not get rid of the double
peaked structure.

We interpret this finding in the following way. First, there are
likely always several $Z_{p}$ that do not have their exact values and
thus draw the solution in different directions, away from its exact
value. Second, the different curve shapes are important so that, even for
weighted sums, the terms $p=1$ and $p=p_{\mathrm{max}}$ are of crucial
importance for the overall fit.

In summary, none of the above alternative ways of fitting
simultaneously for $\beta$ and $\Delta$ provides clearly superior
results to what can be obtained from the straightforward minimization
of Eq.~\ref{eq:to_minimize}. We conclude that, for successful two-parameter
fits of $\beta$ and $\Delta$, highly accurate $Z_{p}$ are a must. A
quantitative estimate of "highly accurate" is given in the next
section.
\begin{figure}[tbp]
\centerline{
\includegraphics[width=8.8cm]{folini_walder_2016_fig_07.eps}
}
\caption{Role of higher order moments. Inclusion of higher
  order ESS scaling exponents (from $p=5$, top left, to $p=10$, bottom right)
  gradually reduces the erroneous peak in the best-fit co-dimension
  (around $C_{\mathrm{f}}=2.5$). Shown are PDFs of $C_{\mathrm{f}}$
  for at most 5\% perturbed $Z_{p}$ values (100\,000 random data sets)
  and exact pair $(\beta=1/3, \Delta=2/3)$. Spacing of the
  $\beta$-$\Delta$-grid for fitting is 0.002.}
\label{fig:d_b_of_p}
\end{figure}
\subsubsection{Required accuracy of $Z_{p}$ for 'good' best-fit  $\beta_{\mathrm{f}}$ and $\Delta_{\mathrm{f}}$}
\label{sec:accu_needed}
We now ask how accurate the $Z_{p}$ have to be in order to reach a
prescribed accuracy of $C_{\mathrm{f}} = \Delta_{\mathrm{f}} /
(1-\beta_{\mathrm{f}})$ via fitting $\beta_{\mathrm{f}}$ and
$\Delta_{\mathrm{f}}$.

We formulate our accuracy goal in terms of only $C_{\mathrm{f}}$, since we
illustrated in Sect.~\ref{sec:min_least_square} that a single
peaked and roughly symmetric distribution of $C_{\mathrm{f}}$ goes
hand in hand with high accuracy, not only of $C_{\mathrm{f}}$ but also
of the underlying two parameter fit, $\beta_{\mathrm{f}}$ and
$\Delta_{\mathrm{f}}$. If the latter is not accurate
enough, a double peaked distribution for $C_{\mathrm{f}}$ results. We
find, as a rule of thumb, a single peak distribution if 2/3 of all
$C_{\mathrm{f}}$ lie within 10\% or better of the exact $C$. We use
Eq.~\ref{eq:to_minimize} for the two parameter fit, as the more
elaborate attempts of Sect.~\ref{sec:alt_best_fits} gave no decidedly
better results. As theoretical arguments suggest $\Delta = 2/3$ or
larger~\citep{dubrulle:94, schmidt-et-al:08}, we concentrate on that
part of the $\beta$-$\Delta$-plane.

In practical terms, we define a grid of exact pairs $(\beta, \Delta)$
via a (nearly) equidistant grid of $C = 0.4,.\,.\,.\,,2$ and $\Delta =
0.7,.\,.\,.\,,0.99$ plus, in addition, $\Delta=2/3$. We equally define some
fixed levels of perturbations: $\delta Z_{p}/Z_{p}$ (in \%) $\in
[0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 15, 20]$. For each exact pair and
each perturbation, we created 1000 perturbed sets of $\tilde{Z}_{p}$,
$p=1,.\,.\,.\,,5$ (see Sect.~\ref{sec:best_fit_beta_delta}). Each perturbed set
is fitted via Eq.~\ref{eq:to_minimize}. For each initial
pair ($\beta$, $\Delta$) and for each prescribed $\delta Z_{p} /
Z_{p}$ this yields 1000 fitted pairs $(\beta_{\mathrm{f}}, \Delta_{\mathrm{f}})$
and derived co-dimensions $C_{\mathrm{f}}$. These can, in principle,
be arranged in histograms, as in Fig.~\ref{fig:accu_vs_codim}.
Finally, we identify the largest $\delta Z_{p}/Z_{p}$ for which 2/3 of
all $C_{\mathrm{f}}$ lie within the demanded accuracy of the exact,
initial $C$.
\begin{figure}[tbp]
\centerline{\includegraphics[width=8.0cm]{folini_walder_2016_fig_08a.eps}}
\centerline{\includegraphics[width=8.0cm]{folini_walder_2016_fig_08b.eps}}
\caption{Required accuracy (color coding) of perturbed $Z_{p}$,
  $p=1,.\,.\,.\,,5$, such that for at least 2/3 of the best-fit pairs
  ($\beta_{\mathrm{f}}$, $\Delta_{\mathrm{f}}$) the associated
  $C_{\mathrm{f}}$ are within 10\% (top) or 5\% (bottom)
  of the $C$ associated with the initial, unperturbed ($\beta$,
  $\Delta$) pair. As can be seen, the required accuracy depends
  crucially on the location within the $\beta$-$\Delta$ plane. Gray
  lines indicate the range within which 2/3 of the actual
  $\beta_{\mathrm{f}}$ (vertical direction) and $\Delta_{\mathrm{f}}$
  (horizontal direction) lie. Black lines indicate constant $C=1$
  (lower line) and $C=2$ (upper line). For details see text.}
\label{fig:accu_c_bd_plane}
\end{figure}

Fig.~\ref{fig:accu_c_bd_plane} illustrates the result. Obviously, the
accuracy of the $\tilde{Z}_{p}$ (colored squares, in \%) that is
needed to get at least 2/3 of best-fit $C_{\mathrm{f}}$ to lie within
10\% (top panel) or 5\% (bottom panel) of the exact (initial) $C$
depends on the position within the $\beta$-$\Delta$-plane. In the
lower right parts, $\delta Z_{p}/Z_{p} \ge 1 \mathrm{\%}$ is
sufficient to get 10\% accurate $C_{\mathrm{f}}$, and $\delta
Z_{p}/Z_{p} \ge 0.5 \mathrm{\%}$ yields 5\% accurate $C_{\mathrm{f}}$.
By contrast, in the upper left parts of the plane ($\beta \ge 0.5$)
one needs $\delta Z_{p}/Z_{p} \le 0.1 \mathrm{\%}$ to get 10\%
accurate $C_{\mathrm{f}}$.  Gray lines in
Fig.~\ref{fig:accu_c_bd_plane} (for clarity only shown for a subset of
the colored squares) indicate the range within which at least 2/3 of
the actually fitted $\beta_{\mathrm{f}}$ and $\Delta_{\mathrm{f}}$
lie.  The range is larger for $\Delta$ (horizontal lines) than for
$\beta$ (vertical lines). This is plausible from
Fig.~\ref{fig:dmap_var}, from the area confined by multiple $Z_{p}$
curves. Repeating Fig.~\ref{fig:accu_c_bd_plane} but demanding 10\%
accuracy for $\Delta_{\mathrm{f}}$ instead of $C_{\mathrm{f}}$ yields
a similar pattern in the $\beta$-$\Delta$-plane (not shown), while
demanding 10\% accuracy for $\beta_{\mathrm{f}}$ gives a much more
homogeneous pattern (0.2\% to 0.5\% accuracy for $\tilde{Z}_{p}$, not
shown).

For highly compressible turbulence, best-fit $\beta_{\mathrm{f}}$ and
$\Delta_{\mathrm{f}}$ are expected to lie in the lower part of the
$\beta$-$\Delta$-plane in Fig.~\ref{fig:accu_c_bd_plane}, roughly
$\beta \le 1/3$ and $\Delta \ge 2/3$~\citep{boldyrev:02,
  padoan-et-al:04, pan-et-al:09}. Here, accuracies of 0.5\%, 0.1\%,
and 0.05\% for the $\tilde{Z}_{p}$ translate into accuracies of 10\%,
2\%, and 1\% for $\beta_{\mathrm{f}}$ and $\Delta_{\mathrm{f}}$ (not
shown). Fits of similar quality require much more accurate
$\tilde{Z}_{p}$ in the mildly compressible regime, where $\beta > 1/3$
(upper part of panels in Fig.~\ref{fig:accu_c_bd_plane}) and
ultimately $\beta = 2/3$ in the incompressible limit. We note that in
practical applications, best-fit values may be further improved by
combining, for example, data from different time slices~\citep{pan-et-al:09}.

In summary, a 2\% (1\%) accurate simultaneous fit for
$\beta_{\mathrm{f}}$ and $\Delta_{\mathrm{f}}$ should be possible in
the highly compressible regime if the $Z_{p}$ are 0.1\% (0.05\%)
accurate.  If no satisfying fit is possible for such $Z_{p}$, this may
indicate that the model is not applicable to the turbulence data under
examination.
\section{Discussion}
\label{sec:discussion}
We address three topics. First, can the necessary accuracy for
the $Z_{p}$ be met in practical applications? Second, if we had this type of
accurate simulation data, what could be learned about the hierarchical
structure model and its applicability or non-applicability to driven,
isothermal, supersonic turbulence in a 3D periodic box?  Third, we
want to briefly revisit the frequently used one parameter fits.

We start with the question whether 0.1\% or even 0.05\% accurate
$Z_{p}$ for $p=1,.\,.\,.\,,5$ are achievable, as are needed to get 2\% (1\%)
accurate fits for $\beta$ and $\Delta$. The answer is probably yes, at
least in the context of 3D periodic box simulations.
\citet{schmidt-et-al:08} estimate the accuracy of their $Z_{p}$,
$p=1,.\,.\,.\,,5$, to 1\% (3D box simulations, $1024^{3}$).
\citet{kritsuk-et-al:07} estimate 1\% accuracy or better for absolute
scaling exponents $\zeta_{p}$, $p=1,.\,.\,.\,,3$ (3D box simulations,
$1024^{3}$). Meanwhile, 3D box simulations with $4096^{3}$
exist~\citep[e.g.][]{federrath:13, beresnyak:14}. A first order
estimate suggests the four times better resolution to translate into
four times (first order scheme) or 16 times (second order scheme) more
accurate structure functions, thus accuracies of 0.25\% or even
0.0625\%. Moreover, if the accuracy of the $Z_{p}$ is good enough to
avoid double peaked histograms as in Fig.~\ref{fig:accu_vs_codim},
accuracy may be further enhanced by exploring multiple time slices.
\citet{pan-et-al:09} (3D box simulations, $1024^{3}$) used data from
nine time slices for their two parameter fit. Their work is comparable
although they rely on dissipation rates instead of velocity structure
functions, since the involved scaling exponents ($\tau_{p}$ for the
dissipation rate) are structurally similar according to Kolmogorov's
refined similarity hypothesis~\citep{kolmogorov:62}, $\zeta_{p} = p/3
+ \tau_{p/3}$.  The tests we carried out using ratios of $\tau_{p}$
instead of $\zeta_{p}$ indeed show a similar behavior. From the quality
of their fit and based on our results, we estimate their $\tau_{p}$ to
be about 1\% accurate, which is plausible given their numerical
resolution. A reliable two parameter fit to 3D periodic box data of
highly supersonic turbulence that is based on velocity structure
functions thus appears feasible with today's data.

What could be learned from better simulation data ($4096^{3}$ or
better) and associated, more accurate $Z_{\mathrm{p}}$?  Each
$Z_{\mathrm{p}}$ defines a curve with associated uncertainty in the
$\beta$-$\Delta$ plane, the curves for different $p$ may or may not
intersect to within uncertainties. If they intersect, the model
by~\citet{she-leveque:94} may indeed carry over to highly compressible
turbulence. It is then interesting to see whether the fitted range for
$\Delta$, currently estimated as ($0.67-0.78$)
by~\citet{pan-et-al:09}, remains compatible with $\Delta=2/3$, the
value tacitly assumed in a large body of literature. It is also
interesting to check whether $\beta=1/3$, as theoretically anticipated
by~\citet{boldyrev:02}.  If indeed $(\beta, \Delta) = (1/3, 2/3)$,
results from one parameter fits that quantify the transition to
incompressible turbulence ($\beta=\Delta=2/3$) with decreasing Mach
number likely apply~\citep{padoan-et-al:04}. If $\Delta \ne 2/3$, two
parameter fits for $\beta$ and $\Delta$ would also be needed in the
mildly compressible regime.  However, the analysis in
Sect.~\ref{sec:accu_needed} suggests that these kind of fits are
likely beyond reach of today's computer resources.

The latter of the above cases, where the $Z_{\mathrm{p}}$ curves do not
intersect to within their uncertainty, would imply that the model
by~\citet{she-leveque:94} does not carry over to 3D periodic box
simulations of driven, isothermal, highly compressible turbulence. A
simple reason here could be that theoretical results are based on the
assumption of an infinite Reynolds number, a criterion clearly
violated by numerical simulations. More importantly,
\citet{she-waymire:95} already pointed out that there is no reason why
only one dimension should be associated with the most dissipative
structure. They argued that in such a large portion of space as is
typically analyzed, a variety of most dissipative structures may
co-exist with different co-dimensions.  \citet{hopkins:13} suggests a
slightly different model based on work by~\citet{castaing:96}, which
is more compatible with not strictly log-normal density PDFs as
observed in isothermal supersonic turbulence. Finally, yet other
models exist,~\cite[e.g.  via multifractals][]{macek:11,
  zybin-sirota:13}, as well as other perspectives on the fractal
character of a turbulent medium~\citep[see e.g.][]{kritsuk-et-al:07}.

Lastly, we briefly come back to the one parameter fits that are often
used in the literature. Fixing the value of $\Delta$ by hand greatly
reduces the impact of uncertainties in the $Z_{p}$ on the accuracy of
the estimated best-fit co-dimension $C_{\mathrm{f}}$.
\citet{folini-et-al:14} found 5\% uncertain $Z_{p}$, $p=1,.\,.\,.\,,5$, to
translate into roughly 10\% uncertainty of the $C_{\mathrm{f}}$ for
fixed $\Delta=2/3$.  Fig.~\ref{fig:dmap} offers a qualitative
understanding of this reduced "error propagation", which suggests some
sensitivity to the specific location in the $\beta$-$\Delta$-plane,
and indicates that fixing $C$ or $\beta$ instead of $\Delta$ has a
similar effect.  Fixing $\beta$ seems questionable at first sight
since there is, to our knowledge, little theoretical understanding of
what numerical value $\beta$ might have~\citep[see
e.g.][]{dubrulle:94}. On the other hand, \citet{she-et-al:01}
presented a theoretical framework that allows for an independent
determination of only $\beta$ from the relative scaling exponents
$Z_{p}$ (see also~\citet{hily-blant-et-al:08}, their Appendix A3). One
could thus imagine breaking the two parameter fit for $\beta$ and
$\Delta$ into a two step procedure: first, fix the value of $\beta$,
then use this value and do a one parameter fit for $\Delta$.
Hopefully this type of a two step approach is more robust against
uncertainties in the $Z_{p}$, but this question is beyond the scope of
the current paper and we are unaware of corresponding attempts in the
literature.

\section{Summary and conclusions}
\label{sec:conc}
This study was motivated by the overarching question of whether or not
the random cascade model~\citep{she-leveque:94, dubrulle:94,
  she-waymire:95, boldyrev:02} applies to simulation data of highly
compressible isothermal turbulence and, if so, with what parameter
values for $\beta$ and $\Delta$. If applicable, the model offers a
theoretical link between observable properties of the turbulence,
namely ratios $Z_{p}$ of scaling exponents of the structure functions,
and non-observable turbulence characteristics, for example the
dimension $D$ of the most dissipative structures. To date,
applicability of the model is assumed in much of the literature with
$\Delta=2/3$, a value just marginally compatible with simulation-based
best estimates~\citep{pan-et-al:09}: $\Delta=0.71$ with an uncertainty
range $\Delta \in (0.67, 0.78)$.

We examine how uncertainties in the $Z_{p}$ translate into
uncertainties of best-fit $\beta$-$\Delta$-pairs and discuss what
best-fits, consequently, seem achievable with today's computer
resources. A Monte Carlo approach is used to mimic actual simulation
data. The results can be summarized in six main points.

\begin{itemize}
\item{Simultaneous fitting of $\beta$ and $\Delta$ to sets of
    substantially (5\%) perturbed (uncertain) $Z_{p}$ yields a "double
    peaked ridge" of best-fit values in the $\beta$-$\Delta$ plane.
    None of the two peaks is co-located with the initial $(\beta,
    \Delta)$ pair.}

\item{The highest and lowest order $p$ are particularly relevant for
    simultaneous fitting of $\beta$ and $\Delta$. A somewhat optimal
    choice is $p=1,.\,.\,.\,,5$. Yet higher order structure functions add
    comparatively little to the quality of the fit, while they tend to
    be afflicted with larger uncertainties in real applications.}

\item{A simultaneous, 2\% (1\%) accurate fit of $\beta$ and $\Delta$
    should be possible if the $Z_{p}$, $p=1,.\,.\,.\,,5$, are~0.1\% (0.05\%)
    accurate and if the (yet to be determined) value of $\beta$ is
    about 1/3 or less.}

\item{Applicability of the model thus may be best tested in the highly
    compressible regime, where $\beta \approx 1/3$ is expected, and
    not in the mildly compressible regime where $\beta$ ultimately
    must approach its incompressible value of 2/3.}

\item{We argue that today's computer resources likely allow to reach
    this accuracy. Existing simulations of
    $4096^{3}$~\citep{federrath:13, beresnyak:14} probably allow for
    at least 2\%, possibly 1\% accurate estimates of $\beta$ and
    $\Delta$.}

\item{Should the ambiguity in the determination of $\beta$ and
    $\Delta$ persist despite such highly accurate $Z_{p}$, this may
    indicate that the notion of~\citet{she-waymire:95} ($\beta$ and
    $\Delta$ take a continuum of values) or~\citet{hopkins:13} (a
    different model for the statistics of the inertial range) is
    correct or that yet a different turbulence model is needed in this
    regime.}
\end{itemize}

While the authors lack the computational resources to really test the
estimates presented here, this study may encourage other groups to
analyze their data in the light of this study.
\begin{acknowledgements}
  RW and DF acknowledge support from the French National Program for
  High Energies PNHE.  We acknowledge support from the P\^{o}le
  Scientifique de Mod\'{e}lisation Num\'{e}rique (PSMN), from the
  Grand Equipement National de Calcul Intensif (GENCI), project number
  x2014046960, and from the European Research Council through grant ERC-AdG
  No. 320478-TOFU.

\end{acknowledgements}
\bibliographystyle{apj}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
\label{intro}
After the Higgs was discovered in 2012 \cite{higgs1,higgs2}, understanding its properties, and looking for any possible deviations of its properties from the SM prediction, becomes a very important task of LHC. This includes measuring the its spin (scalar or tensor), its mass which has important implications for the vacuum stability and cosmology, its flavor couplings to look for Beyond Standard Model (BSM) signatures of flavor changing effect, its exotic decay modes such as $H\to\text{invisible}$, and its CP property which is related to the matter-antimatter imbalance in the universe. If the Higgs is in a pure CP eigen state, we will find if it is CP even or odd. On the other hand, if it is in a CP mixture, we will measure the CP mixing angle. This is a more exciting scenario, as it gives rise to an extra CP violation source, and the current known source (a single complex phase in CKM) is too small to explain the matter-antimatter imbalance.

\section{CP test in the bosonic decays of the Higgs}
\label{boson_decay}

The CP symmetry of the Higgs coupling to vector bosons are parametrized with an effective Lagrangian by the ATLAS collaboration as \cite{EFT1,EFT1b}
\begin{eqnarray}
\begin{array}{ll}
\mathcal{L}_0^V = & \left\{ c_\alpha \kappa_{SM}  \left[ \frac{1}{2} g_{HZZ} Z_\mu Z^\mu +  g_{HWW}W_\mu^+ W^{-\mu} \right] -\frac{1}{4} \frac{1}{\Lambda}  \left[ c_\alpha \kappa_{HZZ} Z_{\mu\nu} Z^{\mu\nu} + s_\alpha \kappa_{AZZ} Z_{\mu\nu} \tilde{Z}^{\mu\nu}  \right] \right.  \\
& \left. -\frac{1}{2} \frac{1}{\Lambda} \left[ c_\alpha \kappa_{HWW}  W_{\mu\nu}^+ W^{-\mu\nu} + s_\alpha \kappa_{AWW}  W^+_{\mu\nu} W^{-\mu\nu} \right] \right\} X_0 ,
\end{array}
\label{eq:eq1}
\end{eqnarray}
where $\alpha$ is the CP mixing angle, $c_\alpha = \cos{\alpha}$, $s_\alpha = \sin{\alpha}$, $\Lambda$ is scale cut-off for dim-6 operators, $\tilde{V}_{\mu\nu}(V=W/Z)$ is the dual tensor of $V_{\mu\nu}$, and $X_0$ is the neutral resonance such as the Higgs. The terms with $1/\Lambda$ coefficients are BSM dim-6 operators, with the $c_\alpha$ ($s_\alpha$) terms representing the BSM CP even (odd) contributions. Based on Eq.~\ref{eq:eq1}, the pure CP states of Higgs can be parametrized as in Tab.~\ref{tab:tab1}. CMS collaboration used a similar expression, and the part with $Z$ boson terms only are \cite{EFT2}
\begin{equation}
\mathcal{L}_{HVV} \sim a_1 \frac{m_Z^2}{2} HZ^\mu Z_\mu - \frac{\kappa_1}{\Lambda_1^2} m_Z^2 H Z_\mu \square Z^\mu - \frac{1}{2} a_2 H Z^{\mu\nu}Z_{\mu\nu} - \frac{1}{2} a_3 HZ^{\mu\nu}\tilde{Z}_{\mu\nu},
\label{eq:eq2}
\end{equation}
\begin{table}[h]
\caption{ The definitions of $HVV$ coupling coefficients for pure Higgs CP states (ATLAS) \cite{EFT1}. }
\label{tab:tab1}
\centering
\begin{tabular}{cccccc}
\hline \hline
$J_P$ & Model & \multicolumn{4}{c}{Choice of tensor couplings} \\
& & $\kappa_{SM}$ &  $\kappa_{HVV}$ &  $\kappa_{AVV}$ & $\alpha$ \\ \hline
$0^+$ & Standard Model Higgs boson & 1 & 0 & 0 & 0 \\
$0_h^+$ & BSM spin-0 CP-even & 0 & 1 & 0 & 0 \\
$0^-$ & BSM spin-0 CP-odd & 0 & 0 & 1 & $\pi/2$ \\ \hline \hline
\end{tabular}
\end{table}
Both ATLAS and CMS have tested the CP states of the Higgs assuming it is in a pure CP state. Both experiments used a Matrix Element (ME) method to test the $0^+$ state against $0^-$.  ATLAS combined three channels in the test: $H\to WW\to e\nu \mu\nu$ 0-jet category, $H\to\gamma\gamma$ and $H\to ZZ\to 4l$. The results are as shown in Fig.~\ref{fig:fig1}. CMS combined $H\to WW, ZZ, Z\gamma, \gamma\gamma$ channels and for the pure CP test \cite{EFT2}.  The CP odd $0^-$ model is excluded at CL better than $99.9\%$ by both experiments.
\begin{figure}[h]
\centering
\includegraphics[width=6.cm,clip]{fig_01.png}
\caption{ The test results of Higgs with different spin and CP states at ATLAS \cite{EFT1}. The measurements are indicated by the black solid bars. }
\label{fig:fig1}
\end{figure}
When the Higgs is in a CP-mixed state, both experiments have made constraints on the coefficients of different couplings in Eq.~\ref{eq:eq1} and \ref{eq:eq2}, as shown in Tab.~\ref{tab:tab2} and Fig.~\ref{fig:fig2}. One can conclude that the non-SM tensor couplings are consistent with zero for both experiments.
\begin{table}[h]
\caption{ The best-fit and excluded ranges for the coefficients of the BSM CP-even and CP-odd terms in Eq.~\ref{eq:eq1} at ATLAS, with $H\to ZZ^*$ and $H\to WW^*$ combined \cite{EFT1}.}
\label{tab:tab2}
\centering
\begin{tabular}{ccccc}
\hline \hline
Coupling ratio & \multicolumn{2}{c}{Best fit value} & \multicolumn{2}{c}{$95\%$ CL Exclusion Regions} \\
Combined & Expected & Observed & Expected & Observed \\ \hline
$\tilde{\kappa}_{HVV}/\kappa_{SM}$ & 0.0 & -0.48 & $(-\infty,-0.55]\cup[4.80,\infty)$ & $(-\infty,-0.73]\cup[0.63,\infty)$ \\
$(\tilde{\kappa}_{AVV}/\kappa_{SM})\cdot\tan\alpha$ & 0.0 & -0.68 & $(-\infty,-2.33]\cup[2.30,\infty)$ & $(-\infty,-2.18]\cup[0.83,\infty)$ \\
\hline \hline
\end{tabular}
\end{table}
\begin{figure}[h]
\centering
\includegraphics[width=12.cm,clip]{fig_03.png}
\caption{ The bounds on the coefficients of the different BSM terms for $H\to ZZ^*$, $H\to WW^*$, $H\to Z\gamma$ and $H\to \gamma\gamma$ decays at CMS \cite{EFT2}. }
\label{fig:fig2}
\end{figure}

Another important way to search for Higgs CP violation is in the Vector-Boson-Fusion (VBF) production of the Higgs. One of the most promising channel for this search is $H\to\tau\tau$ decay. The combined Run-1 sensitivity of ATLAS and CMS has already exceeded $5\sigma$ \cite{Htt}. This channel is special in that it is sensitive to CP in both $HVV$ and $Hf\bar{f}$ couplings. In the 2HDM model, the rate of $H\to\tau\tau$ is often enhanced, and mixed CP for Higgs can be accommodated.

Omitting the BSM CP-even terms, the $HVV$ coupling can be written as \cite{OO}
\begin{equation}
\mathcal{L}_{eff} = \mathcal{L}_{SM} + \tilde{g}_{HAA}H\tilde{A}_{\mu\nu}A^{\mu\nu} + \tilde{g}_{HAZ}H\tilde{A}_{\mu\nu}Z^{\mu\nu} + \tilde{g}_{HZZ}H\tilde{Z}_{\mu\nu}Z^{\mu\nu} + \tilde{g}_{HWW}H\tilde{W}^+_{\mu\nu}W^{-\mu\nu}.
\label{eq:eq3}
\end{equation}
Since the vector bosons in the process are not directly observable, it is easier to treat the CP-odd terms as one by the assumption in Eq.~\ref{eq:eq4} about the relations of the coefficients:
\begin{equation}
\tilde{g}_{HAA} =  \tilde{g}_{HZZ} = \frac{1}{2} \tilde{g}_{HWW} = \frac{g}{2m_W} \tilde{d}, ~~\text{and}~~\tilde{g}_{HAZ}=0.
\label{eq:eq4}
\end{equation}
Traditionally, the signed $\Delta\phi$ between the two tagging in the VBF process are used as the CP sensitive variable \cite{VBF}. ATLAS used a different variable, the Optimal Observable (OO), which is expected to perform better. It is based on the ME of
\begin{equation}
\mathcal{M} = \mathcal{M}_{SM} + \tilde{d} \cdot \mathcal{M} _{CP-odd},  ~~\text{and}~~
\left|\mathcal{M}\right|^2 = \left|\mathcal{M}_{SM}\right|^2 + \tilde{d} \cdot 2 Re\left( \mathcal{M}_{SM}^*\mathcal{M}_{CP-odd} \right) +  \tilde{d}^2 \left| \mathcal{M}_{CP-odd} \right|^2,
\label{eq:eq5}
\end{equation}
and the OO is defined as
\begin{equation}
OO = \frac{2 Re \left(\mathcal{M}_{SM}^* \mathcal{M}_{CP-odd}  \right)}{\left| \mathcal{M}_{SM} \right|^2}.
\label{eq:eq5b}
\end{equation}

With all 4-momenta of the final state particles (the Higgs and two tagging jets) measured, the LO ME of SM and CP-odd can be calculated from HAWK \cite{hawk}, and then the OO can be calculated per event. As Fig.~\ref{fig:fig3} shows, of there is no CP violation, the mean of the OO distribution should be zero. For positive (negative) CP-odd component (determined by $\tilde{d}$), its mean will be shifted to positive (negative) values. This method can be applied to other decays such as $H\to\gamma\gamma$.
\begin{figure}[h]
\centering
\includegraphics[width=6.cm,clip]{fig_04.png}
\caption{ The distributions of OO for the SM ($\tilde{d}=0$) and two other values of $\tilde{d}$ in the VBF $H\to\tau\tau$ signal at ATLAS \cite{OO}. }
\label{fig:fig3}
\end{figure}
To obtain a pure signal sample, a cut is first made on the Multi-Variate-Analysis (MVA) output score. A signal to background ratio of about 0.3 is achieved by this cut. Next, a likelihood fit to the OO distribution is performed to find the best value for $\tilde{d}$. Figure~\ref{fig:fig4} shows the OO distribution in the $H\to\tau\tau\to ll+4\nu$ subchannel, and the increase of the negative-log-likelihood $\Delta$NLL with respect to its minimum in the likelihood scan.  Each point in the plot indicates a $\Delta$NLL calculated with a particular hypothesis template and the data. The $68\%$ CL interval is found by the intersection points of the $\Delta$NLL curve and the horizontal line at $\Delta\text{NLL}=0.5$. The result is $-0.11<\tilde{d}<0.05$. While the result shows that the CP-odd contribution is consistent with zero, it is about 10 times tighter than the one shown in Tab.~\ref{tab:tab2} for the $HVV$ coupling.
\begin{figure}[h]
\centering
\includegraphics[width=5.7cm,clip]{fig_05.png}\hspace{5mm}
\includegraphics[width=6cm,clip]{fig_06.png}
\caption{ The distributions of the OO in the $H\to\tau\tau\to ll+4\nu$ subchannel (left), and the $\Delta$NLL for different $\tilde{d}$ hypotheses (right) at ATLAS \cite{OO}.}
\label{fig:fig4}
\end{figure}
In the HL-LHC phase, the Missing Transverse Energy and the di-tau mass resolution will degrade as the pileup increases. The signal-background discrimination increases with the tracking coverage, as the VBF forward jets are better separated from the pileup jets with the association between the jets and the primary vertex, shown in Fig.~\ref{fig:fig5}. Assuming zero theory uncertainty, the uncertainty on the signal strength $\mu$ at 8-18\% (2-5\%) can be achieved with the ATLAS $\tau_l \tau_h$ alone \cite{HL1} (CMS $\tau\tau$ inclusive \cite{HL2}). The projected uncertainty of $\mu$ for different tracking extensions at ATLAS is listed in Tab.~\ref{tab:tab3}. Thus, we would expect a corresponding increase in the precision for the Higgs CP measurement in the VBF channel due to both increased data and extended tracking coverage.
\begin{figure}[h]
\centering
\includegraphics[width=5.0cm,clip]{fig_07.png}
\caption{ The jet-PV association efficiency as a function of the VBF quark $\eta$ in the $H\to\tau\tau$ channel at CMS. The improvement is evident with the tracking volume extension \cite{HL2}. }
\label{fig:fig5}
\end{figure}
\begin{table}[h]
\caption{ The projected uncertainty on the signal strength for different scenarios of tracking coverage and forward pileup jet rejection at HL-LHC for ATLAS \cite{HL1}.}
\label{tab:tab3}
\centering
\begin{tabular}{c|ccc}
\hline \hline
forward pile-up jet rejection & $50\%$ & $75\%$ & $90\%$ \\ \hline \hline
forward tracker coverage & \multicolumn{3}{c}{$\Delta\mu$} \\ \hline
Run-I tracking volume & \multicolumn{3}{c}{0.24} \\
$\eta<3.0$ & \multicolumn{1}{c| }{0.18} & \multicolumn{1}{c| }{0.15} & \multicolumn{1}{c}{0.14} \\
$\eta<3.5$ & \multicolumn{1}{c| }{0.18} & \multicolumn{1}{c| }{0.13} & \multicolumn{1}{c}{0.11} \\
$\eta<4.0$ & \multicolumn{1}{c| }{0.16} & \multicolumn{1}{c| }{0.12} & \multicolumn{1}{c}{0.08} \\
\hline \hline
\end{tabular}
\end{table}

\section{CP test in the fermionic decays of the Higgs}
\label{fermion_decay}

The CP-odd Yukawa coupling can enter the Lagrangian as a dimension-4 operator as in Eq.~\ref{eq:eq6}, thus the $Hf\bar{f}$, especially the  $H\tau\tau$ coupling, is a sensitive probe of CP at tree-level rather than the loop level as with the dimension-6 operators in the $HVV$ coupling.
\begin{equation}
\mathcal{L} = -g_\tau \left( \cos\phi\bar{\tau}\tau + i\sin\phi\bar{\tau}\gamma_5\tau \right)h,
\label{eq:eq6}
\end{equation}
where $\phi$ is the mixing angle between the CP even and odd terms. The CP of $H\tau\tau$ coupling can be distinguished by the transverse tau spin correlations, as the decay width is proportional to
\begin{equation}
\Gamma(H,A\to\tau^-\tau^+) \sim 1 - s_z^{\tau-}s_z^{\tau+} \pm s_T^{\tau-}s_T^{\tau+} ,
\label{eq:eq7}
\end{equation}
where $s_{z,T}$ are the tau spin components in the longitudinal and transverse directions with respect to the tau momentum. As a result, the angle between the planes spanned by the tau and its charged track is sensitive to the CP. For example, in the $\tau\to\pi\nu$ decay, one can look at the angle between tau decay planes to extract the CP mixing angle $\phi$:
\begin{equation}
\frac{d\Gamma\left(H\to\tau\tau\to\pi^+\pi^-+2\nu \right)}{d\phi_{CP}} \sim 1-\frac{\pi^2}{16}\cos\left(\phi_{CP}-2\phi \right),
\label{eq:eq8}
\end{equation}
where $\phi_{CP}$ is the angle between the tau decay planes in the ditau rest frame.

Using the $H\to\tau\tau$ decay to measure the CP is experimentally challenging because the neutrinos are not reconstructed. There are two main methods to extract the CP \cite{CP1,CP2,CP3,CP4,CP5}. One is by using the impact parameters to approximately reconstruct the tau decay plane from its leading track. This is best done for the $\tau\to\pi\nu$ decay, and the analyzing power is compromised for the other tau decays. The other is by using the $\tau\to\rho\nu\to\pi^\pm\pi^0\nu$ decay. The tau decay plane can be approximately reconstructed by the track and the neutral pion. However, the relative energy of $\pi^\pm$ and $\pi^0$ need to be classified in order to maximize the analyzing power. In order to use the methods, the tau decay modes (or its substructure) need to be well differenciated. In ATLAS \cite{sub1}, after the initial tau reconstruction and identification, the hadronic tau candidate is further analyzed with a particle flow (PF) algorithm to better reconstruct and differenciate the neutral pions. The energy in the calorimeter deposited by $\pi^\pm$ is then subtracted, and the remaining energy cluster is reconstructed as $\pi^0$ and identified with a Boosted Decision Tree (BDT) method. The final decay mode classification is done with another BDT. The reconstruction and identification efficiencies for a single tau at ATLAS is listed in Tab.~\ref{tab:tab4}. The efficiency for separating the different decay modes and the amount of crosstalk is visualized in Fig.~\ref{fig:fig7}. In general, there is a non-negligible fraction of 2 (1) $\pi^0$ reconstructed as 1 (0) $\pi^0$. In Fig.~\ref{fig:fig8}, the tau energy resolution with the PF is compared with the old baseline method, and the reconstructed tau mass and relative abundance of different modes are compared with data. In conclusion, with the substructure, a factor of 2 improvement of tau energy with respect to the calo-based method at the low $p_T$ (the resolution of the neutral $\pi^0$ is about 16\%), and a factor of 5 improvement in the angular resolution can be achieved. For the neutral $\pi^0$, a precision of 0.006 (0.012) can be reached for it $\eta$ ($\phi$) measurement.
\begin{table}[h]
\caption{ The Branching Fractions, acceptance, reconstruction and identification efficiencies for a single tau in different decay modes at ATLAS \cite{sub1}. }
\label{tab:tab4}
\centering
\begin{tabular}{lccc}
\hline
Decay mode & $\mathcal{B}$ [\%] & $\mathcal{A}\cdot\epsilon_{reco}$ [\%] & $\epsilon_{ID}$ [\%] \\ \hline
$h^\pm$ & 11.5 & 32 & 75 \\
$h^\pm\pi^0$ & 30.0 & 33 & 55 \\
$h^\pm\geq 2\pi^0$ & 10.6 & 43 & 40 \\
$3h^\pm$ & 9.5 & 38 & 70 \\
$3h^\pm\geq 1\pi^0$ & 5.1 & 38 & 46 \\ \hline
\end{tabular}
\end{table}
\begin{figure}[h]
\centering
\includegraphics[width=4.7cm,clip]{fig_10.png}\hspace{5mm}
\includegraphics[width=4.7cm,clip]{fig_11.png}
\caption{ The efficiencies for the generated tau decay modes to be reconstructed as different measured modes (left), and the fractions of different reconstructed modes coming from different generated modes indicating the purity of the reconstruction (right), at ATLAS \cite{sub1}. }
\label{fig:fig7}
\end{figure}
\begin{figure}[h]
\centering
\includegraphics[width=4cm,clip]{fig_12.png}\hspace{3mm}
\includegraphics[width=4cm,clip]{fig_13.png}\hspace{3mm}
\includegraphics[width=4cm,clip]{fig_14.png}
\caption{ The resolution of the reconstructed tau energy (left), the mass (middle) and modes (right) of the reconstructed taus with substructure at ATLAS \cite{sub1}. }
\label{fig:fig8}
\end{figure}
CMS also uses the PF constituents (charged and neutral particles) to reconstruct taus \cite{sub2,sub3}. Discriminants are then applied to reject jets and leptons (with a MVA based tau identification). Multiple $\tau_h$ hypotheses for each jet are constructed by a combinatorial approach. The one passing all cuts and with the highest $p_T$ is selected as the decay mode for the tau candidate. The mode reconstruction efficiency, and tau mass with Run-1 and Run-2 data are shown in Fig.~\ref{fig:fig9}. Compared with ATLAS, CMS currently has less decay modes classified, but the crosstalk among different modes is also smaller than ATLAS. Good agreement with Run-1 and Run-2 data is achieved, although the $W$+jets and QCD background are visibly higher in the latter.
\begin{figure}[h]
\centering
\includegraphics[width=4cm,clip]{fig_15.png}\hspace{3mm}
\includegraphics[width=4.0cm,clip]{fig_16.png}\hspace{3mm}
\includegraphics[width=3.8cm,clip]{fig_17.png}
\caption{ The efficiencies for the generated tau decay modes to be reconstructed as different measured modes (left), the tau mass distributions with Run-1 (middle) and Run-2 (right) data, at CMS \cite{sub2,sub3}. }
\label{fig:fig9}
\end{figure}
Both ATLAS and CMS used the embedding technique to estimate the $Z\to\tau\tau$ background with reduced systematics. This is done by replacing muons from the $Z\to\mu\mu$ data by simulated taus and decaying the taus, and then embedding the simulated tau decay products into the original real data event. No spin correlation information is conserved by this procedure. The Tauola and TauSpinner    packages \cite{tsp} could be used to ensure the correct $Z$ boson polarization and to restore the di-tau spin correlation.

\section{Conclusion}
\label{conclusion}
Testing the CP nature of Higgs is one of the important tasks after its discovery. With the diboson channels ($H\to\gamma\gamma, ZZ, WW$), and assuming pure CP state for the Higgs, the CP even state is favored by data. The tensor structure of the $HVV$ coupling can also be probed for mixed CP scenarios. No significant CP mixing effect is observed and limits are set on the CP-odd terms in the effective Lagrangian. The $H\to\tau\tau$ channel is ideal for probing CP through both the $HVV$ coupling through VBF, and the $H\tau\tau$ Yukawa coupling. Combined sensitivity of ATLAS and CMS already exceeded $5\sigma$ with the Run-1 data. The former has results with VBF $H\to\tau\tau$ based on the Run-1 data. The CP-mixing parameter in the VBF production is currently consistent with zero. The CP mixing in $H\to\tau\tau$ decay can be large in theory, but experimentally very challenging. Both ATLAS and CMS have refined the tau reconstruction with substructure information, which is essential for the CP study in the tau decays. Looking forward, the CP test precision will improve with the HL-LHC data, in which just a few percents uncertainty on the $H\to\tau\tau$ signal strength is expected.


text:
\section{Introduction}
The fingerprinting method is a well-known technique generally used for determining the identity of a person. The basic concept is illustrated in Fig.~\ref{fig:fingerprint_principle}. The overall process can be decomposed into three different steps: \textit{(i)} a fingerprint recording of an ensemble of subjects, \textit{(ii)} the creation of a database (also called dictionary) where fingerprint images are associated with person identities, and \textit{(iii)} a recognition process where a numerical search algorithm finds the closest database element to the fingerprint of an unknown subject. Assuming that fingerprints are different for each person, a mapping between fingerprints and persons can be defined, making possible the identification protocol.
\begin{figure}
	\includegraphics[width=1.0\linewidth]{fig1}
	\caption{(Color online) The fingerprinting principle: A physical system is identified by using a pattern recognition between a fingerprint measurement and the elements of a database.}\label{fig:fingerprint_principle}
\end{figure}
This idea can be generalized to any system which has unique properties that can be revealed by a measurement process. This approach can be applied in a static setting, but also in a dynamical one where the system is subjected to an external control field. In this latter case, each element of the database corresponds to the time evolution of some observables under the action of the field, thus increasing the complexity of the fingerprints and the precision of the estimation. This idea has been recently adapted to Magnetic Resonance Imaging (MRI) for the identification of tissue parameters \cite{naturemri}. This initial investigation led to an impressive number of studies in this domain (see e.g. \cite{naturemri2,jacobpaper} to cite a few). A crucial issue in this fingerprinting process is the design of the excitation field. A simple approach using a time-dependent random field was proposed in \cite{naturemri} in order to limit the correlations between the different fingerprints. However, this approach does not incorporate information about the system dynamics and the recognition process and is therefore not expected to reach its precision limit. To overcome this fundamental difficulty, we propose in this paper to combine a standard fingerprinting process with recently developed optimal control techniques in quantum control \cite{rabitzreview,quaintreview,grape}. This Optimal Fingerprinting Process (OFP) allows us to maximize the efficiency of the identification and to minimize the errors made in the parameter estimation. As an illustrative example, this method will be used to identify the relaxation parameters of a spin- 1/2 particle. In this case, the estimation is made from a series of free-induction decay signals induced by impulsive excitations of different intensities. Note also that the measurements we consider are classical as they result from a continuous measurement of a large number of quantum systems. Nevertheless, the same type of processes \cite{haidong_metrology,haidong_metrology2} could be carried out with quantum measurements in quantum metrology \cite{quantummetrology1,quantummetrology2}. Closely related but different concepts using external fields to estimate the parameters of quantum systems have also been developed in the past few years \cite{geremia1,geremia2,turinici,schirmer1,schirmer2,schirmer3,li1,li2}. The combination of optimization and fingerprinting techniques has been explored in recent works in MRI~\cite{jacobpaper,maidens}.

The paper is organized as follows. Section~\ref{sec2} describes the theoretical framework of the method. The technique is applied in Sec.~\ref{sec3} to estimate the relaxation parameters of a spin system. The stability of the estimation in presence of noise is discussed and analyzed numerically. The efficiency of the optimal fingerprinting process is demonstrated experimentally in Sec.~\ref{sec4} on a spin 1/2 particle by using techniques of Nuclear Magnetic Resonance. Conclusion and prospective views are given in Sec.~\ref{conc}. A mathematical description of the method and numerical results are reported respectively in the Appendices \ref{appa}, \ref{appb} and \ref{appc}.
\section{Theoretical framework}\label{sec2}
We start the analysis with a general presentation of the method on an abstract system. We refer the reader to the Appendix~\ref{appa} for a detailed mathematical description of this technique. The state of the system is given at a time $t$ by $\Psi(t)\in\mathcal{H}$ ($\mathcal{H}$ is generally a Hilbert space), whose dynamics are governed by the following differential equation:
\begin{equation}\label{eq1}
\dot{\Psi}=\hat H (\vec{S},u(t))\Psi,
\end{equation}
where $\hat H$ is a linear operator (the Hamiltonian for quantum systems), $\vec{S}\in\mathbb{R}^p$ defines the $p$ unknown parameters to estimate and $u(t)$ is the control field.
For a generic control field, the resulting time evolution $\Psi(t)$ will be different for each system characterized by different $\vec{S}$. The experimental system returns a specific response $g(\vec{S}_{0},t)$. This response defines the fingerprint of the system and is assumed to be unique for a given vector $\vec{S}_{0}$. In this measurement process, note that experimental noise has to be accounted for. This point will be investigated in the section about the case study in NMR.

The database is built from numerical simulations of the time evolution of $N$ systems characterized by specific values of the $\vec{S}$ parameters, denoted $\{\vec{S}_n\}_{n=1,\cdots, N}$. The database is defined as a set $\{f_n(t)\}_{n=1,\cdots ,N}$ of $N$ real functions associated with each $\vec{S}_n$.
The vector $\vec{S}_{0}$ is determined from the best match between $g(t)$ and one of the elements $f_n(t)$ of the dictionary. This leads to the estimation $\vec{S}_{0}\simeq \vec{S}_k$ for the experimental system, where $k\in \{1,\cdots , N\}$. The match is performed by minimizing the functional $D$ (called the \textit{recognition map}) over the elements of the dictionary. $D$ is defined as follows:
\begin{equation}\label{cost1}
D[f_n,g]=\left\Vert \frac{f_n}{||f_n||} - \frac{g}{||g||}\right\Vert^2.
\end{equation}
In Eq.~\eqref{cost1}, the two vectors are divided by their norm to eliminate a possible scaling factor between the experimental data and the theoretical model.

At this point, the method can be performed for any control field that distinguishes the elements of the dictionary. OFP is defined by introducing a
figure of merit which is aimed at maximizing the distance between the elements of the dictionary, and thus improving the recognition process and the precision of the method. This functional $C$ can be expressed as:
\begin{equation}
C_N = \frac{1}{2N^2}\sum_{m,n} \mu_{m n} D\left[f_m,f_n \right],
\end{equation}
where the $\mu_{m n}$ are some weight factors. The parameter $C_N$ is the normalized average distance between the $N$ elements of the dictionary. The normalization factor ensures that $C_N\leq 1$ when all weights $\mu_{m n}=1$ (see Appendix~\ref{appa}). Maximizing this quantity allows the minimization of the overlap of $g$ with the other functions of the dictionary and thus the error made in the estimation procedure. As described in Appendix~\ref{appb}, a generalized version of the optimal control GRAPE algorithm \cite{grape,trackinggrape} can be used to numerically generate the control field. We stress that this extension is not trivial since the maximization of the distance between the different systems is not made at one point but for the whole time evolution. OFP is schematically described in Fig.~\ref{fig:map}.
\begin{figure}[h]
	\includegraphics[width=1.0\linewidth]{fig2}
	\caption{(Color online) The OFP is composed of two different loops. The first loop (yellow or light gray arrows) is the standard fingerprinting process. A control field $u(t)$ is designed at the starting point of the loop. This field is applied to a physical system which returns a specific response $g(t)$. On the other side, the response is computed numerically for an ensemble of physical systems with different values of the parameters. These simulations define a \textit{dictionary} of functions $f_n(t)$. The recognition process allows us to find the best match between elements of the dictionary $f_n$ and the result of the measurement $g$ (see the text for details). The second loop (orange or dark gray arrows) describes the dictionary optimization. The optimization is performed for an ensemble of $N$ systems with different values of the parameters. An optimal control algorithm is used to maximize numerically the figure of merit $C_N$.\label{fig:map}}
\end{figure}

\section{A case study in Nuclear Magnetic Resonance}\label{sec3}
As an illustrative example, we investigate the estimation of the relaxation parameters of a spin system by OFP in NMR~\cite{ernstbook,levittbook}. We consider an inhomogeneous ensemble of spin- 1/2 particles with different resonance offsets $\omega$ and radio-frequency inhomogeneities $\alpha$ whose dynamics are ruled by the Bloch equations. In the rotating frame, the equation of motion for each isochromat is given by~\cite{levittbook,kozbar}:
\begin{equation}\label{eq:bloch}
\frac{d}{dt} \vec M
    =
    \left(
    \begin{array}{cccc}
    1 & 0&0&0\\\
    0             & -\frac{1}{T_2} & -\omega     & \alpha\omega_y(t) \\
    0             & \omega     & -\frac{1}{T_2} &  -\alpha\omega_x(t) \\
    \frac{1}{T_1} & -\alpha\omega_y(t)  & \alpha\omega_x(t) & -\frac{1}{T_1}
    \end{array} \right)
    \vec M
\end{equation}
where $\vec M = (1,M_x,M_y,M_z)^t$ is the extended Bloch vector (the radius of the Bloch ball is normalized to 1) and $M_{x,y,z}^{(\omega)}$ its coordinates along the $x,y,z$- directions. In this case, note that $\mathcal{H} = \mathbb{R}^4$. The relaxation times $T_1$ and $T_2$ are assumed to be the same for all the isochromats of the sample. The control amplitudes are given by $\omega_x(t)$ and $\omega_y(t)$. $\omega$ is the resonance offset and the parameter $\alpha$ describes the experimental scaling of the radio-frequency field applied to the sample~\cite{kozbar}. In this example, OFP is used to estimate the parameters $T_1$ and $T_2$, and the measured signal results from the average magnetization of the spins with different values of $\omega$ and $\alpha$. The functions of the dictionary are given by
\begin{equation}\label{eqdicspin}
f_n(t) = \big(\bar M_x^{(n)}(t),\bar M_y^{(n)}(t) \big),
\end{equation}
where $\bar M_\mu$, $\mu=x,y$, is the average of $M_\mu$ over the sample.
The averaging procedure is defined by using a probability distribution in $\omega$ and $\alpha$ which can either be known before the optimization of the database or adjusted during the recognition process (see below for an example).

The control field is a sequence of short pulses (with respect to $T_1$ and $T_2$), modeled by Dirac distributions:
\begin{align}
\omega_\mu(t) = \sum_{k=1}^{N_p} \omega_{\mu,k} \delta (t-kT), && \mu = x,y,
\end{align}
where $\omega_{\mu,k}$ is the amplitude of the $k$-th pulse, $N_p$ the number of pulses and $T$ the time between each pulse.
This approximation leads to a straightforward time discretization of the dynamics of the system. The measured signal corresponds to the average magnetization just after each $\delta$- pulse, $\big(\bar M_x^{(n)}(t=kT),\bar M_y^{(n)}(kT) \big)$ with $k=1,\cdots ,N_p$.

\begin{figure}[h]
	\includegraphics[width=\linewidth]{fig3}
	\caption{(Color online) The dictionary is composed of 4 elements regularly distributed in the interval $T_1\in[0.1,0.5]$~s. Giving an optimal (d) and a random (h) control fields (the black and the red (dark gray) lines represent respectively $\omega_x$ and $\omega_y$), we can compute the associate square modulus of dictionary functions (c) and (g). The dictionary functions are dimensionless. The efficiency can be checked with the  recognition maps $(T_1(m),T_1(n))\mapsto D[f_m, f_n]$ in (b) and (f). The panels (a) and (e) show the distance between the elements of the dictionary and the system to identify ($T_1 = 0.3~s$). The parameter $k$ refers to the number of pulses in the control process.}
	\label{fig:recognition_map}
\end{figure}
 %
We first analyze the ideal situation of a homogeneous ensemble of spin- 1/2 particles irradiated on resonance, which is described by Eq.~\eqref{eq:bloch} with $\omega=0$ and $\alpha=1$. We assume that $T_2=0.2$~s is perfectly known, the goal being to estimate the value of $T_1=0.3$~s. To simplify the presentation of the results, we consider a simple database associated with four values of $T_1$: 0.1, 0.233, 0.366 and 0.5~s. Following the general procedure of OFP, we compute the optimal field for this dictionary in the case where all the coefficients $\mu_{mn}$ are set to 1. The time $T$ is set to 10~ms. The efficiency of the optimal solution is benchmarked against a time-dependant random field as shown in Fig.~\ref{fig:recognition_map}, which displays the recognition map $(T_1(m),T_1(n))\mapsto D[f_m,f_n]$ for the two databases and the time evolution of the different elements of the dictionary. The contrast of Fig.~\ref{fig:recognition_map} gives a first quantitative measure of the precision of the recognition process. In this example, $C_N$ is equal to 0.06 for the optimal fields and 0.03 under the random fields. The minimum values of the recognition maps are respectively 0.019 and 0.001.

A first estimation of the value of $T_1$ can be made directly with the colorbars of Fig.~\ref{fig:recognition_map} and leads to $T_1\simeq 0.366$~s. Better accuracy of the fingerprinting process can be obtained by increasing the size of the dictionary. However, this procedure has a limit in terms of computational time, in particular to find the global optimum of the problem since the complexity of the control landscape increases rapidly with $N$. These numerical difficulties inherent to OFP can be avoided by using curve fitting in the post-measurement lookup stage. The fit is made with a minimization of $D$ based on a descent gradient algorithm with respect to the parameters $\vec{S}$ (here $T_1$ and $T_2$). In this case, the control field is fixed and a discrete derivative is used to compute the gradient $(\partial_{T_1}D[f(T_1,T_2),g],\partial_{T_2}D[f(T_1,T_2),g])$. Numerical simulations reveal that this approach converges after 50 or 100 iterations. Note that this concept is close to the Levenberg-Marquardt Method~\cite{Marquardt}, which is included in most of the curve fitting codes. For the ideal system, we obtain $T_1=0.3$~s both for the optimal and the random fields.

\begin{figure}[h]
\includegraphics[width=\linewidth]{fig4}
\caption{(Color online) Time evolution of the experimental signal (in the $y$- direction) during a sequence of 500 $\delta$- pulses. The signal is expressed in arbitrary units. The insert is a zoom showing the fluctuations of the signal. The parameter $k$ refers to the number of the pulse in the control process.}
	\label{fig:noise}
\end{figure}
We investigate the stability in presence of noise of this approach. An experimental example is displayed in Fig.~\ref{fig:noise}, where we observe the fluctuations of the signal around a mean value. The experimental setup is modeled by considering a simulated noise added to the response of the system:
\begin{equation}
g(t)=\bar{g}+\varepsilon\mathcal{N}(t),
\end{equation}
where $\bar{g}$ is the mean value of $g$ over many measurements, $\varepsilon$ the standard deviation and $\mathcal{N}$, a gaussian noise centered in 0 with a variance of 1. Since the radius of the Bloch ball is normalized to 1, $\varepsilon$ can be interpreted as a percent deviation. Using the optimal and random fields of Fig.~\ref{fig:recognition_map}, we optimize the parameter $T_1$ for different responses $g(t)$. The algorithm converges towards different values of $T_1$ for each response of the system. Figure~\ref{fig:uncertainty_fx_noise} displays the mean value and the width of the $T_1$- distribution (denoted $\Delta T_1^{\textrm{OPT}}$ and $\Delta T_1^{\textrm{RAND}}$ for the optimal and random fields respectively) as a function of $\varepsilon$. For each value of $\varepsilon$, we consider 30 measurements  $g(t)$ and the widths are determined by assuming a Gaussian distribution. This width can be interpreted as the accuracy of the corresponding estimation process. We observe in Fig.~\ref{fig:uncertainty_fx_noise} that the gain can be very large with the optimization procedure, a factor of the order of 100 for $\varepsilon=0.001$ is obtained. The random field fails to predict $T_1$ accurately, even for low noise amplitude. In a standard experiment, the amplitude of the noise is generally of the order of 1\% of the maximum of the signal. This correspond here to $\varepsilon = 0.01$. Similar results have been obtained for the parameter $T_2$.

\begin{figure}[h]
\includegraphics[width=\linewidth]{fig5}
\caption{(Color online) (top) Width of the distribution of the estimated $T_1$ values by the fingerprinting process (dark gray or blue - optimal, light gray or red - random) as a function of the noise amplitude $\varepsilon$, which is dimensionless. The dashed lines depict the mean values of the two distributions. The horizontal solid line is the value of the $T_1$ parameter. (bottom) Plot of the ratio $\Delta T_1^{\textrm{RAND}}/\Delta T_1^{\textrm{OPT}}$ of the width of the two distributions as a function of $\varepsilon$.}
	\label{fig:uncertainty_fx_noise}
\end{figure}

\section{Experimental results}\label{sec4}
We study experimentally the simultaneous estimation of the relaxation time $T_2$ and the distribution parameters of the offset inhomogeneities, while the $\alpha$- parameter can be set to 1 with good accuracy (see Eq.~\eqref{eq:bloch}). The offset distribution $\rho(\omega)$ is assumed to be Lorentzian:
\begin{equation}
\rho(\omega)\propto \left( 1+\frac{4(\omega - \bar \omega)^2}{\Delta \omega ^2}\right) ^{-1},
\end{equation}
where $\Delta \omega$ is the full width at half maximum (FWHM) and $\bar \omega$ the center of the distribution. The parameter $T_1$ was previously estimated to be 87~ms by inversion recovery \cite{levittbook}. The estimation of the parameter $T_2$ is a challenging issue because $T_2$ and $\Delta\omega$ are both responsible for the decay of the measured transverse magnetization. An effective transverse relaxation time $T_2^*$ defined by the relation
\begin{equation}\label{eqT2star}
\frac{1}{T_2^*}=\frac{1}{T_2} + \frac{\Delta \omega}{2}
\end{equation}
is usually introduced in magnetic resonance to account for the two physical effects \cite{levittbook}.
\begin{figure}
\includegraphics[width=\linewidth]{fig6}
\caption{(Color online) (top) Optimized control field along the $x$- direction ($\omega_y = 0$) to estimate the parameter $T_2$. (bottom) Time evolution of the simulated trajectories $M_x(t)$ (red or light gray) and $M_y(t)$ (blue or dark gray). The experimental data correspond to the gray areas around the numerical solutions, which give an estimation of the accuracy of the measurement (see the text for details).}
	\label{fig:experiment}
\end{figure}

A specific optimal pulse sequence sensitive to $T_2$ for an ensemble of spins with an average value of $\Delta\omega=20$~rad.s$^{-1}$ was designed. Note that only one control field along the $x$- direction was used to improve the convergence of the algorithm. Experiments were performed at room temperature on a Bruker Avance 600~MHz spectrometer. We used the $^1$H spins of H$_2$O with D$_2$O (99.9\%) as a solvent in a Shigemi tube. CuSO$_4$ was added as a $T_1$-shortening agent. The control field is a sequence of $N_p=500$ $\delta$-pulses separated by a time $T=10$~ms. The control field and the time evolution of the transverse magnetization are plotted in Fig.~\ref{fig:experiment}.
A reasonable match is found between the theoretical and the experimental results, which can be compared with the experimental error made in the measurement of the Bloch vector, as shown in Fig~\ref{fig:experiment}. Independent measurement based on a spin echo sequence leads to $T_2 = 60.5 \pm 0.5$~ms and $\Delta\omega=28.5 $~rad.s$^{-1}$. If we assume that the value of $\Delta \omega$ is known then OFP gives $T_2=60.4\pm 3.6$~ms and $\bar \omega = 0.1 \pm 0.6$~rad.s$^{-1}$. In the general case, due to the correlations between $\Delta \omega$ and $T_2$, it was not possible to estimate precisely the two parameters. As displayed in Fig.~\ref{fig:experiment2}, we observe that the figure of merit $D$ is almost the same for $\Delta\omega\in [20,38]$~rad.s$^{-1}$. On this interval, the value of $T_2^*$ is constant and in agreement with the experiment. Additional information would be required to estimate $T_2$ independently of $\Delta\omega$. From a computational point of view, it seems difficult to include different values of the bandwidth in the definition of the dictionary for improving the accuracy of the estimation. The same analysis was performed with several random sequences and we were not able to recover the right values of $T_2$ or $T_2^*$, showing thus the efficiency of OFP.

\begin{figure}
\includegraphics[width=\linewidth]{fig7}
\caption{(Color online) (Top) Minimum distance $D$ (defined in Eq.~\eqref{cost1}, $D$ is dimensionless) between the simulated and the experimental curves as a function of $\Delta \omega$. The dashed rectangle indicates the interval where $D$ is minimum. (Middle) Evolution of $T_2$ and $\bar{\omega}$ as a function of $\Delta\omega$. In the different panels, the blue (black) and the red/yellow (dark gray/ light gray) curves represent respectively the mean value of the signal and the upper/lower bounds of uncertainty. (Bottom) Plot of $1/T_2$ as a function of $\Delta \omega$. In the interval $[20,38]$~rad.s$^{-1}$, the parameter $T_2^*$ is constant as shown by the dotted line.} \label{fig:experiment2}
\end{figure}

\section{Conclusion}\label{conc}
We have introduced in this work the principles of OFP with an application to spin dynamics. The optimization procedure provides a method to approach the physical limits of the protocol in terms of sensitivity. FP has several advantages over the conventional methods. It allows a quantitative estimation of multiple parameters at the same time (e.g. the times $T_1$ and $T_2$), while only information about a single parameter is traditionally achieved. This advantage must be tempered if several parameters are correlated. This aspect has been illustrated in Sec.~\ref{sec4} with the offset terms and the $T_2$ relaxation time. The repeated acquisitions of data for the standard techniques are replaced by a single-shot measurement process in FP, which can drastically reduce the overall time of the experiment~\cite{naturemri}. Finally, FP is expected to be less sensitive to experimental imperfections and to the presence of noise. All these aspects are improved by the optimization procedure proposed in this paper. As shown in Sec.~\ref{sec3}, the better stability of OFP against noise perturbation is illustrated in a model example. A quantitative comparison with the standard Inversion Recovery approach~\cite{levittbook} is also made in App.~\ref{appc} for the estimation of the $T_1$ parameter. This theoretical comparison shows the better accuracy of OFP in this case.

This analysis paves the way for further investigations in MRI and realistic in vivo experiments~\cite{MRIbook,MRIconolly,lapertcontrast} in which the standard version of the fingerprinting process with random pulses has been applied with success~\cite{naturemri}. The concept of OFP could also be transferred to other domains such as quantum optics and atomic and molecular physics. An example is given by the control of molecular alignment and orientation in which pulse shaping techniques have been applied with success \cite{stapelfeldt:2003,seideman:2006}. The measure of the alignment could be used to estimate molecular parameters such as, e.g., the collisional relaxation rates \cite{ramakrishna:2005,vieillard:2013,vieillard:2008}. Another aspect could be to explore the applicability of this approach in a dynamical feedback framework where the control field would be adjusted in real time according to the results of the measurements. Finally, it seems promising to combine this technique with other methods of data analysis such as filter-diagonalization \cite{filter}, Bayesian estimation \cite{schirmer1} or Fisher  information~\cite{jacobpaper}.

\noindent\textbf{ACKNOWLEDGMENT.}\\
S.J. Glaser acknowledges support from the DFG (Gl 203/7-1). D. Sugny and S. J. Glaser acknowledge support from the ANR-DFG research program Explosys (ANR-14-CE35-0013-01; DFG-Gl 203/9-1). This work was done with the support of the Technische Universit\"at M\"unchen – Institute for Advanced Study, funded by the German Excellence Initiative and the European Union Seventh Framework Programme under grant agreement 291763. Experiments were performed at the
Bavarian NMR center at TU M\"unchen.