text:
\section{Introduction}\label{Intro}
Strong-field ionization of atoms has been proposed as a means
of estimating
 the peak intensity of  laser pulses having an intensity of $I>10^{20}$ W/cm$^2$ \cite{Hetzheim2009,BaukeHetzheimetal2011,Ciappinaetal2019,CiappinaPopruzhenko_2020,CiappinaBulanovetal2020}.
Such ultra-intense laser pulses are becoming available at several facilities worldwide, two examples being the ELI Beamlines \cite{ELIbeamlines2015}
and the Shanghai super-intense ultrafast laser facility \cite{Guo2018}.
In \cite{Ciappinaetal2019,CiappinaPopruzhenko_2020}, it was shown that the relative yields of different ionic charge states are highly sensitive to the peak intensity of the laser
pulse and that the final charge state distribution
 can be used to estimate the laser
field
 intensity at
the
beam focus. Moreover, it was concluded in \cite{CiappinaPopruzhenko_2020}
that the highest accuracy can be achieved in the determination of the laser
field intensity when the ionization yields of He- and H-like ions are adopted.

In order to
relate quantitatively
 the measured ion yield
 to
  the laser field intensity, it is necessary to simulate the ion yields using
  reliable
   formulas for the tunneling and
   over-the-barrier (OTB)
   static-field
    ionization rates.
 Ionization rates
for H-like ions are well established through the
tunneling-rate formula originally derived by Smirnov and Chibisov \cite{SMIRNOVCHIBISOV1966},
which has also been referred to as the
 Perelemov-Popov-Terent'ev (PPT) formula   \cite{PPT_1966,Popov_2004,V_Popruzhenko_2014} or the Ammosov-Delone-Krainov (ADK) \cite{Ammosov1986} formula,
providing
 reliable estimates of the ionization rates in the tunneling regime.
 Within the single-active electron approximation,
the ADK theory has been extended to molecular systems
\cite{TongZhaoLin2002,Kjeldsen_2004,Zhao2011,KornevZon2015,KornevChernovZon2017}.
Attempts have also been made to include many-electron effects in static-field tunneling ionization
\cite{FisherMaronPitaevskii1998,Brabecetal2005}.
When the field intensity is
sufficiently weak, the weak-field asymptotic theory (WFAT)
\cite{TolstikhinMorishitaMadsen2011} provides a rigorous expansion of the pre-exponential
factor in the tunneling rate in powers of the field strength.
 The
WFAT has been extended to many-electron systems \cite{TolstikhinMadsenMorishita2014,TolstikhinaMorishitaTolstikhin2014,Trinh_2015,Trinh_2016} and has been applied to molecular ionization \cite{YueBauchMadsen2017,DnestryanTolstikhinJensenMadsen2019}.
In \cite{Trinh_2016}, the WFAT  was applied to He-like ions for
an arbitrary value of
the nuclear charge number $Z$.

In order to estimate a static-field ionization rate
 in the  OTB regime,
 we need to adopt numerical methods because
  there are no analytical formulas available.
Static-field ionization
rates of
H \cite{NicolaidesThemelis1992,Karlsson_1992,TrinhTolstikhinMadsenMorishita2013},
He \cite{NicolaidesThemelis1993,ScrinziGeisslerBrabec1999,Scrinzi2000} and
H$_2$ \cite{Saenz2000} have been studied numerically by the complex scaling methods,
which have also been applied to
 many-electron atoms \cite{Jagau2016} and
 diatomic
 and triatomic molecules (CO, N$_2$, O$_2$, CO$_2$) \cite{Majety_2015,MajetyScrinzi2015}.
 In Ref.~\cite{ParkerArmstrongBocaTaylor2009}, static-field ionization rates of He were obtained by real-time integration of
 the time-dependent Schr\"odinger equation (TDSE).
 However, to the best of our knowledge, no numerical studies
 have been reported so far on the tunneling and OTB ionization rates of
  He-like ions.
  It has been awaited that reliable methods for estimating ionization rates for He-like ions valid in both the tunneling and the OTB regimes are established, so that the scheme for estimating the laser field intensity proposed in \cite{Ciappinaetal2019} will be realized. The advantage of treating He-like ions is that the effect of the electron-electron correlation in the static-field ionization processes can be varied systematically by the change in the nuclear charge number $Z$.

 In the present study, we
  investigate the static-field tunneling and OTB ionization of He-like ions with
  the
  nuclear charge number $Z$ in the range of $2\le Z\le 36$. We
  compare
    (i) ionization rates obtained by
  the PPT formula, (ii)
  ionization rates obtained by
   a PPT-based fitting formula proposed by Tong and Lin  \cite{TongLin_2005},
 (iii)
 ionization rates calculated by a perturbative expansion in powers of $1/Z$, (iv) ionization rates calculated by the single-active electron (SAE) approximation, and (v) reference
 ionization
 rates obtained by
 the  multiconfiguration time-dependent Hartree-Fock (MCTDHF) method \cite{ZanghelliniScrinzi2003,KatoKono2004}, in which
 the electron-electron interaction is fully accounted for.
We compare the  deviation of the
approximate rates from the rates calculated by the MCTDHF method, and find that both the
approach (ii) and
the
approach (iii) result in rates which deviate by about 10\% from the MCTDHF rates.
Considering that the approach (ii) is easier
to be applied than the approach (iii), we recommend the approach (ii) to be used in
 the evaluation of static-field ionization rates of He-like ions.

\section{Theory}\label{Sec:Theory}
In this section, we describe the theoretical methods
that we employ in the present study.
Atomic units (a.u.) are used throughout the
paper unless
 otherwise indicated. For the conversion from a.u.\ to SI units, we have $1\;{\rm a.u.} = 4.134\times 10^{16}\;{\rm s}^{-1}$ for ionization rates, and $1\;{\rm a.u.} = 5.142\times 10^{11}$ V/m for electric fields.

In
 the calculation of the static-field ionization rates,
we employ the nonrelativistic   TDSE. Although we will present
results for the ionization of He-like ions with
the large charge number $Z$ up to $Z = 36$, where the intensity of a laser field
required for ionization reaches
 $10^{22}$ W/cm$^2$, it has been shown that the ionization process itself is essentially  nonrelativistic
for ions with $Z<20$ \cite{MilosevicKrainovBrabec2002,Ciappinaetal2019}.
As long as the ionization potential $I_p$ of an atom is much smaller than the rest energy of an electron,
$I_p\ll c^2$, where $c\approx 137$ a.u.,
the relativistic
effects in the ground state wave function can be neglected.
Using $I_p\approx Z^2/2$ for large $Z$, we see that $I_p\ll c^2$ is satisfied even for $Z=36$ because $I_p/c^2=36^2/(2c^2)\approx 0.03$.

The static-field ionization  rates obtained in this study can be used to  estimate the
ionization probability in a low-frequency laser field
as long as the static-field ionization rates can be  used as
instantaneous ionization rates.
This procedure is valid if the Keldysh
parameter $\gamma=\sqrt{2I_p}\omega/\Efield$, where $\omega$ is the
angular frequency of the laser field and $\Efield$ is the peak field strength, is
smaller than 1 \cite{Keldysh1964,Ciappinaetal2019}.
In this study,
we consider values of $\Efield$ so that $\gamma<1$.
For example,  for He ($Z=2$), we consider $\Efield\approx 0.4$ a.u., which
gives $\gamma=0.2$ assuming $\omega=0.057$ a.u.\ (corresponding to a wavelength of 800 nm), and for He-like Kr ($Z=36$), we consider $\Efield\approx 2\times 10^3$ a.u.,
which gives
$\gamma=0.001$ at $\omega=0.057$ a.u.

\subsection{$Z$-scaling}\label{subsec:Zscaling}
The TDSE for a He-like ion with
the
 charge number $Z$ in a static electric field $\bm{\Efield}=\Efield\hat{\bm{z}}$ reads
\begin{equation}
\label{TDSE_unscaled}
i\frac{\partial}{\partial t}\Psi(\bm{r}_1,\bm{r}_2,t)=
\left[H+\Efield(z_1+z_2)\right]
\Psi(\bm{r}_1,\bm{r}_2,t),
\end{equation}
where $\bm{r}_i=(x_i,y_i,z_i)$ is the coordinate of the $i$th electron, and  the field-free Hamiltonian is defined as
\begin{equation}
H=\sum_{i=1,2}\left(-\frac{1}{2}\frac{\partial^2}{\partial\bm{r}^2_i} -\frac{Z}{r_i} \right)
+\frac{1}{r_{12}},
\end{equation}
with $r_{12}=|\bm{r}_1-\bm{r}_2|$.
 After scaling  the spatial coordinates and time according to \cite{ScherrKnight1963}
\begin{equation}
\bm{r}_i=\frac{2\bm{r}'_i}{Z},\qquad t = \frac{4t'}{Z^2},
\end{equation}
and rewriting the TDSE \eqref{TDSE_unscaled} in the scaled coordinates, we obtain
\begin{equation}
\label{TDSE_scaled}
i\frac{\partial}{\partial t'}\Psi'(\bm{r}'_1,\bm{r}'_2,t')=
\left[H'+\Efield'(z'_1+z'_2)\right]
\Psi'(\bm{r}'_1,\bm{r}'_2,t'),
\end{equation}
where
\begin{equation}\label{scaled_H0}
H'=\sum_{i=1,2}\left(-\frac{1}{2}\frac{\partial^2}{\partial\bm{r}'^2_i} -\frac{2}{r'_i} \right)
+\frac{2}{Z}\frac{1}{r'_{12}},
\end{equation}
and $\Psi'(\bm{r}'_1,\bm{r}'_2,t')=(2/Z)^3\Psi(\bm{r}_1,\bm{r}_2,t)$.
The scaling of the amplitude of the electric  field
is
\begin{equation}\label{FieldScaling}
\Efield'=\frac{8}{Z^3}\Efield.
\end{equation}

The transformation described above implies that a He-like atomic ion with
the
charge number $Z$ is described by the same TDSE as He
and that
 the electron-electron interaction is scaled by $2/Z$,
 indicating that the relative magnitude  of the electron-electron interaction
 compared to the electron-nucleus interaction decreases as $Z$ increases.

For
the ionization rate,
 we have the following scaling relation,
\begin{equation}\label{scalingrelation_rates}
\W(\Efield)=\frac{Z^2}{4}\W'(\Efield'),
\end{equation}
where $\W(\Efield)$ is the rate
we obtain
 by solving Eq.~\eqref{TDSE_unscaled} at the field
 $\Efield$ and
$\W'(\Efield')$ is the rate
we obtain
 by solving Eq.~\eqref{TDSE_scaled} at the field $\Efield'= 8 \Efield/Z^3$.

In the following part of the paper,  symbols with a prime
 always refer to
the quantities in
the transformed system \eqref{TDSE_scaled}
and symbols without a prime
  refer to
those in
the original system \eqref{TDSE_unscaled}.

\subsection{MCTDHF}\label{subsec:MCTDHF}
In order to solve the  scaled TDSE \eqref{TDSE_scaled} numerically, we employ the MCTDHF method
\cite{ZanghelliniScrinzi2003,KatoKono2004,HochstuhlHinzBonitz2014
IshikawaSato_Review2015,Loetstedt2017,SatoOrimoetal2018}. In the MCTDHF method, the  two-electron  wave function is expanded in terms   of time-dependent Slater determinants
$\|\phi'_i(t')\alpha\,\phi'_j(t')\beta\|$,
\begin{equation}\label{MCTDHF_def}
\Psi'(t')=\sum_{i,j=1}^M C'_{ij}(t')\|\phi'_i(t')\alpha\,\phi'_j(t')\beta\|,
\end{equation}
where  $\phi'_i(\bm{r}', t')$ is a time-dependent spatial orbital, $M$ is the total number of orbitals, $\alpha$ (spin up) and
$\beta$ (spin down) are spin functions, and $C'_{ij}(t')$ is a time-dependent  configuration-interaction (CI) coefficient. The equations of motion for the spin-orbitals and the CI coefficients are derived by applying the time-dependent variational principle \cite{ZanghelliniScrinzi2003,KatoKono2004}.
The MCTDHF method has previously been applied to He
in the calculation of
two-photon ionization probabilities \cite{Hochstuhl2011}, high-harmonic spectra \cite{SatoIshikawaetal2016}, above-threshold photoelectron spectra \cite{Orimoetal2018},
and strong-field excitation probabilities \cite{LotstedtSzidarovszkyetal2020}.

It should be noted here
 that the MCTDHF method is different from other commonly employed close-coupling approaches
\cite{Zhang_1995,SMYTH1998,LaulanBachau2003,ParkerArmstrongBocaTaylor2009} in that products of
single-electron orbitals are used
in the expansion of
 the time-dependent wave function, while  in \cite{Zhang_1995,SMYTH1998,LaulanBachau2003,ParkerArmstrongBocaTaylor2009} the two-electron wave function is expanded
in terms of two-electron radial basis functions multiplied with two-electron spherical harmonics.

We solve the MCTDHF equations of motion in the same way as described in \cite{LotstedtKatoYamanouchi2017,LotstedtSzidarovszkyetal2020}.
The details of the numerical implementation of the MCTDHF method are
summarized in
Appendix ~\ref{App_MCTDHF}.

In order to obtain the
ionization
 rate $\W$, we solve
the scaled TDSE \eqref{TDSE_scaled} starting from the ground state, $\Psi'(t'=0)=\Psi'_0$, and calculate the time-dependent population of the ground state, $p_0(t')=|\langle \Psi'_0|\Psi'(t)\rangle|^2$.
For
$t' < t'_0$
($t'_0\approx 50$ a.u.),
 $p_0(t')$ oscillates due to the transitions to
the
excited states induced by the sudden turn-on of the
static
field. For $t'>t'_0$, $p_0(t')$ decreases exponentially
and can be fitted  to an exponential function $e^{-\W' t' -q'}$  with
the
constant  $\W'$, representing
 the scaled ionization rate,
 and $q'$ \cite{Scrinzi2000}.
 After the simulations with
  different values of the
  static
  field $F'$, the scaled field-dependent ionization rate $\W'(\Efield')$
  is obtained, from which  the
 ionization rate of a He-like ion is
 calculated
 by the scaling of Eq.~\eqref{scalingrelation_rates}.
 However, if the ionization rate is very large or very small, it becomes difficult to obtain the ionization rate by this approach because $p_0(t')$ cannot be fitted well to
 $e^{-\W' t' -q'}$ when it changes too slowly or too rapidly.
Consequently, the scaled fields need to be limited practically in the range of $0.1\lesssim\Efield' \lesssim 1$ a.u.

\subsection{Perturbation theory}\label{subsec:PertTheory}
In the limit
of
 $Z\to\infty$, the scaled Eq.~\eqref{TDSE_scaled} decouples and becomes an equation for two independent He$^+$ ions. Therefore, in the large  $Z$ limit,  the scaled ionization rate for a He-like ion can be expressed as
\begin{equation}\label{ExactScaling_WH}
\W'(\Efield')\underset{Z\to \infty}{=}8\WH(\Efield),
\end{equation}
where $\WH(\Efield)$ is the ionization rate of H evaluated at
the
unscaled field of $\Efield=\Efield'/8$.

For finite but large values of $Z$, we can treat the electron-electron interaction term $2/(Zr'_{12})$ in the scaled Hamiltonian \eqref{scaled_H0} as a perturbation, and  evaluate corrections
using  perturbation theory. By inserting the expansions
\begin{equation}
\Psi'(\bm{r}'_1,\bm{r}'_2)=\sum_{n=0}^\infty \frac{\chi'_n(\bm{r}'_1,\bm{r}'_2)}{Z^n}
\end{equation}
and
\begin{equation}
\epsilon'=\sum_{n=0}^\infty \frac{\epsilon'_n}{Z^n}
\end{equation}
into the
eigen equation,
\begin{equation}
H'(e^{i\Theta}\bm{r}'_1,e^{i\Theta}\bm{r}'_2)\Psi'(\bm{r}'_1,\bm{r}'_2)=\epsilon' \Psi'(\bm{r}'_1,\bm{r}'_2),
\end{equation}
 where $H'(e^{i\Theta}\bm{r}'_1,e^{i\Theta}\bm{r}'_2)$ is the
 scaled complex-rotated
  Hamiltonian  \eqref{scaled_H0}  \cite{Reinhardt1982}, we derive   equations relating $\chi'_n$ and $\epsilon'_n$ at different orders of $1/Z$.   Throughout our calculations we use
the complex rotation parameter $\Theta =0.2$. The equations up to
the third order
 are given in Appendix~\ref{PertTheoryAppendix_subsec1}.
 After the
field-dependent correction terms $\epsilon'_n(\Efield')$
for
 the eigenvalue of the ground state
are calculated,
 the unscaled ionization rate
 corrected up to
 the
 order $N$ at a given value of $Z$ is given by
 \begin{equation}\label{perttheory_W}
 \Wpert^N(\Efield') = \frac{Z^2}{4}\sum_{n=0}^N\frac{\lambda'_n(\Efield')}{Z^n},
 \end{equation}
 with
 \begin{equation}\label{pertCoeffs}
\lambda'_n(\Efield') =-2 \imag\left[\epsilon'_n(\Efield')\right].
 \end{equation}
 It should be noted here  that the correction terms $\epsilon'_n$ contain corrections both to the eigenenergy (the real part of $\epsilon'_n$) and to the ionization rate (the imaginary part of $\epsilon'_n$) of the ground state
 and that the rate given by Eq.~\eqref{perttheory_W} contains contributions
from all orders in the field $\Efield'$, because the interaction with the field is included in the zeroth-order Hamiltonian.

\subsection{Single-active electron approximation}\label{subsec:SAE}
In the single-active electron approximation (SAE), the SAE Hamiltonian reads
\begin{equation}\label{SAEHamiltonian}
H'_{\rm SAE}(\bm{r}')=-\frac{1}{2}\frac{\partial^2}{\partial\bm{r}'^2} +V(\bm{r}') +\Efield' z',
\end{equation}
where $V(\bm{r}')$ is an effective potential.
We consider below the SAE approximation to the scaled TDSE \eqref{TDSE_scaled}.

A commonly used effective potential for He-like ions reads \cite{Scrinzi2000}
\begin{equation}\label{V1}
V_{1}(\bm{r}')=-\frac{2-\frac{2}{Z}+\frac{2}{Z}e^{-\zeta' r'}}{r'},
\end{equation}
where the parameter $\zeta'$ is adjusted so that the absolute value of the ground state energy
becomes equal to the ionization potential of the corresponding two-electron system.
For He, we additionally consider a polarization potential in order to model the effect of the electric field on the inner electron,
that is, the electron left bound after the ejection of the other electron. The polarization potential is thus defined as
\begin{equation}\label{Vpol}
V_{\rm pol}(\bm{r}) = -\Efield \alpha_{{\rm He}^+} \cos\theta \frac{1-e^{-\eta r}\left(1+\eta r+\frac{1}{2}\eta^2 r^2\right)}{r^2},
\end{equation}
where $\theta$ is the polar angle, $\alpha_{{\rm He}^+}=9/2^5\approx 0.28$ a.u.\ is the polarizability of He$^+$ \cite{Kotani1951},  and $\eta$ is a  parameter
describing the range of the polarization potential.
 The polarization potential
is constructed
so
 that $V_{\rm pol}(r\to\infty)=-\Efield \alpha_{{\rm He}^+} \cos\theta/r^2$ and $V_{\rm pol}(r\to0)=0$.
Finally, the full effective SAE potential including $V_{\rm pol}(\bm{r})$ reads
\begin{equation}\label{V2}
V_2(\bm{r})=V_{1}(\bm{r})+V_{\rm pol}(\bm{r}).
\end{equation}
The parameter $\eta$ is fixed
so that two
 times the polarizability of the single active electron matches
 the polarizability  of He.

For He, we have also investigated a more complicated  SAE  potential of the form
$V_3(r)=-[1+e^{-\zeta_1r}+\beta r(r-q)e^{-\zeta_2r}]/r$,
where the parameters $\zeta_1$, $\zeta_2$, $\beta$, and $q$ are determined by the minimization of  the difference between the ground state single-electron density obtained by the SAE approximation and the ground state single-electron density obtained by the MCTDHF method.
The difference between
 the ionization rates obtained by using
the potential $V_3(r)$ and the  ionization rates obtained by  the potential $V_1(r)$ defined in Eq.~\eqref{V1} was found to be small.
Indeed, the relative difference
is smaller than 4\% for fields $0.2\le \Efield\le 0.6$ a.u. We therefore conclude that the ionization rate is not sensitive to
a small variation of the form of the SAE potential as long as the ionization potential is well reproduced.

The ionization rates in the SAE approximation are  obtained by calculating the
lowest energy
 eigenvalue $\epsilon'_{\rm SAE}$ of  the complex-rotated Hamiltonian
$H'_{\rm SAE}(e^{i\Theta}\bm{r}')$ \cite{Reinhardt1982},
\begin{equation}\label{SAE_eigenvalueEq}
H'_{\rm SAE}(e^{i\Theta}\bm{r}')\psi'_{\rm SAE}(\bm{r}') = \epsilon'_{\rm SAE} \psi'_{\rm SAE}(\bm{r}'),
\end{equation}
with $\Theta =0.2$. From the complex-valued eigenvalue $\epsilon'_{\rm SAE}$,
we obtain the scaled tunneling rate for a two-electron atom as
\begin{equation}\label{WSAE}
\WSAE' = -4\imag(\epsilon'_{\rm SAE}),
\end{equation}
and the
unscaled original
  rate by the scaling relation \eqref{scalingrelation_rates}.
Note that there is an additional factor of 2 in Eq.~\eqref{WSAE}
originating from
 the indistinguishability of the two electrons.

In addition to the SAE approximation,
we also
 discuss results obtained by the frozen core (FC) model of He \cite{Scrinzi2000},
\begin{align}\label{SAEfrozencoreHe1s2}
\Psi_{\rm FC} (\bm{r}_1,&\bm{r}_2,t)=
c_0(t)\Psi^{\rm He}_0(\bm{r}_1,\bm{r}_2)
\nonumber
\\
&+\frac{1}{\sqrt{2}}\left[\psi_0^{{\rm He}^+}(\bm{r}_1)\phi(\bm{r}_2,t)+\phi(\bm{r}_1,t)\psi_0^{{\rm He}^+}(\bm{r}_2)\right],
\end{align}
 where $\Psi^{\rm He}_0(\bm{r}_1,\bm{r}_2)$ is the ground state wave function of He, and $\psi_0^{{\rm He}^+}(\bm{r})$ is the ground state wave function of He$^+$.
We can derive equations
 of motion for $c_0(t)$ and $\phi(\bm{r},t)$  by inserting $\Psi_{\rm FC} (\bm{r}_1,\bm{r}_2,t)$ into the TDSE \eqref{TDSE_unscaled}.
  We obtain the ionization rate $\WFC$ by
 fitting the decreasing ground-state norm $|c_0(t)|^2$ to an exponential function $e^{-\WFC t-q}$.
 In the context of
 static-field ionization of He, the ansatz \eqref{SAEfrozencoreHe1s2} was introduced by Scrinzi \cite{Scrinzi2000}, and similar models have been successfully used
 in the numerical simulations of
   many-electron molecules in strong fields \cite{SpannerPatchkovskii2009,SPANNER201310,MajetyScrinzi2015,MajetyZielinskiScrinzi2015}.
 The wave function \eqref{SAEfrozencoreHe1s2} describes a single active electron with the ejected part of the wave function given by $\phi(\bm{r},t)$,
 and includes the exchange interaction
 of the ejected electron with the He$^+$ core,
 differently from the SAE model.
\section{Results}\label{Sec:Results}
\subsection{Comparison of MCTDHF rates and PPT rates}

\begin{figure}
\includegraphics[width=\figwidth]{Fig1}
\caption{\label{Fig1}Ionization rates $\W=Z^2\W'/4$ calculated with the MCTDHF method  for (a) He, (b) He-like Be ($Z=4$), and (c) He-like Ar ($Z=18$) are shown with filled circles. The ionization rates are shown as a function of the scaled field $F'=8F/Z^3$.
Comparisons are made with the
 ionization rates calculated
 using
  the PPT formula [solid line; see Eq.~\eqref{PPTrate}]
  and
   the PPT-based fitting formula
suggested in \cite{TongLin_2005} [dotted line; see Eq.~\eqref{modPPT}]. In panel (a), the dash-dotted line shows the results obtained by Scrinzi, Geissler and Brabec in
\cite{ScrinziGeisslerBrabec1999}. The scaled barrier suppression field $F'_{\rm BS}$ is indicated with a
vertical solid line. The blue shaded area indicates the interval where the ionization probability during one laser cycle
$P_{\rm ion}$
[defined in Eq.~\eqref{Pioncycle}]
is in the range between  10\% and 90\%.
In terms of peak  laser field intensities, $F'=0.4$ a.u.\ corresponds to $5.6\times 10^{15}$ W/cm$^2$ for $Z=2$,
$3.6\times 10^{17}$ W/cm$^2$ for $Z=4$, and $3.0\times 10^{21}$ W/cm$^2$ for $Z=18$.}
\end{figure}

In Fig.~\ref{Fig1}, we show the comparison
between
the
ionization rates $\W$ calculated with the MCTDHF method and the ionization rates $\WPPT$ calculated with the PPT formula
\cite{PPT_1966,V_Popruzhenko_2014,CiappinaPopruzhenko_2020}. The PPT formula
employed in this study is given by   Eq.~\eqref{PPTrate} in Appendix ~\ref{PPTformula}. Note that the PPT rate
is defined with an additional  factor of 2,
originating from
 the fact that two electrons contribute to the total ionization rate.
For He
($Z=2$),  a comparison is made with the rates $\W_{\rm SGB}$ calculated by Scrinzi, Geissler and Brabec \cite{ScrinziGeisslerBrabec1999}
using the complex scaling method for a correlated two-electron wave function.

Also shown in Fig.~\ref{Fig1} are
the results obtained by the
PPT-based fitting  formula for the ionization rate
  suggested by Tong and Lin \cite{TongLin_2005},
\begin{equation}\label{modPPT}
\WTong(\Efield) = \WPPT(\Efield) \Lambda_\alpha(F) ,
\end{equation}
where
\begin{equation}\label{empirical_correctionPPT}
\Lambda_\alpha(\Efield)=\exp\left[-\alpha\frac{(Z-1)^2}{I_p}\frac{\Efield}{(2I_p)^{\frac{3}{2}}}\right]
\end{equation}
is an empirical correction factor, and $I_p$ is the ionization potential. For  the $I_p$ of He-like ions, we use the expression given in Eq.~\eqref{Ip_def} in Appendix \ref{PPTformula}.
We note that the expression \eqref{modPPT} should be viewed as a formula with which the ionization rates obtained
by ab initio methods can be fitted.
In \cite{TongLin_2005}, the parameter $\alpha$ in the  empirical correction factor
$\Lambda_\alpha(F)$  was suggested to be
$\alpha=7$ for He, but the values for other He-like ions were not given. We derive $\alpha$ by a
least-squares fit of the
PPT-based fitting formula $\WTong(\Efield)$ to the rate obtained by the MCTDHF method, using the ionization rates in the range
of  $0.2\le \Efield'\le 0.4$ a.u.
We obtain $\alpha=7.12$ for $Z=2$, $\alpha=  7.14$ for $Z=4$, and $\alpha= 5.62$ for $Z=18$.
The
applicability  of the PPT-based fitting formula \eqref{modPPT} has been confirmed
previously in the case of   strong-field ionization of H,
where
 the PPT-based fitting formula \eqref{modPPT} with $\alpha=6$ differs from accurate numerical
ionization rates  by less  than 10\% for
the unscaled fields of
 $\Efield<0.125$ a.u.\ \cite{TongLin_2005}.

Because $\Lambda_\alpha(\Efield)$
becomes exponentially small
as $\Efield$ increases,
the formula \eqref{modPPT}
cannot be valid when $\Efield$ becomes too large,
 but, as can be seen in Fig.~\ref{Fig1}, it is remarkably accurate
as long as the field
strength  is in the tunneling regime as well as in the OBT regime below the field strength
at
$\Efield\approx2\EBS$, where the barrier-suppression field $\EBS$ is
defined as
\begin{equation}\label{EBS_def}
\EBS = \frac{I_p^2}{4(Z-1)}.
\end{equation}
The quantity $Z-1$
represents
 the residual charge
with which  the ejected electron interacts.
We regard that the field strength is in the tunneling regime when
$\Efield<\EBS$ and that the field strength is in the OTB regime when $\Efield>\EBS$.
We show $\EBS$ in Fig.~\ref{Fig1} as a thick vertical line.
 The scaled barrier-suppression field $\Efield'_{\rm BS}=8\Efield_{\rm BS}/Z^3$ is defined according to the scaling relation, Eq.~\eqref{FieldScaling}.

In \cite{ZhangLanLu2014}, an improved
empirical correction factor
$\Lambda_{\alpha_1\alpha_2\alpha_3}(\Efield) =
\exp(-\alpha_1\Efield^2/\EBS^2-\alpha_2\Efield/\EBS-\alpha_3)$ to the PPT formula with the three fitting  parameters $\alpha_1$, $\alpha_2$, $\alpha_3$  was proposed.
However, we
 employ the original correction factor introduced in \cite{TongLin_2005},
 in which only one fitting parameter, $\alpha$, is used, for  its  simplicity.

In order to relate  the static-field rates
to the probability of ionization in an oscillating laser field, we evaluate the ionization probability
$P_{\rm ion}$ during one optical cycle of a laser field with
the
peak amplitude $\Efield$ and the
angular frequency $\omega$ as
\begin{equation}\label{Pioncycle}
P_{\rm ion}(\Efield)=1-\exp\left[-\frac{1}{\omega}\int_0^{2\pi} \W(|\Efield \sin s|) \rmd s  \right],
\end{equation}
where $\omega=0.057$ a.u., corresponding to a wavelength of 800 nm,
is chosen.
The expression \eqref{Pioncycle} for the ionization probability of  a He-like ion in a laser field  is valid under the assumption  of  a slowly varying field.
In Fig.~\ref{Fig1}, the interval where $P_{\rm ion}$ fulfills $0.1\le P_{\rm ion}\le 0.9$ is indicated by
the area shaded with
light blue color.

We can see in Fig.~\ref{Fig1} that the standard PPT
formula overestimates the ionization rates compared with the reference MCTDHF rates, even though the extent of the overestimation becomes smaller
 for
 the
 larger values of $Z$.
 On the other hand, the rates obtained by the  PPT-based fitting formula agree well with the MCTDHF rates.
For $Z=2$, the relative
deviation
 takes the values of
  $|\WTong-\W|/\W<0.1$
  when   $\Efield'\le 0.4$ a.u.\ and $|\WTong-\W|/\W<0.5$
  when
    $\Efield'\ge 0.4$ a.u.
 For $Z=4$, we have  $|\WTong-\W|/\W<0.1$
 when
  $\Efield'\ge 0.25$ a.u. For $Z=18$, $|\WTong-\W|/\W<0.1$ when $0.25\le \Efield'\le 0.5$ a.u.\ and
 $|\WTong-\W|/\W<0.15$ when $\Efield'> 0.5$ a.u.
 This shows that the PPT-based fitting formula has a relative accuracy of 10\%
 or better in the electric field  range where the ionization probability is substantially large
 but below
 the
 saturation
 ($0.1<P_{\rm ion}<0.9$), corresponding to
  the range indicated by light blue color in Fig.~\ref{Fig1}.
 An accurate formula for the ionization rate in this
 field strength
  range is crucial
 for the simulation of the ion yields required in the intensity estimation scheme proposed in
  \cite{Ciappinaetal2019,CiappinaPopruzhenko_2020}.

 In the case of  He, the MCTDHF rates and the rates $\W_{\rm SGB}$ calculated by Scrinzi, Geissler, and Brabec \cite{ScrinziGeisslerBrabec1999} are in good agreement with each other,
 which confirms the validity of the MCTDHF method for calculating static-field ionization rates. The relative
 deviation takes the values of
$|\W_{\rm SGB}-\W|/\W_{\rm SGB}<0.05$ when $\Efield'\le 0.5$ a.u.\ and $|\W_{\rm SGB}-\W|/\W_{\rm SGB}<0.08$ when $\Efield'>0.5$ a.u.
The small discrepancy between our results and the  results obtained in \cite{ScrinziGeisslerBrabec1999} can be ascribed to the difference in the method used
in the calculations.
In \cite{ScrinziGeisslerBrabec1999},  the ionization rates were obtained by the complex scaling
of
 a two-electron wave function expanded in terms of explicitly correlated basis functions.
The rates for He obtained by the MCTDHF method are also in good agreement with the
ionization rates reported by Parker et al.\ in \cite{ParkerArmstrongBocaTaylor2009}, where the TDSE was integrated in real time and the ionization rate was obtained by  monitoring the decay of the norm of the wave function in a sphere
whose radius is 18 a.u.\ around the atomic core. For comparison,
numerical values of the ionization rates for He obtained by the MCTDHF method are listed
together with the ionization rates
reported in  Ref.~\cite{ScrinziGeisslerBrabec1999} and
Ref.~\cite{ParkerArmstrongBocaTaylor2009} in Appendix \ref{Sec:Herates_comparison}.

\begin{figure}
\includegraphics[width=\figwidth]{Fig2}
\caption{\label{Fig2}(a) Ionization rate $\W=Z^2\W'/4$ at $\Efield=\EBS$ calculated with the MCTDHF method (blue solid circles), the rate $\WTong$ calculated with the
PPT-based fitting formula (red triangles) and a fit of $\W$ to a quadratic function $\afit Z^2$ with $\afit=1.4\times 10^{-3}$ (thick solid gray line)  are shown as a function of $Z$. (b) Ionization rates evaluated at $\Efield=(2I_p)^{3/2}\times 0.05$, such that the reduced field $\Fred=0.05$ a.u. The coefficient $\afit$ in the quadratic fit is $\afit=1.3\times 10^{-4}$.
}
\end{figure}

To further assess the accuracy of the PPT-based fitting formula, we evaluate the ionization rate as a function of $Z$ at two
different
values of $\Efield$.
In Fig.~\ref{Fig2}(a), we show the ionization rate evaluated at the barrier-suppression field $\EBS$  for
 $Z\le 36$. We compare the rate $\W(\Efield=\EBS)$
 calculated with the MCTDHF method and the rate $\WTong(\Efield=\EBS)$ calculated with the PPT-based fitting formula  defined
  by Eq.~\eqref{modPPT}.  For each value of $Z$, $\alpha$ is obtained by fitting
  $\WTong(\Efield')$ to the ionization rates
 $\W(\Efield')$ calculated by the MCTDHF method  for
 the  scaled fields in the range
 of
  $0.2\le\Efield'\le0.4$ a.u. The numerical values of $\alpha$ are given in
 Table~\ref{alphaTable}.
 For  reference, we provide
 numerical values of the ionization rates in the Supplemental material
 \footnote{See Supplemental
 material at [link inserted by PRA] for a text file containing the static-field ionization rates obtained by the MCTDHF method
 for $2\le Z\le 36$}.

 The relative difference between $\W$ and $\WTong$,
 defined as $|\WTong(Z)-\W(Z)|/\W(Z)$,
  is less than 6\% for all $Z\ge2$.
  It can also be seen in Fig.~\ref{Fig2}(a)
  that   the $Z$-dependence of the ionization rate $\W$ at $\Efield=\EBS$
  is quadratic  to a good approximation.
 This means that the
 scaled rate $\W'(\EBS')$ is almost independent of $Z$ because
 $\W$ is scaled by $Z^2/4$ with respect to $\W'$ as shown in the relation \eqref{scalingrelation_rates}.
  A fit of $\W$ to the function $\W_{\rm fit}=\afit Z^2$
 gives $\afit \approx 1.4\times 10^{-3}$.
 As can be seen  in Fig.~\ref{Fig2}(a),  $\W_{\rm fit}=\afit Z^2$ agrees well with the rate obtained by the MCTDHF method.
 For $Z$ in the range $2\le Z\le10$,
 the relative error of the quadratic fit is
 $|\W_{\rm fit}(Z)-\W(Z)|/\W(Z)<0.15$, and for $Z>10$,  it becomes $|\W_{\rm fit}(Z)-\W(Z)|/\W(Z)<0.05$.

 \begin{table}
\caption{\label{alphaTable}Parameter $\alpha$
adopted  in the PPT-based fitting formula \eqref{modPPT}.}
\begin{ruledtabular}
\begin{tabular}{m{\colw}m{\colw} | m{\colw}m{\colw} | m{\colw}m{\colw} | m{\colw}m{\colw}}
$Z$ & $\alpha$ & $Z$ & $\alpha$ & $Z$ & $\alpha$ & $Z$ & $\alpha$\\
\hline
    2&      7.12&        8  &     6.23 &        20&      5.57&       32&      5.41\\
    3&      7.52&        10&      6.00  &      22&      5.53&       34&       5.40\\
    4&      7.14&        12 &     5.88  &      24 &     5.49&       36&      5.39\\
    5&      6.84&        14 &     5.77  &      26 &     5.47 &      &         \\
    6&      6.58&        16  &    5.68  &      28 &     5.45 &       &         \\
    7&      6.38&        18  &    5.62 &       30  &    5.43 &        &         \\
\end{tabular}
\end{ruledtabular}
\end{table}

 In Fig.~\ref{Fig2}(b), we show the ionization rates evaluated at
 $\Efield=(2I_p)^{3/2}\times0.05$ so
  that the reduced field $\Fred$, defined as
 \begin{equation}\label{defFred}
\Fred=\frac{\Efield}{(2 I_p)^{3/2}},
\end{equation}
equals $0.05$ a.u. The reduced field $\Fred$ appears in the PPT rate as
$\exp[-2/(3\Fred)]$ [see Eq.~\eqref{PPTrate}], and therefore,
 the PPT rate  varies rapidly with $\Fred$ in the tunneling regime. The particular value of
$\Fred=0.05$ a.u.\ was empirically proposed in  \cite{Ciappinaetal2019} as a threshold for  saturated ionization, beyond which atoms are completely ionized. The condition
of
$\Fred=0.05$ a.u.\ means that the
field strength is  in the tunneling regime rather than in the over-barrier regime.
For example,
for $Z=2$, $\Efield=(2I_p)^{3/2}\times 0.05=0.12\;{\rm a.u.}<\EBS=0.20$ a.u., and, for $Z=10$,
$\Efield=(2I_p)^{3/2}\times 0.05=41\;{\rm a.u.}<\EBS=54$ a.u.
As shown in Fig.~\ref{Fig2}(b), the PPT-based fitting formula works well also at the reduced  field  of $\Fred=0.05$ a.u.
 Indeed, the deviations are sufficiently small in the entire range of $Z$, that is,
$|\WTong-\W|/\W<0.12$ for  $Z\le 12$ and $|\WTong-\W|/\W<0.03$ for  $Z> 12$.
However, it should be noted that the quadratic function fits well only
in
the large $Z$ above
$\sim 20$, that is, (i) $|\W_{\rm fit}-\W|/\W>0.4$ for $Z<10$,
(ii) $0.4>|\W_{\rm fit}-\W|/\W>0.1$ for $10\le Z<20$, and
(iii) $|\W_{\rm fit}-\W|/\W<0.1$ for $Z\ge20$.

\subsection{Ionization rates by perturbation theory}

\begin{figure}
\includegraphics[width=\figwidth]{Fig3}
\caption{\label{Fig3}Ionization rates  for (a) $Z=10$ and (b) $Z=18$ obtained by the MCTDHF method
(filled circles) and those obtained by perturbation theory as
 a function of the
scaled field $\Efield'=8\Efield/Z^3$. The perturbative rate
$\Wpert^N$ contains corrections up to
the
order $1/Z^N$. In panel (c), we show the rates as a function of
$Z$ at
 a fixed value $\Efield'=0.4$ a.u.\ of the scaled field.}
\end{figure}

In Fig.~\ref{Fig3}, we compare the ionization rates $\W$ calculated using the MCTDHF method and
 the ionization rates $\Wpert^N$ [see Eq.~\eqref{perttheory_W}] calculated by perturbation theory. The numerical values of the expansion coefficients $\lambda'_n$
 given by Eq.~\eqref{pertCoeffs}
  are shown in Appendix \ref{PertTheoryAppendix_subsec2}.
Our numerical approach for solving the
equations of perturbation theory,
 described in Appendix \ref{PertTheoryAppendix_subsec2}, limits
the range of the
 scaled field strengths to $\Efield'\ge 0.25$ a.u.\ for
 the orders of $n=0$ and 1, and to
$\Efield'\ge 0.35$ a.u.\ for the orders of $2$ and 3.

We can see in Figs.~\ref{Fig3}(a) and (b) that the rate $\Wpert^{N=3}$ including  the third-order corrections agrees well with the rate calculated with the
MCTDHF method.
For both $Z=10$ and $Z=18$, the relative
deviations are as small as  $|\Wpert^{N=3}-\W|/\W<0.07$ in the scaled field range of
 $\Efield'\ge 0.4$ a.u.
The values of $Z=10$ and $Z=18$ are selected to demonstrate the performance of
perturbation theory for
the
large values of $Z\ge 10$, where
the perturbation  to the third order
  is enough to
obtain ionization rates with a sufficiently small deviation from the reference MCTDHF rates.
It can be seen in Figs.~\ref{Fig3}(a) and (b) that
  perturbation theory gives relatively accurate results for the ionization rate in both the tunneling and the OTB regimes (for $Z=10$, $\EBS'=0.43$ a.u.,
and for $Z=18$, $\EBS'=0.46$ a.u.). In the tunneling regime, where the ionization rate is exponentially sensitive to
the variation in   the ionization potential, it is necessary to include
sufficiently large orders in perturbation theory
so
 that an accurate value of the ground state energy is obtained. The results shown in Figs.~\ref{Fig3}(a) and (b) suggest that $N=3$ is large enough to obtain ionization rates
with a relative deviation from the MCTDHF rates below 7\%  both in the tunneling and the OTB regimes.

In Fig.~\ref{Fig3}(c), we show the ionization rates as a function of $Z$, evaluated
at the scaled field of  $\Efield'=0.4$ a.u. The deviation of the ionization rate calculated by the third order perturbation theory
becomes
 smaller as $Z$ increases, that is,
$|\Wpert^{N=3}-\W|/\W=0.29$ for $Z=2$,
$|\Wpert^{N=3}-\W|/\W=0.09$ for $Z=8$, and
 $|\Wpert^{N=3}-\W|/\W<0.07$ for $Z\ge 10$. However,  the third-order rate
 $\Wpert^{N=3}$ exhibits  a minimum in the rate
at  $Z=5$ in a similar manner as  the MCTDHF rate exhibiting a minimum at $Z=6$,
as can be seen in the curves plotted with filled circles and filled squares in Fig.~\ref{Fig3}(c).
This minimum in the ionization rate appearing at a certain $Z$ value can be explained by  the asymptotic behavior of the
PPT rate. At the fixed scaled field of $\Efield'=8\Efield/Z^3=0.4$ a.u., the PPT rate increases as $\WPPT\propto Z^2$ for large $Z$  ($Z>30$).  For the smaller values of $Z\le 7$, where the ionization rate is in the tunneling regime, $\WPPT$ decreases as $Z$  increases because  the  reduced field
$\Fred\approx \frac{\Efield'}{8} \left(1-\frac{5}{4Z}\right)^{-\frac 3 2}$ decreases as
$Z$ increases. Therefore, in the intermediate range, there must be a minimum in
the ionization rate, which is expected to appear  approximately at the value of $Z$ where
$\Efield' \approx \EBS'$. In the case of $ \Efield'=0.4$ a.u.,  $\Efield'<\EBS'$ when $Z\le 7$.

\subsection{Comparison of SAE approximation rates and MCTDHF rates}\label{subsec:SAE_MCTDHF_comp}
\begin{figure}
\includegraphics[width=\figwidth]{Fig4}
\caption{\label{Fig4}Ionization rates $\W=Z^2\W'/4$ as a function of $\Efield'$ for (a) He and (b)  He-like Ne ($Z=10$), calculated with the MCTDHF method
and  the
 SAE model, using potential $V_{1}$ [see Eq.~\eqref{V1}]. In (a) we additionally show the rates obtained with the SAE model including the polarization potential
[see Eq.~\eqref{V2}], and the FC  model including a correlated He ground state wave function [see Eq.~\eqref{SAEfrozencoreHe1s2}]
In (c), we show the relative
deviation $\Delta\W/\W =(\WSAE-\W)/\W$ of the rate $\WSAE$
 obtained with the SAE approximation  from the MCTDHF rate $\W$
  as a function of $Z$,
 evaluated at $\Efield'=0.4$ a.u.\
 and
  $\Efield=\EBS$. Also shown are fits to the function $b/Z$, with $b(\Efield'=0.4)=0.68$, and $b(\Efield=\EBS)=0.84$. In the case of
 $\Efield=\EBS$, only $\Delta\W/\W$ values for $Z\ge 10$ were included in
 the fit. }
\end{figure}

In Fig.~\ref{Fig4}, we show
the comparison between the
results obtained by the SAE approximation
described in Sec.~\ref{subsec:SAE} and those obtained by the MCTDHF method. The numerical values of the  parameter $\zeta'$  are
$\zeta'=2.1325$ for $Z=2$ and  $\zeta'=2.4489$ for $Z=10$.
For the polarization potential defined in Eq.~\eqref{Vpol}, we employ the value  $\eta=4.8713$ for He ($Z=2$),
so that  twice the polarizability $\alpha_{\rm SAE}$ of
the active electron equals the polarizability
of He, $\alpha_{\rm He}=1.3832$ a.u.\ \cite{Yangetal2017}.  The numerical calculation of the SAE rates and the  calculation of the  SAE polarizability are described in
Appendix \ref{AppendixSAE}.

We can see in Fig.~\ref{Fig4}(a) that the SAE rate
for He is larger than  the MCTDHF rate by around 30\% (31\% for $\Efield'= 0.4$ a.u.\ and 27\% for $\Efield'= 0.5$ a.u.).
The agreement is not significantly improved when the  polarization potential is
included in the  effective SAE potential, as defined in Eq.~\eqref{V2}.
 On the other hand, the results obtained by the  FC model defined in  Eq.~\eqref{SAEfrozencoreHe1s2}
agree very well with the MCTDHF results
($|\WFC-\W|/\W<0.06$ for $\Efield'\le 0.4$ a.u.).
Because  the SAE model
should give  the correct ionization potential
in its  construction and the
FC model,
in which the correlated ground state is included  in the wave function ansatz \eqref{SAEfrozencoreHe1s2}, should also give
 the correct ionization potential, the better agreement
 of the ionization rate given by
  the FC model cannot be attributed to  the difference in
   the ionization potential.
   In addition, because the
 core polarization is not included in the FC model, the polarization of the He$^+$
 core
 cannot be responsible
  for the difference between the SAE results and the FC results. The difference between the FC model and the SAE model
 can be ascribed to
  the exchange interaction between the ejected electron and the He$^+$ core
  that is not included in the SAE model but in the FC model.
  Indeed, the FC wave function \eqref{SAEfrozencoreHe1s2} is symmetric with respect to
  the exchange of the two electron coordinates. Therefore, the disagreement between the SAE rate and the MCTDHF rate
  is considered to originate from the exchange interaction.
  The overestimation of the ionization rate given by the SAE approximation can be ascribed to the increase in the ionization probability originating from the absence of the exchange interaction. The importance of
   the exchange interaction
   was reported before
    in the alignment dependence of the tunneling ionization rates in
    CO$_2$, O$_2$, and N$_2$ \cite{MajetyScrinzi2015}.
In order to further prove that the exchange interaction can affect the ionization process, we have included in Appendix \ref{AppendixSAE_HF} a discussion
on two different frozen-core Hartree-Fock models, that is,
one in which the exchange symmetry is fulfilled and the other in which no exchange symmetry is considered.

Figure~\ref{Fig4}(b), in which the  MCTDHF rate and the SAE rate
are compared for $Z=10$, shows that the relative difference is less than 10\%
in the range of $0.25\le\Efield'\le 0.6$ a.u.
In the case of $Z=18$ (not shown in Fig.~\ref{Fig4}), the relative
difference is smaller than 6\% for $0.25\le\Efield'\le 0.6$ a.u.,
 suggesting that the
effect of the exchange interaction decreases
as  $Z$ increases.

In order to quantify the effect of the exchange interaction, we define the
relative deviation of the SAE rate from the MCTDHF rate
as
\begin{equation}
\frac{\Delta\W}{\W}=\frac{\WSAE-\W}{\W},
\end{equation}
where $\W$ denotes the ionization rate calculated by the MCTDHF method,
and show $\Delta\W/\W$ as a function of $Z$ in Fig.~\ref{Fig4}(c). We evaluate the relative
deviation at two different field strengths, $\Efield'=0.4$ a.u., and $\Efield=\EBS$.
In both cases, the relative deviation $\Delta\W/\W$
decreases approximately as $1/Z$ as $Z$ increases, and the
decrease  proportional to $1/Z$ is shown clearly
 by the best fit curves having  the form of $b/Z$ with a constant factor of $b$.
 The deviation of the SAE rate from the MCTDHF rate caused by the exchange interaction is approximately 30\% in He, and
 the deviation
 decreases as $1/Z$ to about 2\% in He-like Kr.

\section{Summary and conclusions}
We have performed an
investigation of the static-field  ionization rates of He-like ions for
the
nuclear charge numbers $Z$ up to $Z=36$, in both the tunneling and OTB regimes. We have
made two major findings. The first is that both
the PPT-based fitting formula by Tong and Lin \cite{TongLin_2005} [see Eq.~\eqref{modPPT} and Table~\ref{alphaTable}] and the third-order perturbation theory in $1/Z$ [see Eq.~\eqref{perttheory_W}]
give small deviations from the rates calculated by the MCTDHF method. In the case of
the PPT-based fitting formula,
the
relative deviation from the reference rates obtained by the MCTDHF method is
 about 10\%
 when  $Z>2$ and  $\Efield\le0.075 Z^3$ a.u.,
 corresponding to scaled fields $\Efield'\le 0.6$ a.u.
  In the case of the ionization rates calculated with third-order perturbation theory,
 the relative deviation from the MCTDHF rates is 7\%
  for  $Z\ge 10$.
 In practice, the PPT-based fitting formula with the
 parameter $\alpha$ whose values are  given  in Table~\ref{alphaTable}
 provides reasonably good estimates of the ionization rate and will help
  researchers
 estimate the ion yields in
  the intensity measurement scheme proposed in \cite{Ciappinaetal2019,CiappinaPopruzhenko_2020}.

 The second finding is that the exchange interaction in
 the
 tunneling and OTB barrier ionization reduces the ionization rate by an amount proportional to $1/Z$, ranging from a reduction of  30\% for He
 ($Z=2$) to about 2\% for He-like Kr
 ($Z=36$).
 This demonstrates the importance of the exchange interaction
 in the
 strong-field ionization of He-like ions, and shows how the relative strength of the exchange
 interaction decreases as the nuclear charge number
 $Z$ increases.

\begin{acknowledgments}
This research was supported by  JSPS KAKENHI grants no.\ JP18K05024 and
no.\ JP15H05696,
a JSPS Invitational Fellowship for Research in Japan (no.\ S18108),
 and a travel grant  from the Czech Academy of Sciences (no.\ VAJVA-19-04).
M. F. C. acknowledges the project Advanced research using high intensity laser produced photons and particles (CZ.02.1.01/0.0/0.0/ 16\_019/0000789) from European Regional Development Fund (ADONIS), the Spanish Ministry MINECO (National Plan 15 Grant: FISICATEAMO No.\ FIS2016-79508-P, SEVERO OCHOA No.\ SEV-2015-0522, FPI), European Social Fund, Fundaci\'o Cellex, Generalitat de Catalunya (AGAUR Grant No.\ 2017 SGR 1341 and CERCA/ Program), ERC AdG NOQIA, and
the National Science Centre, Poland-Symfonia Grant No.\ 2016/20/W/ST4/ 00314.
\end{acknowledgments}


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

Ancillaries are loved or hated, accepted or rejected, but typically
ignored. Recent approximate ancillary methods (e.g., \citep{reid}) give a decomposition of the sample space rather than providing
statistics on the sample space (e.g., \citep{cox1,mcc1}). As a result, continuity gives the contour along which the
variable directly measures the parameter and then gives the subcontour
that provides measurement of a parameter of interest. This, in turn,
enables the high accuracy of cumulant generating function
approximations \citep{dan,bar1} to extend to cover a wide
generality of statistical models.

Ancillaries initially arose (see \citep{fis1}) to examine the accuracy
of the maximum likelihood estimate, then (see \citep{fis2}) to
calibrate the loss of information in the use of the maximum likelihood
estimate and then (see \citep{fis3}) to develop a key instance
involving the configuration statistic. The configuration of a sample
arises naturally in the context of sampling a location-scale model,
where a standardized coordinate $z = (y - \mu)/\sigma$ has a fixed and
known error distribution $g(z)$: the $i$th coordinate of the response
thus has $f(y_i; \mu, \sigma) = \sigma^{-1} g\{(y_i - \mu)/\sigma\}$.
The configuration $a(y)$ of the sample is the plug-in estimate of the
standardized residual,
\begin{equation}
\label{hat}
a(y) = \hat z =  \biggl({y_1 - \hat\mu\over\hat\sigma}, \ldots, {y_n -
\hat\mu\over\hat\sigma}  \biggr)^\prime,
\end{equation}
where $(\hat\mu, \hat\sigma)$ is the maximum likelihood value for $(\mu
, \sigma)$ or is some location-scale equivalent. Clearly, the
distribution of $\hat z$ is free of $\mu$ and $\sigma$ as the
substitution $y_i = \mu+ \sigma z_i$ in (\ref{hat}) leads to the
cancellation of dependence on $\mu$ and $\sigma$. This supports a
common definition for an ancillary statistic $a(y)$, that it has a
parameter-free distribution; other conditions are often added to seek
sensible results.

More generally, the observed value of an ancillary identifies a sample
space contour along which parameter change modifies the model, thus
yielding the conditional model on the observed contour as the
appropriate model for the data. The ancillary method is to use directly
this conditional model identified by the data.

One approach to statistical inference is to use only the observed
likelihood function $L^0(\theta) = L(\theta; y^0)$ from the model $f(y;
\theta)$ with observed data $y^0$. Inference can then be based on some
simple characteristic of that likelihood. Alternatively, a weight
function $w(\theta)$ can be applied and the composite $w(\theta)L(\theta
)$ treated as a distribution describing the unknown $\theta$; this
leads to a rich methodology for exploring data, usually, but
unfortunately, promoted solely within the Bayesian framework.

A more incisive approach derives from an enriched model which is often
available and appropriate. While the commonly cited model is just a set
of probability distributions on the sample space, an enriched model can
specifically include continuity of the model density function and
continuity of coordinate distribution functions. An approach that
builds on these enrichments can then, for example, examine the observed
data $y^0$ in relation to other data points that have a similar shape
of likelihood and are thus comparable, and can do even more. For the
location-scale model, such points are identified by the configuration
statistic; then, accordingly, the model for inference would be $f\{y
\mid
a(y) = a^0 ;\theta\}$, where $a(y)$ is the configuration ancillary.

Exact ancillaries as just described are rather rare and seem limited to
location-type models and simple variants. However, extensions that use
approximate ancillaries
(e.g., \citep{fra5,fra9}) have recently been broadly fruitful,
providing approximation in an asymptotic sense. Technical issues can
arise with approximate values for an increasing number of coordinates,
but these can be managed by using ancillary contours rather than
statistics; thus, for a circle, we use explicitly a contour
$A = \{(x,y)= (a^{1/2}\cos t, a^{1/2}\sin t) \dvt  t \mbox{ in }[0, 2\curpi)\}$
rather than using implicitly a statistic $x^2 + y^2 = a$.

We now assume independent coordinate distribution functions that are
continuously differentiable with respect to the variable and the
parameter; extensions will be discussed separately. Then, rather than
working directly with a coordinate distribution function $u_i = F_i
(y_i; \theta)$, we will use the inverse, the quantile function $y_i =
y_i(u_i; \theta)$ which presents a data value $y_i$ in terms of a
corresponding $p$-value $u_i$. For additional advantage, we could use a
scoring variable $x$ in place of the $p$-value, for example, $x = \Phi
^{-1}(u)$ or $x = F^{-1}(u; \theta_0)$, where $\Phi(\cdot)$ is the
standard Normal distribution function. We can then write $y = y(x;
\theta)$, where a coordinate $y_i$ is presented in terms of the
corresponding scoring variable $x_i$.

For the full response variable, let
$
y = y(x; \theta) =
\{y_1(x_1; \theta),\ldots, y_n(x_n; \theta) \}'
$
be the quantile vector expressing $y$ in terms of the reference or
scoring variable $x$ with its given distribution: the quantile vector
records how parameter change affects the response variable and its
distribution, as prescribed by the continuity of the coordinate
distribution functions.

For an observed data point $y^0$, a convenient reference value $\hat
x^0$ or the fitted $p$-value vector is obtained by solving the equation
$ y^0 = y(x; \hat\theta^0)$
for $x$, where $\hat\theta^0$ is the observed maximum likelihood value;
for this, we assume regularity and asymptotic properties for the
statistical model.
The contour of the second order ancillary through the observed data
point as developed in this paper is then given
as the trajectory of the reference value,
\begin{equation}\label{obs}
A^0= \{ y(\hat x^0;t) \dvt  t \mbox{ in } \mathbb{R}^p\},
\end{equation}
to second order under parameter change,
where $p$ here is the dimension of the parameter. A~sample space point
on this contour has, to second order, the same estimated $p$-value
vector as the observed data point and special
properties for the contours are available to second order.

The choice of the reference variable with given data has no effect on
the contour: the reference variable could be Uniform, as with the
$p$-value; or, it could be the response distribution itself for some
choice of the parameter, say $\theta_0$.

For the location-scale example mentioned earlier, we have the
coordinate quantile function $y_i = \mu+ \sigma z_i$, where $z_i$ has
the distribution $g(z)$. The vector quantile function is
\begin{equation} \label{quan}
y(z; \mu, \sigma) = \mu1+\sigma z ,
\end{equation}
where $1 = (1, \ldots, 1)^\prime$ is the `one vector.' With the data
point $y^0$, we then have the fitted $\hat z^0 =(y^0 - \hat\mu^0 1)/
\hat\sigma^0$. The observed ancillary contour to second order is then
obtained from (\ref{obs}) by substituting $\hat z^0$ in the quantile
(\ref{quan}):
\begin{equation} \label{lsanc}
A^0 = \{y(\hat z^0; t)\} = \{m1 + s\hat z^0; (m, s) \mbox{ in } \mathbb{R}\times\mathbb{R}^+\} =\mathcal{L}^+(1; \hat z^0)
\end{equation}
with positive coefficient for the second vector.
This is the familiar exact ancillary contour $a(y)=a^0$ from (\ref{hat}).

An advantage of the vector quantile function in the context of the
enriched model mentioned above
is that it allows us to examine how parameter change modifies the
distribution and thus how it moves data points as a direct expression
of the explicit continuity.
In this sense, we define the velocity vector or vectors as
$v(x; \theta) = (\partial/ \partial\theta)y(x; \theta) =\partial y/
\partial\theta$.
In the scalar $\theta$ case, this is a vector recording the direction
of movement of a point $y$ under $\theta$ change; in the vector $\theta
$ case, it is a $1 \times p$ array of such vectors in $\mathbb{R}^n$,
$V(x; \theta) = \{v_1(x_1; \theta), \ldots, v_p(x_p; \theta)\},
$
recording the separate effects from the parameter coordinates $\theta
_1, \ldots, \theta_p$.
For the location-scale example, the velocity array is
$V(z; \mu, \sigma) = (1, z),$
which can be viewed as a $1 \times2$ array of vectors in $\mathbb{R}^n$.

The ancillary contour can then be presented using a Taylor series about
$y^0$ with coefficients given by the velocity and acceleration $V$ and $W$.
For the location-scale example,
the related acceleration vectors are equal to zero.

For more insight, consider the general scalar $\theta$ case and the
velocity vector $v(x; \hat\theta^0)$. For a typical coordinate, this
gives the change $\mathrm{d}y=v(x; \hat\theta^0)\,\mathrm{d}\theta$ in the variable as
produced by a small change $\mathrm{d}\theta$ at $\hat\theta^0$. A re-expression
of the coordinate variable can make these increments equal and produce
a location model; the product of these location models is a full
location model $g(y_1 - \theta, \ldots, y_n - \theta)$ that precisely
agrees with the initial model to first derivative at $\theta= \hat
\theta^0$ (see \citep{fra7,and}). This location model then, in
turn, determines a full location ancillary with configuration $a(y) =
(y_1 - \bar y, \ldots, y_n - \bar y)$. For the original model, this
configuration statistic has first-derivative ancillarity at $\theta=
\hat\theta^0$ and is thus a first order approximate ancillary; the
tangent to the contour at the data point is just the vector $v(\hat
x^0; \hat\theta^0)$. Also this contour can be modified to give second
order ancillarity.

In a somewhat different way, the velocity vector $v(y^0; \theta)$ at
the data point $y^0$ gives information as to how data change at $y^0$
relates to parameter change at various $\theta$ values of interest.
This allows us to examine how a sample space direction at the data
point relates to estimated $p$-value and local likelihood function
shape at various $\theta$
values; this, in turn, leads to quite general default priors for
Bayesian analysis (see \citep{fra8}).

In the presence of a cumulant generating function, the saddle-point
method has produced highly accurate third order approximations for
density functions (see \citep{dan}) and for distribution functions (see
\citep{lug}). Such approximations are available in the presence of
exact ancillaries \citep{bar1} and extend widely in the presence of
approximate ancillaries (see \citep{fra5}). For third order accuracy,
only second order approximate ancillaries are needed, and for such
ancillaries, only the tangents to the ancillary contour at the data
point are needed (see \citep{fra5,fra6}). With this as our
imperative, we develop the second order ancillary for statistical inference.

Tangent vectors to an ancillary at a data point
give information as mentioned above concerning a location model approximation
at the data point. For a scalar parameter, these provide a vector field
and integrate quite generally to give a unique approximate ancillary to
second order
accuracy. The resulting
conditional model then provides definitive $p$-values by available
theory; see, for example,
\citep{fra9}.
For a vector parameter, however, the multiple vector fields may not satisfy
the Frobenius conditions for integrability and thus may not define a function.

Under mild conditions, however, we show that such tangent vectors do generate
a surface to second order without the Frobenius conditions holding.
We show this in several steps. First, we obtain the coordinate quantile
functions
$
y_i = y_i(x_i; \theta).
$
Second, we Taylor series expand the full vector quantile
$y=(y_1,\ldots,y_n)$ in terms
of the full reference variable $x=(x_1,\ldots,x_n)$ and the parameter
$\theta=(\theta_1,\ldots,\theta_p)$ about data-based values, appropriately
re-expressing coordinates and working to second order.
Third, we
show that this generates a partition with second order ancillary properties
and the usual tangent vectors. The seeming need for the full Frobenius
conditions is bypassed
by finding that two integration routes need not converge to each other, but
do remain on the same contour, calculating, of course, to second order.

This construction of an approximate ancillary is illustrated in Section~\ref{sec2} using the familiar example, the Normal-on-the
circle from \citep{fis4}; see also \citep{cox2,bar2,fra7,fra3}. The example, of course, does have an exact
ancillary and the present procedure gives an approximation to that ancillary.
In Section~\ref{examples}, we consider various examples that have exact and
approximate ancillaries, and then in Sections~\ref{sec4} and~\ref{sec5}, we present the
supporting theory.
In particular, in Section~\ref{sec4}, we develop notation for a $p$-dimensional
contour in $\mathbb{R}^n$,
$
A= \{y(x_0;t)\dvt t \mbox{ in } \mathbb{R}^p\},
$
and use velocity and acceleration vectors to present a Taylor series
with respect to $t$. Then, in Section~\ref{sec5}, we consider a regular
statistical model with asymptotic properties and
use the notation from Section~\ref{sec4} to develop the second order ancillary contour
through an observed data point $y^0$. The re-expression
of individual coordinates, both of the variable and the parameter,
plays an essential role in the development; an asymptotic analysis is
used to establish the second order approximate ancillarity.
Section~\ref{sec6} contains some discussion.

\section{Normal-on-the-circle}\label{sec2}

We illustrate the second order approximate ancillary with a simple
nonlinear regression model,
the \textit{Normal-on-the-circle} example (see \citep{fis4}). The model
has a well-known exact ancillary.
Let $y = (y_1, y_2)^\prime$ be Normal on the plane with
mean $(\rho\cos\theta, \rho\sin\theta)'$ and variance matrix
$I/n$ with $\rho$ known. The mean is
on a circle of fixed radius $\rho$ and the distribution has rotationally
symmetric error with variances $n^{-1}$, suggesting an
antecedent sample size $n$ for an asymptotic approach. The full
$n$-dimensional case is examined as Example~\ref{ex2} in
Section~\ref{examples} and the present case derives by routine conditioning.

The distribution is a unit probability mass centered at $(\rho\cos\theta, \rho\sin\theta)'$ on the circle with radius $\rho$. If
rotations about the origin are applied to $(y_1, y_2)'$, then the
probability mass rotates about the origin, the mean moves on the circle
with radius $\rho$ and an element of
probability at a distance $r$ from the origin moves on a circle of radius $r$.
The fact that the rotations move
probability along circles but not between circles of course implies
that probability on any circle about the origin
remains constant: probability flows on the ancillary contours.
Accordingly, we have that the radial distance $r = (y_1^2
+ y_2^2)^{1/2}$ has a fixed $\theta$-free distribution and is thus ancillary.

The statistic $r(y)$ is the Fisher exact ancillary for this problem and
Fisher recommended
that inference be based on the conditional model, given the observed
ancillary contour. This conditional approach has a long but uneven
history; \citep{fra4} provides an overview and \citep{fra10} offer
links with asymptotic theory.
We develop the approximate second order ancillary and examine how it
relates to the Fisher exact ancillary.

The model for the Normal-on-the-circle has independent coordinates, so
we can invert the coordinate
distribution functions and obtain the vector
quantile function,
\[
\pmatrix{
y_1\cr
y_2
}
=\rho
\pmatrix{
\cos \theta\cr
\sin \theta
}
+
\pmatrix{
x_1\cr
x_2
}
,
\]
where the $x_i=\Phi^{-1}(u_i)/n^{1/2}$ are independent normal variables
with means $0$ and variances $n^{-1}$, and $\Phi$ is the standard
Normal distribution function.
We now examine the second order ancillary contour $A^0$ given by (\ref{obs}).

Let $y^0 = (y_1^0, y_2^0)' = (r^0\cos a^0, r^0\sin a^0)$ be the
observed data point where $r^0$, $a^0$ are the corresponding polar
coordinates; see Figure~\ref{circle}. For this simple nonlinear normal
regression model, $\hat\theta^0 = a^0$ is the angular direction of the
data point. The fitted reference value $\hat x^0$ is the solution of
the equation
$y^0 = y(x; \hat\theta^0) =\rho(\cos a^0, \sin a^0)' + (x_1, x_2)$,
giving $\hat x^0 = (\hat x_1^0, \hat x_2^0)'= y^0 - \rho(\cos a^0,
\sin a^0)' =y^0 - \hat y^0$,
where $\hat y^0=\rho(\cos a^0, \sin a^0)'$ is the fitted value, which
is the projection of the data point $y^0$ onto the circle. The observed
ancillary contour is then
\[
A^0 =  \left\{\rho
\pmatrix{
\cos\theta\cr
\sin\theta
}
+ y^0 - \hat y^0\dvt  \theta \mbox{ near } a^0\right \}
= y^0 -\hat y^0+  \left\{\rho
\pmatrix{
\cos(a^0 + t)\cr
\sin (a^0 + t)
}
\dvt  t \mbox{ near } 0\right \}.
\]
Figure~\ref{circle} shows that $A^0 = \{y(\hat x^0; t)\dvt  t \mbox{ near }a^0\}$ is a translation, as shown by the arrow of a segment $S$ of the
solution contour, from the fitted point $\hat y^0$ to the data point
$y^0$.
\begin{figure}

\includegraphics{248f01.eps}

\caption{The regression surface $S$ is a circle of radius $R$; the
local contour of
the approximate ancillary $A^0$ is a circle segment of $S$ moved from
$\hat y^0$ to $y^0$;
the exact ancillary contour is a circle segment of radius $r^0$ through
the data point $y^0$.}
\label{circle}
\end{figure}

The second order ancillary segment at $y^0$ does not lie on the exact
ancillary surface $r(y_1, y_2)=r^0$. The tangent vector at the data
point $y^0$ is
$v = (\partial y / \partial t)|_{t = a^0} =
(-\rho\sin a^0, \rho\cos a^0)'$, which
is the same as the tangent vector for the exact ancillary and which
agrees with
the usual tangent vector $v$ (see \citep{fra9}). However, the
acceleration vector
is $w =(\partial^2 / \partial t^2)y|_{t = a^0} = (-\rho\sin a^0,
-\rho\cos a^0)'$,
which differs slightly from that for the exact ancillary: the
approximation has radius of curvature $\rho$, as opposed to $r^0$ for
the exact, but the difference in moderate deviations about $y^0$ can be
seen to be small and is second order.

The second order ancillary contour through $y^0$ can also be expressed
in a Taylor series as $A^0=\{y^0 + tv+wt^2/2\dvt  t \mbox{ near } 0\}$;
here, the acceleration vector $w$ is
orthogonal to the velocity vector $v$.
Similar results hold in wide generality when $y$ has dimension $n$ and
$ \theta$ has dimension $p$; further examples are discussed in the next
section and the general development follows in Sections~\ref{sec4} and~\ref{sec5}.

\section{Some examples}\label{examples}

\begin{example}[(Nonlinear regression, $\sigma_0$ known)]\label{ex1}
Consider a nonlinear regression model $y=\eta(\theta) + x$ in $\mathbb{R}^n$,
where the error $x$
is $\operatorname{Normal}(0;\sigma^2_0 I)$ and the regression or solution surface $S=\{
\eta(\theta)\}$ is smooth with
parameter $\theta$ of dimension, say, $r$. For given data point $y^0$,
let $\hat\theta^0$ be the
maximum likelihood value. The fitted value is then $\hat y^0=\eta(\hat
\theta^0)$ and the fitted reference value is $\hat x^0=y^0 - \eta(\hat
\theta^0) =y^0-\hat y^0$. The model as presented is already in quantile
form; accordingly,
$V=(\partial\eta/ \partial\theta)|_{\hat\theta^0}, W=(\partial^2\eta
/\partial\theta^2)|_{\hat\theta^0}$
are the observed velocity and acceleration arrays, respectively, and
the approximate ancillary
contour at the data point $y^0$ is
$
A^0=\{y^0 + Vt+ t'Wt/2+\cdots\dvt t \mbox{ in } \mathbb{R}^r\},
$
which is just a $y^0-\hat y^0$ translation of the solution surface
$
S=\{\hat y^0 + Vt+ t'Wt/2+\cdots\dvt t \mbox{ in }  \mathbb{R}^r\}.
$
For this, we use matrix multiplication to linearly combine the elements
in the arrays $V$ and $W$.
\end{example}

\begin{example}[(Nonlinear regression, circle case)]\label{ex2}
As a special case, consider the regression model where the solution surface
$S=\{\eta(\theta)\}$ is a circle of radius $\rho$ about the origin;
this is the full-dimension version
of the example in Section~\ref{sec2}. For notation, let $C=(c_1,\ldots, c_n)$ be an
orthonormal
basis with vectors $c_1, c_2$ defining the plane that includes $S$. Then
$\tilde y=C'y$ provides rotated coordinates and $\tilde\eta(\theta
)=C'\eta(\theta)
=(\rho\cos\theta, \rho\sin\theta,0,\ldots,0)$ gives
the solution surface in the new coordinates.

There is an exact ancillary given by $r=(\tilde y_1^2+\tilde
y_2^2)^{1/2}$ and $(\tilde
y_3,\ldots, \tilde y_n)$; the corresponding ancillary contour through
$\tilde y^0$ is a
circle
of radius $r^0$ through the data point $y^0$ and lying in
the plane $\tilde y_3=\tilde y_3^0,\ldots,\tilde y_n=\tilde y_n^0$. The
approximate
ancillary contour is a segment of a circle of radius $\rho$ through the
data point $y^0$ and
lying in
the same plane. This directly agrees with the simple
Normal-on-the-circle example
of Section~\ref{sec2}.

For the nonlinear regression model, Severini (\citep{sev}, page 216)
proposes an approximate ancillary by using the obvious pivot $y-\eta
(\theta)$ with the plug-in maximum likelihood value
$\theta=\hat\theta$; we show that
this gives a statistic $ A(y)=y-\eta(\hat\theta)$ that can be
misleading. In the rotated coordinates, the statistic $ A(y)$ becomes
\begin{eqnarray*}
\tilde A(y)&=&(r\cos\hat\theta, r\sin\hat\theta,\tilde y_3,\ldots,\tilde
y_n)'-(\rho\cos\hat\theta, \rho\sin\hat\theta,0,\ldots,0)'\\
&=& \{(r-\rho)\cos\hat\theta, (r-\rho)\sin\hat\theta,\tilde y_3,\ldots
,\tilde
y_n\}',
\end{eqnarray*}
which has observed value $\tilde A^0 =\{(r^0-\rho)\cos\hat\theta^0,
(r^0-\rho)\sin\hat\theta^0,\tilde
y_3^0,\ldots,\tilde y_n^0\}'.$

If we now set the proposed ancillary equal to its observed value,
$\tilde A=\tilde A^0$, we obtain $\tilde y_3=\tilde y_3^0,\ldots, \tilde
y_n= \tilde
y_n^0$
and also obtain $r=r^0$ and $\hat\theta=\hat\theta^0$. Together, these
say that
$y=y^0$, and
thus that the proposed approximate ancillary is exactly equivalent to the
original response variable, which is clearly not ancillary. Severini
does note
``\ldots it does not necessarily follow that $a$ is a second-order
ancillary statistic since
the dimension of $a$ increases with $n$.'' The consequences of using the plug-in
$\hat\theta$ in the pivot are somewhat more serious: the plug-in
pivotal approach for this example does
not give an approximate ancillary.
\end{example}
\begin{example}[(Nonlinear regression, $\sigma$ unknown)]\label{ex3}
Consider a nonlinear regression model $y=\eta(\theta) + \sigma z$ in
$\mathbb{R}^n$, where
the error $z$ is $\operatorname{Normal}(0;I)$ and the solution surface $S=\{\eta
(\theta)\}$ is
smooth with surface dimension $r$ (see \citep{fra11}). Let $y^0$ be the
observed data point and
$(\hat\theta^0, \hat\sigma^0)$ be the corresponding maximum likelihood
value. We then have the fitted regression $\hat y^0$, the fitted residual
$\hat x^0 = y^0 -\hat y^0$, and the fitted reference value $\hat
z^0=\hat x^0/\hat\sigma^0$
which is just the standardized
residual.

Simple calculation gives the velocity and acceleration arrays
\[
\bar V=(V   \hat z^0),\qquad
\bar W=
\pmatrix{
W & 0 \cr
0 & 0
}
\]
using $V$ and $W$ from Example~\ref{ex1}. The approximate ancillary contour at
the data point
$y^0$ is then
\begin{eqnarray*}
\tilde A^0 &=& \{y^0 + VT +t'Wt/2 + \cdots+s\hat z^0\dvt  t  \mbox{ in } \mathbb{R}^r, s \mbox{ in } \mathbb{R}^+ \} \\
&=& \{\eta(t) + s \hat z^0 \dvt  t \mbox{ in } \mathbb{R}^r, s \mbox{ in } \mathbb{R}^+\} \\
&=& A^0 + \mathcal{L}^+(\hat z^0) ,
\end{eqnarray*}
where $A^0$ is as in Example~\ref{ex1}. This is the solution surface from
Example~\ref{ex1}, translated from $\hat y^0$ to $y^0$ and then
positively radiated in the $\hat z^0$ direction.
\end{example}

\begin{example}[(The transformation model)]\label{ex4}
The transformation model (see, e.g., \citep{fra1}) provides a
paradigm for exact ancillary conditioning. A typical continuous
transformation model for a variable $y=\theta z$ has parameter $\theta$
in a smooth transformation group $G$ that
operates on an $n$-dimensional sample space for $y$;
for illustration, we assume here that the group acts coordinate by coordinate.
The natural quantile function for the $i$th coordinate is $y_i=\theta
z_i$, where $z_i$ is a coordinate reference variable with a fixed
distribution; the linear regression model with known and unknown error
scaling are simple examples.
With observed data point $y^0$, let $\hat\theta^0$ be the maximum likelihood
value and $\hat z^0$ the corresponding reference value satisfying
$y^0=\hat\theta^0 \hat z^0$. The second order approximate ancillary is
then given as $\{\theta\hat z^0\}$, which is just the usual transformation
model orbit $G \hat z^0$. If the group does not apply separately to independent
coordinates, then the present quantile approach may not be immediately
applicable; this raises issues for the construction of the trajectories
and also
for the construction of default priors (see, e.g.,
\citep{ber}). Some discussion of this in connection with
curved parameters will be reported separately. A modification achieved
by adding structure to the transformation model is given by the
structural model \citep{fra1}. This takes the reference distribution
for $z$ as the primary probability space for the model
and examines what events on that space are identifiable from an
observed response; we do not address here this alternative modelling approach.
\end{example}

\begin{example}[(The inverted Cauchy)]
Consider a location-scale model centered at $\mu$ and scaled by $\sigma
$ with error given by the standard Cauchy; this gives the
statistical model
\[
f(y;\mu,\sigma)=\frac{1}{\curpi\sigma\{1+ (y-\mu)^2/\sigma^2\}}
\]
on the real line.
For the sampling version, this location-scale model is an example of
the transformation model discussed in the
preceding Example~\ref{ex4} and the long-accepted ancillary contour is the
half-plane (\ref{lsanc}).

McCullagh \citep{mcc2} uses linear fractional transformation results
that show that
the inversion $\tilde y=1/y$ takes the Cauchy ($\mu, \sigma$) model for $y$
into a Cauchy ($\tilde\mu, \tilde\sigma$) model for $\tilde y$, where
$\tilde\mu= {\mu}/(\mu^2+\sigma^2), \tilde\sigma= {\sigma}/(\mu
^2+\sigma^2)$.
He then notes that the usual location-scale ancillary for the derived
model does not map back to give the usual location-scale
ancillary on the initial space and would thus
typically give different inference results for the parameters; he indicates
``not that conditioning is a bad idea, but that the usual mathematical
formulation is
in some respects ad hoc and not completely satisfactory.''

We illustrate this for $n=2$ in Figure~\ref{an}.
For a data point in the upper-left portion of the plane in part (b)
for the inverted Cauchy, the observed ancillary contour is shown as a
shaded area; it is a half-plane subtended by $\mathcal{L} (1)$. When this
contour is mapped back to the initial plane in part (a), the contour
becomes three disconnected segments
with lightly shaded edges indicating the boundaries; in particular, the
line with
marks 1, 2, 3, 4, 5, 6 becomes three distinct curves again with
corresponding marks 1, 2, 3, 4, 5, 6, but
two points $(0,1), (1,0)$ on the line have no back images. Indeed, the
same type of singularity,
where a point with a zero coordinate cannot be mapped back, happens for any
sample size $n$.
Thus the proposed sample space is not one-to-one continuously
equivalent to the given sample space: points are left out and points
are created. And the quantile function used
on the proposed sample space for constructing the ancillary
does not exist on the given sample space: indeed, it is not defined at
points and is thus not continuous.

\begin{figure}

\includegraphics{248f02.eps}

\caption{(a) The location-scale Cauchy model for the inverted
$\tilde{y}_1 = 1/y_1$, $\tilde{y}_2 = 1/y_2$ has an ancillary contour
given by the shaded area in (b). When interpreted back for the
original $(y_1, y_2)$ the connected ancillary contour becomes three
unconnected regions, shown in (a). A line $\tilde{y}_2 = \tilde{y}_1
+1$ on the contour in (b) is mapped back to
three curved segments in (a)
and numbered points in sequence on the line are mapped back to the numbered
points on the unconnected ancillary contour.}
\label{an}
\end{figure}

The Cauchy inversion about $0$ could equally be about an arbitrary
point, say $a$,
on the real line and would lead to a corresponding ancillary. We would
thus have a wealth of competing ancillaries and a corresponding wealth
of inference procedures, and all
would have the same lack of one-to-one continuous equivalence to the
initial sample space.
While Fisher seems not
to have explicitly specified continuity as a needed ingredient for
typical ancillarity,
it also seems unlikely that he would have envisaged ancillarity without
continuity.
If continuity is included in the prescription for developing the
ancillary, then the proposed ancillary for the inverted Cauchy would
not arise.

Bayesian statistics involves full conditioning on the observed data and
familiar frequentist inference avoids, perhaps even evades,
conditioning. Ancillarity, however, represents
an intermediate or partial conditioning and, as such, offers a partial
bridging of the two extreme approaches to inference.
\end{example}

\section{An asymptotic statistic}\label{sec4}
For the Normal-on-the-circle example, the exact ancillary contour was
given as the observed
contour of the radial distance $r(y_1, y_2)$: the contour is described
implicitly. By contrast, the approximate ancillary was given
as the trajectory of a point $y(\hat x^0; t)$ under change of an index
or mathematical
parameter $t$: the contour is described explicitly. For the
general context, the first approach has serious difficulties, as found
even with nonlinear
regression, and these difficulties arise with an approximate statistic
taking an approximate
value; see Example~\ref{ex2}. Accordingly, we now turn to the second, the explicit
approach, and develop the needed notation and expansions.

Consider a smooth one-dimensional contour through
some point $y_{0}$.
To describe such a contour in the implicit manner requires
$n - 1$ complementary statistics. By contrast, for the explicit method,
we write $y=y(t)$, which maps a
scalar $t$ into the sample space $\mathbb{R}^n$. More generally, for a $p$-dimensional
contour, we have $y=y(t)$ in $\mathbb{R}^n$,
where $t$
has dimension $p$ and the mapping is again into $\mathbb{R}^n$.

For such a contour, we define the row array $
V(t)= (\mathrm{d}/{\mathrm{d}t'}) y(t)=\{v_1(t),\ldots,v_p(t)\}$
of tangent vectors, where the vector
$v_\alpha(t) =(\mathrm{d}/\mathrm{d}t_\alpha)y(t)$
gives the direction or gradient of $y(t)$
with respect to change in
a coordinate $t_\alpha$.
We are interested in such a contour near a particular point $y_0 =
y(t_0)$; for convenience, we often choose $y_0$ to be the observed data
point $y^0$
and the $t_0$ to be centered
so that $t_0=0$.
In particular,
the array $V=V(t_0)$ of tangent vectors at a particular data point
$y_0$ will
be of special interest. The vectors in $V$ generate a
tangent plane $\mathcal{L} (V)$ at the point $y_0$
and this plane provides a linear
approximation to the contour. Differential geometry gives
length
properties of such vectors as the first fundamental form:
\[
V'V =
\pmatrix{
v_1\cdot v_1&\cdots& v_1\cdot v_p \cr
\vdots& & \vdots\cr
v_p\cdot v_1 & \cdots& v_p\cdot v_p
}
=
\pmatrix{
v_1^\prime v_1 & \cdots& v_1^\prime v_p\cr
\vdots& & \vdots\cr
v_p^\prime v_1 & \cdots& v_p^\prime v_p
}
;
\]
this records the matrix of inner products for the vectors $V$
as inherited from the inner product on $\mathbb{R}^n$. A
change in the parameterization $\tilde t=t(t)$ of the contour will
give different tangent vectors $V$, the same tangent plane $\mathcal{L}(V)$
and a different, but corresponding, first fundamental form.

Now, consider the derivatives of the tangents $V(t)$ at $t_0$:
\[
\label{second}
W =\frac{\mathrm{d}}{\mathrm{d}t'}V(t) \bigg|_{t=t_0}=
\pmatrix{
w_{11} &\cdots& w_{1p}\cr
\vdots & & \vdots\cr
w_{p1} & \cdots& w_{pp}
}
,
\]
where
$w_{\alpha\alpha'}=(\partial^2/\partial t_\alpha\,\partial t_{\alpha
'})y(t)|_{t=t_0}$
is an acceleration or curvature vector relative
to coordinates $t_\alpha$ and $t_{\alpha^\prime}$ at $t_0$.
We regard the array $W$ as a $p\times p$ array
of vectors in $\mathbb{R}^n$. We
could have used tensor notation,
but the approach here has the advantage
that we can write
the second degree Taylor expansion of $y(t)$ at $t_0=0$ as
\begin{equation}
\label{tayloranc}
y(t)= y_0+Vt + t'Wt/2+ \cdots,
\end{equation}
which uses matrix multiplication for linearly combining the vectors in the
arrays $V$ and $W$.
Some important characteristics of the quadratic term in (\ref{tayloranc}) are
obtained by orthogonalizing the elements of $W$ to the tangent plane
$\mathcal{L} (V)$, to give residuals
\[
\tilde w_{\alpha\alpha'}= \{I-V(V'V)^{-1}V'\}w_{\alpha\alpha'} =
w_{\alpha\alpha'} -
Pw_{\alpha\alpha'};
\]
this uses the regression analysis projection matrix $P = V(V^\prime
V)^{-1}V^\prime$. The
full array $\tilde W$ of such vectors $\tilde w_{\alpha\alpha'}$
is then written
$
\tilde W = W - PW = W - VH,
$
where $H = (h_{\alpha\alpha'})$ is a $p \times p$ array of elements
$h_{\alpha\alpha'} =
(V^\prime V)^{-1}V^\prime w_{\alpha\alpha'}$; an element $h_{\alpha
\alpha'}$ is a $p
\times1$ vector, which records the regression coefficients of
$w_{\alpha\alpha'}$ on the vectors
$V$.

The array $\tilde W$ of such orthogonalized curvature vectors $\tilde
w$ is the second
fundamental form for the contour at the expansion point. Consider the
Taylor expansion
(\ref{tayloranc}) and substitute $W = \tilde W + VH$:
\begin{eqnarray*}
y(t) & = & y_0 + Vt + t^\prime(\tilde W + VH)t/2 + \cdots\\
& = & y_0 + V(t + t^\prime Ht/2) + t^\prime\tilde W t/2 + \cdots
,
\end{eqnarray*}
where we note that $t$ and $t^\prime$ are being applied to the $p
\times p$ arrays $H$
and $\tilde W$ by matrix multiplication, but the elements are $p \times
1$ vectors for $H$
and $n \times1$ vectors for $\tilde W$, and these are being combined
linearly. We can
then write $ y(t) = y_0 + V\tilde t + \tilde t \tilde W \tilde t^\prime
/2 + \cdots$
and thus have the contour expressed in terms of orthogonal curvature
vectors $\tilde w$
with the reparameterization $\tilde t = t + t^\prime H t/2 + \cdots$.
When we use this in the asymptotic setting, we will have standardized
coordinates and the
reparameterization will take the form $\tilde t = t + t^\prime H
t/2n^{1/2} + \cdots$.

\section{Verifying second order ancillarity}\label{sec5}
We have used the Normal-on-the-circle example to illustrate the
proposed second order
ancillary contour $\{y(\hat x^0; t)\}$. Now, generally, let $f(y; \theta
)$ be a
statistical model with regularity and asymptotic properties as the data
dimension $ n$ increases: we
assume that the vector quantile $y(x; \theta)$ has independent scalar
coordinates and is
smooth in both the reference variable $x$ and the parameter $\theta$;
more general
conditions will be considered subsequently. For the verification, we
use a Taylor
expansion of the quantile function in terms of both $x$ and $\theta$,
and work from theory
developed in \citep{cak} and \citep{and}. The first steps involve the
re-expression of individual coordinates of $y$, $x$, and $\theta$, and
show that the
proposed contours establish a partition on the sample space; the
subsequent steps
establish the ancillarity of the contours.
\begin{longlist}[(1a)]
\item[(1a)] \textit{Standardizing the coordinates.}
Consider the statistical model in moderate deviations about $(y^0,
\hat\theta^0)$ to order $\mathrm{O}(n^{-1})$. For this, we work with coordinate
departures in
units scaled by $n^{-1/2}$. Thus, for the $i$th coordinate, we write
$y_i = \hat y_i^0 +
\tilde y_i/n^{1/2}$, $x_i = \hat x_i^0 + \tilde x_i/n^{1/2}$ and $\theta
_\alpha=
\hat\theta^0_\alpha+ \tilde\theta_\alpha/n^{1/2}$; and for
a modified $i$th quantile coordinate $\tilde y_i = \tilde y_i(\tilde
x_i, \hat\theta)$,
we Taylor expand to the second order, omit the subscripts and tildes
for temporary
clarity, and obtain $
y = x + V\theta+ (ax^{2} + 2xB\theta+ \theta'W\theta)/2n^{1/2}$,
where $V$ is the $1 \times p$ gradient of $y$ with respect to $\theta$,
$B$ is the $1\times p$ cross Hessian with respect to $x$ and $\theta$,
$W$ is the $p \times p$ Hessian with respect to $\theta$
and vector--matrix multiplication is used for combining $\theta$ with
the arrays.

\item[(1b)] \textit{Re-expressing coordinates for a nicer expansion.}
We next re-express an $x$ coordinate, writing $\tilde x = x
+ ax^2/2n^{1/2}$,
and then again omit the tildes to obtain the simpler expansion
\begin{equation}
\label{scalar}
y = x + V\theta+ (2xB\theta+ \theta'W\theta)/2n^{1/2} + \cdots,
\end{equation}
to order $\mathrm{O}(n^{-1})$ for the modified $y$, $x$ and $\theta$, now in
bounded regions about
0.

\item[(1c)] \textit{Full response vector expansion.}
For the vector response $y = (y_1,\ldots, y_n)$ in quantile form, we can
compound the
preceding coordinate expansions and write
$y = x + V\theta+ (2x : B\theta+ \theta'W\theta)/ 2n^{1/2} +
\cdots,$
where $y$ and $x$ are now vectors in $\mathbb{R}^{n}$, $V = (v_1,\ldots, v_p) =
(v_\alpha)$ and $B =
(b_1,\ldots, b_p) = (b_\alpha)$ are $1 \times p$ arrays of vectors in $\mathbb{R}^{n}$,
$W = (w_{\alpha\alpha^\prime})$ is a $p\times p$ array of vectors in $\mathbb{R}^{n}$
and $x \mbox{:} B$ is a $1 \times p $ array of vectors $x \mbox{:} b$, where
the $i$th element of
the vector $x \mbox{:} b$ is the product $x_ib_i$ of the $i$th elements of
the vectors $x$
and $b$.

\item[(1d)] \textit{Eliminate the cross Hessian: scalar parameter case.}
The form of a Taylor series depends heavily on how the function and the
component
variables are expressed. For a particular coordinate of (\ref{scalar})
in (1b),
if we re-express the coordinate $ y=\tilde y + c\tilde y^2/2n^{1/2}$ in
terms of a
modified $\tilde y$, substitute it in (\ref{scalar}) and then, for
notational ease, omit the tildes, we obtain
$y+c(x+v\theta)^2/2n^{1/2}=x+v\theta+(2xb\theta+\theta^2w)/2n^{1/2}.$
To simplify this, we take the $x^2$ term over to the right-hand side
and combine it with $x$ to give a re-expressed $x$, take the $\theta x$
term over to the right-hand
side and choose $c$ so that $cv=b$ and, finally, combine the $\theta^2$
terms giving a new $w$. We then obtain
${y(x; \theta) =x+v\theta+\theta^2w/2n^{1/2}}$
with the cross Hessian removed; for this, if $v = 0$, we ignore the
coordinate as being ineffective for $\theta$.
For the full response accordingly, we then have
$
y(x; \theta) = x + v\theta+ w\theta^2/2n^{1/2}+ \cdots$
to the second order in terms of re-expressed coordinates $x$ and $y$. The
trajectory of a point $x$ is
$A(x) = \{y(x; t)\} = \{x + vt + wt^2/2n^{1/2}+\cdots\}$
to the second order as $t$ varies.

\item[(1e)] \textit{Scalar case: trajectories form a partition.}
In the standardized coordinates, the initial data point is $y^0 = 0$
with corresponding
maximum likelihood value $\hat\theta^0 = 0$; the corresponding
trajectory is
$A(0) = \{vt + wt^2/2n^{1/2} + \cdots\}.$
For a general reference value $x$, but with $\hat\theta(x) = 0$, the
trajectory is
$A(x) = \{x + vt + wt^2/2n^{1/2} +\cdots\} \nonumber= x + A(0)$.
The sets $\{A(x)\}$ with $\hat\theta(x) = 0$ are all translates of
$A(0)$ and thus
form a partition.

Consider an initial point $x_0$ with maximum likelihood value $\hat
\theta(x_0) = 0$ and let $y_1= x_0 + vt_1 + wt_1^2/2n^{1/2} + \cdots$
be a point in
the set $A(x_0) = x_0 + A(0)$.
We calculate the trajectory $A(y_1)$ of $y_1$ and show that it lies on
$A(x_0)$; the
partition property then follows and the related Jacobian effect is constant.
From the quantile function $y = x + v\theta+ w\theta^2/2n^{1/2}$, we
see that the $y$
distribution is a $\theta$-based translation of the reference
distribution described by
$x$. Thus the likelihood at $y_1$ is $l(y_1 - v\theta- w\theta
^2/2n^{1/2})$, in terms of
the log density $l(x)$ near $x_0$. It follows that $y_1 = x_0 + vt_1 +
wt_1^2/2n^{1/2}$
has maximum likelihood value $\hat\theta(y_1) = t_1$.

Now, for the trajectory about $y_1$, we calculate derivatives
\[
\frac{\mathrm{d}y}{\mathrm{d}\theta} =v+w\theta/n^{1/2}, \qquad \frac{\mathrm{d}^2y}{\mathrm{d}\theta^2}=w/n^{1/2},
\]
which, at the point $y_1=vt_1 + wt^2_1/2n^{1/2}$ with $\theta=\hat\theta
(y_1)$, gives
\[
V(y_1) =v+wt_1/n^{1/2}, \qquad  W(y_1)=w/n^{1/2},
\]
to order $\mathrm{O}(n^{-1})$.
We thus
obtain the trajectory of the point $y_1$:
\begin{eqnarray*}
A(y_1) &=& \{x_0 + vt_1+
wt_1^2/2n^{1/2}+(v+wt_1/n^{1/2})t+wt^2/2n^{1/2}\} \\
&=& \{x_0 + vT +wT^2/2n^{1/2}\}
\end{eqnarray*}
under variation in $t$.
However, with $T=t_1+ t$, we have just an arbitrary point on the
initial trajectory. Thus the mapping $y \to A(y)$ is well defined and
the trajectories generate a partition, to second order in moderate
derivations in
$\mathbb{R}^n$. In the standardized coordinates, the Jacobian effect is constant.

\item[(1f)] \textit{Vector case: trajectories form a partition}.
For the vector parameter case, we again use
standardized coordinates and
choose a parameterization that
gives orthogonal curvature vectors $w$
at the observed data point $y^0$. We then
examine scalar parameter change on some line through $\hat\theta(y^0)$.
For this, the results above give a trajectory and
any point on it reproduces the trajectory under that scalar parameter.
Orthogonality ensures that the vector maximum likelihood value is on
the same line
just considered. These trajectories are,
of course, part of the surface defined by $\{Vt+t'Wt/2n^{1/2}\}$.
We then use the partition property of the individual trajectories
as these apply perpendicular to the surface; the surfaces are thus part
of a partition. We can then write the trajectory of a point $x$ as a set
\begin{equation}
\label{taylor3}
A(x)=\{ x + Vt + t'Wt/ 2n^{1/2} + \cdots\dvt t\}=x + A(0)
\end{equation}
in a partition to the second order in moderate deviations.

\item[(2a)] \textit{Observed information standardization}.
With moderate regularity, and
following \citep{fra5} and \citep{fra10}, we have a limiting Normal
distribution conditionally on
$ y^0+\mathcal{L}(V)$. We then rescale the parameter
at $\hat\theta^0$ to give identity observed information and thus an identity
variance matrix for the Normal distribution to second order.
We also have a limiting
Normal distribution conditionally on $ y^0+\mathcal{L}(V,W)$; for this,
we linearly modify the vectors in $W$ by rescaling and regressing on
$\mathcal{L}(V)$
to give distributional orthogonality to $\hat\theta$ and identity
conditional variance
matrix to second order.

\item[(2b)] \textit{The trajectories are ancillary: first derivative parameter change}.
We saw in the preceding section that key local properties of a statistical
model were summarized by the tangent vectors $V$ and the curvature vectors
$W$, and that the latter can, to advantage, be taken to be orthogonal
to the
tangent vectors. These vectors give local coordinates for the model and
can be replaced by
an appropriate subset if linear dependencies are present.

First, consider the conditional model given the
directions corresponding to the span $y^0+\mathcal{L}\{V, W\}$. From the ancillary
expansion (\ref{taylor3}), we have that change of $\theta$ to the
second order
moves points within the linear space $y^0+\mathcal{L}\{V, W\}$; accordingly,
this conditioning is ancillary.
Then, consider the further conditioning to an alleged ancillary
contour, as described by (\ref{taylor3}). Also, let $y_0$
be a typical point having
$\hat\theta(y_0)=\hat\theta^0$ as the corresponding maximum likelihood value;
$y_0$ is thus on the observed maximum likelihood contour.

Now, consider a rotationally symmetric Normal distribution on the $(x,
y)$ plane
with mean $\theta$ on the $x$ axis and let $a=y+cx^2/2$ be linear in $y$
with a quadratic adjustment with respect to~$x$. Then $a=a(x,y)$ is
first-derivative
ancillary at $\theta=0$.
For this, we assume,
without loss of generality, that the standard deviations
are unity. The marginal density for $a$ is then
\[
f(a;\theta)=\int_{-\infty}^{\infty}\phi(x-\theta)\phi(a-cx^2/2)\,\mathrm{d}x ,
\]
which is symmetric in $\theta$; thus
$(\mathrm{d} / \mathrm{d}\theta)f(a;\theta) |_{\theta=0}=0,$
showing
that the distribution of $a$ is first-derivative ancillary
at $\theta=0$ or, more intuitively, that the amount of probability on
a contour of $a$ is first-derivative free of $\theta$ at $\theta=0$.
Of course, for this, the $y$-spacing between contours of $a$ is constant.

Now, more generally, consider an asymptotic distribution for $(x,y)$
that is first order
rotationally symmetric Normal with mean $\theta$ on the $y=0$ plane;
this allows
$\mathrm{O}(n^{-1/2})$ cubic contributions. Also, consider an $s$-dimensional variable
$a=y+Q(x)/2n^{1/2}$ which is a quadratic adjustment of $y$. The
preceding argument
extends to show that $a(y)$ is first-derivative ancillary: the two $\mathrm{O}(n^{-1/2})$
effects are zero and the combination is of the next order.

\item[(2c)] \textit{Trajectories are ancillary: parameter change in moderate deviations.}
Now, consider a statistical model $f(y;\theta)$ with data point $y^0$ and
assume regularity, asymptotics and smoothness of the quantile
functions. We examine the parameter trajectory $\{y(\hat x^0; t)\}$
in moderate deviations under change in $t$. From the preceding
paragraph, we then have
first-derivative
ancillarity at $\theta=\hat\theta=0$. But this holds for each expansion
in moderate
deviations
and we thus have ancillarity in moderate deviations. The key here
has been to use the expansion form about the point that has $\hat\theta$
equal to the parameter value being examined.
\end{longlist}
\section{Discussion}\label{sec6}

\begin{longlist}[(iii)]
\item[(i)] \textit{On ancillarity.}
The \hyperref[intro]{Introduction} gave a brief background on ancillary statistics
and noted that an ancillary is typically viewed as a statistic with a
parameter-free distribution; for some recent discussion, see \citep{fra4}.
Much of the literature is concerned with difficulties that can arise
using this third
Fisher concept, third after sufficiency and likelihood: that maximizing
power given size
typically means not conditioning on an ancillary; that shorter
on-average confidence intervals typically mean ignoring ancillary
conditioning; that techniques that are conditional on an ancillary are
often inadmissible; and more. Some of the difficulty may hinge on whether
there is merit in the various optimality criteria themselves. However,
little in the literature seems
focused on the continued evolution and development of this Fisher
concept, that is, on what
modifications or evolution can continue the exploration initiated in
Fisher's original papers (see \citep{fis1,fis2,fis3}).

\item[(ii)] \textit{On simulations for the conditional model.}
The second order ancillary in moderate deviations has contours that
form a partition, as shown in the preceding section. In the modified or
re-expressed coordinates, the contours
are in a location relationship and, correspondingly, the Jacobian
effect needed for the
conditional distribution is constant. However, in the original coordinates,
the Jacobian effect would typically not be constant and its effect
would be needed
for simulations. If the parameter is scalar, then the effect is
available to the second order
through the divergence function of a vector field; for some discussion
and examples, see \citep{fra2}. For a vector parameter, generalizations
can be implemented, but we do not
pursue these here.

\item[(iii)] \textit{Marginal or conditional.}
When sampling from a scalar distribution
having variable $y$ and moderate regularity,
the familiar central limit theorem gives a limiting Normal distribution
for the sample average $\bar y$ or
sample sum $\sum y_i$.
From a geometric view, we have probability in
$n$-space and contours determined by
$\bar y$, contours that are planes perpendicular to the $1$-vector.
If we then collect the probability on a contour, plus or
minus a differential, and deposit it, say, on the intersection
of the contour with the span
$\mathcal{L} (1)$ of the $1$-vector, then we obtain a limiting Normal
distribution on $\mathcal{L} (1)$, using
$\bar y$ or $\sum y_i$ for location on that line.

A far less familiar Normal limit result applies in the same general
context, but with a totally different geometric decomposition. Consider lines
parallel to the $1$-vector, the affine cosets of $\mathcal{L} (1)$. On these
lines, plus or minus a differential, we then obtain a limiting Normal
distribution for location say $\bar y$ or $\sum y_i$.
In many ways, this conditional, rather than marginal, analysis
is much stronger and more useful.
The geometry, however, is different,
with planes perpendicular to $\mathcal{L} (1)$ being replaced
by points on lines parallel to $\mathcal{L} (1)$.

This generalizes giving a limiting conditional Normal distribution
on almost arbitrary smooth contours in a partition and it has wide application
in recent likelihood inference theory. It also provides third order accuracy
rather than the first order accuracy associated with the usual geometry.
In a simple sense, planes are replaced by lines or
by generalized contours and much stronger, though
less familiar, results are obtained.
For some background based on Taylor expansions of log-statistical
models, see \citep{cak,che} and \citep{and}.
\end{longlist}

\section*{Acknowledgements}

This research was supported by the Natural Sciences and Engineering
Research Council of
Canada. The authors wish to express deep appreciation to the referee
for very incisive comments. We also offer special thanks to Kexin Ji
for many contributions and support with the manuscript and the diagrams.

redpajama_set_name: RedPajamaArXiv