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

Jets from young stars and their associated Herbig-Haro objects (HH)
are now the most studied and best understood of all astrophysical jets.
HH jets are collimated bipolar ejections from low, intermediate and
sometimes high mass protostars
or young stars (see, e.g., the review of Frank et al. 2014). The interaction between HH jets and the surrounding environment generates a main working surface known as the ``head''. This structure (formed by a bow shock/jet shock pair) is a clear sign of the ``turning on'' of the jet, which tipically occurred at times of $\sim 10^3\to 10^4$~yr ago in observed HH jets.

This type of jets also have a structure of emitting knots (some of them resembling the
jet head, and others forming aligned chains of more compact emission peaks)
in the region between the outflow source and the jet head
(see, e.g., Reipurth \& Heathcote 1990 and the review of Reipurth \&
Bally 2001). Several possibilities
of how to model these knots have been explored (see, e.g., Raga \& Kofman, 1992, and Micono et al. 1998),
but presently the favoured explanation is that they are the result of a time-variability
in the ejection, which leads to the formation of ``internal working surfaces''
within the jet beam (see, e.g., Raga et al 1990). Some of the ``classical HH objects''
of the catalogue of \citet{H74} are associated with either jet heads or knots.

A problem that has received little attention is the effect
of a ``slow turning on'' of the outflow on the jet head. As far as we
are aware, the only simulations studying the dynamics of the head of
a ``slow turning on'' jet (as opposed to a jet which is instantaneously
``switched on'' at full velocity) are the ones of \citet{LETAL02}.
These authors focus their study on the survival of H$_2$
molecules in the resulting, accelerating
jet head,  presenting 1D hydrodynamic simulations including an H$_2$
formation/destruction chemical network.

Of course, there are a number of numerical studies of the ejection of
MHD jets from star+accretion disk systems (e.g., Ouyed et al. 2003;
Hiromitsu \& Shibata 2005; Zanni et al. 2007; Ahmed \& Shibata 2008;
Romanova et al. 2018), in which a sudden switch on of the jet is
definitely not imposed (the ``switching on'' of the outflow being
a result of the simulated accretion+ejection flow). These models
produce a highly superalfv\'enic leading head that evolves consistently
with the computed time-dependence of the ejection.

In the present paper, we revisit the ``slow turn on jet'' problem. We study the
dynamics of the head of a cylindrical jet with a linear ramp of increasing
ejection velocity as a function of ejection time. We combine this ejection velocity
with two forms for the ejection density variability:
\begin{itemize}
  \item a constant ejection density,
  \item a constant mass loss rate (so that the ejection density is proportional
  to the inverse of the ejection velocity).
\end{itemize}
We obtain analytical models of this flow by assuming that the material in the
jet beam is free streaming (before reaching the jet head), and that the motion
of the working surface of the jet head can be described:
\begin{itemize}
\item by using a ram-pressure balance condition,
\item by assuming that it coincides with the motion of the center of
  mass of the material that has entered the working surface.
\end{itemize}
The first of these assumptions (see, e.g., Raga \& Cant\'o 1998) is correct for
a ``massless'' working surface which instantaneously ejects most of the material sideways (into a jet cocoon), while the latter assumption
(see Cant\'o et al. 2000) is appropriate for the ``mass conserving'' case in which
the shocked material mostly stays within the working surface.

We also carry out axisymmetric numerical simulations of the flow, and compare
the results with the ram-pressure balance and center of mass analytic models.
This is the first time that a comparison between numerical simulations
and the two analytic approximations (described above) has been attempted.

The paper is organized as follows. In \S2 we present the center of mass and the
ram-pressure balance solutions for a uniformly accelerated jet evolving in a
homogeneous interstellar medium. In \S3 we present the axisymmetric numerical
simulations, and a comparison of the results with the analytic models is done
in \S4. A discussion of the results is held in \S5, and concluding remarks are
given in \S6.

\section{Analytic model}
\label{sec:a-m}

We consider a jet with an ejection velocity of the form
\begin{equation}
    u_0(\tau) = \left\{\begin{array}{lr}
         0 \hspace{0.1cm}, & \tau<0  \\
         a\tau \hspace{0.1cm} , & \tau \ge 0
    \end{array} \right.
    \label{ejec_vel}
\end{equation}
where $a$ is a constant acceleration and $\tau$ the ejection time.
We assume a cylindrical flow, with a jet cross section $\sigma$ which
is circular and constant. We also assume that the jet moves into an
environment of uniform density $\rho_a$.

We consider two different forms for the ejection density:
\begin{enumerate}
\item a constant mass loss rate per unit area ${\dot m}$,
  such that the ejection density is given by
    \begin{equation}
        \rho_0(\tau) = \frac{\dot{m}}{u_0(\tau)},
        \label{eq:rhotau}
    \end{equation}
    \item a time-independent ejection density,
    \begin{equation}
        \rho_0=\mathrm{const.}
        \label{eq:rhocte}
    \end{equation}
 \end{enumerate}

We assume that the flow in the jet beam is free-streaming, satisfying
the condition:
\begin{equation}
    u(x,t)=\frac{x}{t-\tau}=u_0(\tau)\,,
    \label{position}
\end{equation}
where $u(x,t)$ is the velocity as a function of distance $x$ from
the source at an evolutionary time $t$ and $u_0(\tau)$ is the ejection velocity
(at the ejection time $\tau<t$). For an ejection time $\tau<0$, $x=0$,
i.e., no material is ejected. Also, for a cylindrical, free-streaming flow,
the density is given by:
\begin{equation}
    \rho(x,t) = \frac{\rho_0\,u_0}{u_0 - (t - \tau){\dot u}_0},
\label{density}
\end{equation}
where $\rho_0(\tau)$ is the ejection density (see equation \ref{eq:rhotau} or \ref{eq:rhocte})
and ${\dot u}_0=du_0/d\tau$ is the derivative of the ejection velocity
with respect to the ejection time (for a derivation of this equation,
see \citeauthor{RK92}, 1992). The relationship between the evolutionary time $t$ and the ejection time $\tau$ can be obtained from equations (\ref{ejec_vel}) and (\ref{position}). Later, in section 2.1, we explicitly show this relationship.

When the flow starts (at $t=0$), a working surface is formed at $x=0$. This ``jet
head'' then travels away from the outflow source for increasing times. In order
to describe the time evolution of this working surface, we use two analytic
approximations:
\begin{enumerate}
\item \textit{center of mass}- we assume that the position of the working surface
  coincides with the center of mass of the ejected, free-streaming material,
\item \textit{ram-pressure balance}- we assume that the motion of the working
  surface is determined by the jet/environment ram-pressure balance condition.
\end{enumerate}
The first of these approximations is appropriate for a working surface in which
the shocked gas mostly remains within the jet head, and the latter approximation
is valid for a working surface that ejects most of the matter sideways (into
a jet cocoon).

We therefore develop four analytic models, with the two ejection densities
and the two analytic approximations described above. The models all share
the linearly accelerating ejection velocity given by equation (\ref{ejec_vel}).

 \subsection{Center of mass equation of motion}

 We first consider a ``mass conserving'' working surface. For this case,
 we assume that the material going through the bow and jet shocks
 mostly stays within the working surface, with only a small fraction
 of this material being ejected sideways into the jet cocoon.
 The working surface position should then coincide with the center of mass
 of the free-streaming fluid parcels that have piled up within it. This
 center of mass position is given by:
 \begin{equation}
     x_{cm} = \frac{\int xdm}{\int dm},
     \label{center_mass}
 \end{equation}
with the differential element of mass given by:
    \begin{equation}
        dm = \sigma\, \rho_0(\tau)\, u_0(\tau)\, d\tau+\sigma\, \rho_a(x)\, dx.
        \label{dm}
    \end{equation}
    In this equation, the first term on the right-hand-side corresponds
    to the ejected material,
    and the second term to the swept-up, stationary environment entering the
    working surface. In
    equation (\ref{dm}), $\sigma$ is the jet cross section, $\rho_0(\tau)$ is
the ejection density, $u_0(\tau)$ the ejection velocity and $\rho_a(x)$
the (possibly position-dependent) environmental density.

Using equations~(\ref{position}-\ref{dm}) and considering that the ejection
time $\tau'$ is integrated from 0 to $\tau$, we obtain:
\begin{multline}
     x_{cm} \left[\int_0^\tau \rho_0(\tau')\,u_0(\tau')\,d\tau' + \int_0^{x_{cm}} \rho_a\,dx \right] = \\ \int_0^\tau x_j\,\rho_0(\tau')\,u_0(\tau)\,d\tau' + \int_0^{x_{cm}}x\,\rho_a\,dx,
     \label{center_mass_2}
\end{multline}
where we have set $\rho_a=const.$, and
$x_j$ is the position that the fluid parcels would have if they were
still free-streaming:
\begin{equation}
  x_j = (t-\tau')\,u_0(\tau')\,,
  \label{xjttau}
\end{equation}
see equation (\ref{position}).

Also, from equations (\ref{ejec_vel}) and (\ref{position}) we find:
\begin{equation}
        t=\frac{x_{cm}}{a\tau}+\tau\,.
        \label{timet}
\end{equation}
Equation (\ref{timet}) can be inverted to obtain
\begin{equation}
        \tau = \frac{1}{2}\left(t + \sqrt{t^2 - \frac{4x_{cm}}{a}} \right).
        \label{timetau}
\end{equation}

Combining equations (\ref{center_mass_2}-\ref{timet}), we obtain:
\begin{multline}
  x_{cm} \left[\int_0^\tau \rho_0(\tau')\,a\tau'\,d\tau' + \frac{\rho_a\,x_{cm}}{2}
    \right] = \\
  \int_0^\tau (t-\tau')a^2\,\tau'^2\,\rho_0(\tau')\,d\tau'\,.
     \label{center_mass3}
\end{multline}

In order to carry out the remaining integrals, we have to specify
$\rho_0(\tau')$. We consider two forms for
the ejection density (see equation \ref{eq:rhotau}):

\begin{enumerate}[a)]
\item { \it Constant mass loss rate }

equation (\ref{center_mass3}) takes the form:
      \begin{multline}
        x_{cm}\left(\dot{m}\tau + \frac{\rho_a\,x_{cm}}{2} \right) = \dot{m}a\tau^2\left(\frac{t}{2} - \frac{\tau}{3} \right).
        \label{ctemassloss}
    \end{multline}
   Using equations (\ref{timet}) and (\ref{ctemassloss}), we obtain
    \begin{equation}
        x_{cm}^2 + \frac{\dot{m}\tau}{\rho_a}\,x_{cm} - \frac{\dot{m}a\tau^3}{3\rho_a} = 0,
        \label{solution_mlc}
    \end{equation}
this is a quadratic equation with one positive solution, given by
    \begin{equation}
        x_{cm}(\tau) = \frac{\dot{m} \tau}{2\rho_a} \left[-1 + \sqrt{1+\frac{4a\rho_a \tau}{3\dot{m}}} \right].
        \label{mdotcte_solution}
    \end{equation}

   \begin{figure*}[ht]
   \centering
   \includegraphics[scale=0.55]{bowshock_111.eps}
   \includegraphics[scale=0.55]{bowshock_222.eps}

   \caption{\footnotesize{Jet  head working surface in a reference frame at rest with respect to the outflow source (left) and in a reference frame moving along with the working surface (right). The bow shock is in red and the jet shock is in blue.}}
   \label{bowshock}%
   \end{figure*}

    In the low density limit $\rho_a \rightarrow 0$, equation (\ref{solution_mlc}) takes a particularly simple form that leads to finding the solution
    \begin{equation}
        x_{cm}=\frac{a\tau^2}{3}.
    \end{equation}
This solution can also be obtained by carrying out a first-order
    Taylor series expansion of equation (\ref{mdotcte_solution}).

   By substituting equation (\ref{timetau}) into (\ref{ctemassloss}), we can also obtain $x_{cm}$ as an explicit function of $t$, giving
    \begin{equation}
        x_{cm} = \frac{8}{9}x_c \left[\left(1+\frac{3}{4} \frac{t}{t_c} \right)^{3/2} - \left(1+ \frac{9}{8} \frac{t}{t_c} \right) \right],
        \label{sol_xcm_t}
    \end{equation}
   with, $
        t_c \equiv \dot{m}/a\rho_{a}$ and $ x_c \equiv \dot{m}^2/a\rho_a^2$.

        In the $t \ll t_c$ limit, a second order Taylor series
          expansion of equation (\ref{sol_xcm_t}) gives:
    \begin{equation}
        x_{cm} \approx \frac{3}{16} at^2.
    \end{equation}

    \vspace{0.5cm}
  \item { \it Constant ejection density }

With a constant $\rho_0$, equation (\ref{center_mass3}) takes the form:
    \begin{multline}
        x_{cm}\left(\frac{\rho_0\,a\tau^2}{2} + \frac{\rho_a\,x_{cm}}{2} \right) = \\ \rho_0\,a^2\tau^3 \left(\frac{t}{3} - \frac{\tau}{4} \right)\,,
        \label{eqmotion_rhoc}
    \end{multline}
where $t$ and $\tau$ are related through equation (\ref{timet}).

Substituting $t$ as a function of $\tau$ we obtain:
    \begin{equation}
      x_{cm}^2 + \frac{\rho_0\,a\tau^2}{3\rho_a}x_{cm} -\frac{\rho_0\,a^2\tau^4}{6\rho_a}=0,
      \label{xxxcm}
    \end{equation}
    which has a positive solution
    \begin{equation}
        x_{cm}(\tau) = \frac{\rho_0\,a\tau^2}{6\rho_a}\left[-1 + \sqrt{1+\frac{6\rho_a}{\rho_0}} \right].
        \label{solution_rhoc}
    \end{equation}
    In the $\rho_a \rightarrow 0$ limit, equation (\ref{xxxcm}) (or,
    alternatively, equation \ref{solution_rhoc}) takes the form
    \begin{equation}
        x_{cm}=\frac{a\tau^2}{2},
    \end{equation}
    with a quadratic dependency on the ejection time.

    We can also find a solution to (\ref{eqmotion_rhoc}) as a function of
    evolutionary time $t$ (using equation \ref{timetau}), giving:
    \begin{equation}
        x_{cm} = \frac{a}{9}\beta_0 \left[\frac{\beta_0 (\beta_0^2 -18) + (\beta_0^2 + 6)^{3/2}}{(\beta_0^2 - 2)^2} \right] t^2,
        \label{solution_timet}
    \end{equation}
    where,
    \begin{equation}
        \beta_0 = \sqrt{\frac{\rho_0}{\rho_a}},
        \label{beta0}
    \end{equation}
    and therefore the position $x_{cm}$ has a quadratic dependency on
    the evolutionary time $t$.

    Equation (\ref{solution_timet}) takes the form of a constant acceleration condition, with the acceleration given by
    \begin{equation}
        g = \frac{2a}{9}\beta_0 \left[\frac{\beta_0 (\beta_0^2 -18) + (\beta_0^2 + 6)^{3/2}}{(\beta_0^2 - 2)^2} \right],
    \end{equation}
    which in the $\rho_a \rightarrow 0$ low density limit becomes,
    \begin{equation}
        g \approx \left(\frac{2}{3} \right)^{1/2} a\, \beta_0.
    \end{equation}

\end{enumerate}

\subsection{Ram-pressure balance equation of motion}

Figure \ref{bowshock} shows a schematic diagram of a jet head propagating
through a stationary environment. This leading working surface has
a structure consisting of two shocks: an upstream shock known as the
jet shock (or Mach disk), which is slowing down the jet material, and a downstream
bow shock which is accelerating the surrounding material. The two-shock
working surface structure travels away from the outflow source at a velocity
$v_{ws}$.

The right panel of Figure \ref{bowshock} shows the situation seen in
a reference frame at motion with the jet head. For a hypersonic flow, the
two working surface shocks are strong, so that the post-shock gas pressures
are given by:
\begin{equation}
    P_{bs}=\frac{2}{\gamma+1}\rho_a v_{ws}^2,
\end{equation}
for the bow shock, and
\begin{equation}
    P_{js}=\frac{2}{\gamma+1}\rho_j\,(v_j - v_{ws})^2,
\end{equation}
for the jet shock (assuming that the jet and the environment have
the same specific heat ratio $\gamma$). In this equation,
$v_j$ and $\rho_j$ are the velocity (with respect to the outflow
  source) and the density of the material that is presently entering
  the jet shock.

If the working surface is moving with a constant velocity then the condition
\begin{equation}
     P_{bs} = P_{js}\; \Rightarrow\;  \rho_a v_{ws}^2 = \rho_j\,(v_j - v_{ws})^2,
     \label{rampressure_balance}
\end{equation}
is satisfied. This is called the
\textit{ram-pressure balance} condition. This condition is also valid for a
working surface moving with a variable velocity as long as the inertia of
the material between the bow shock and the jet shock is negligible, which is
the case if most of the material is ejected sideways from the working surface
into a jet cocoon.

From equation (\ref{rampressure_balance}) we find that
\begin{equation}
  v_{ws} = \frac{\beta v_j}{1+\beta}\hspace{0.2cm},\;\hspace{0.5cm}
  {\rm with} \hspace{0.5cm}\beta \equiv \sqrt{\frac{\rho_j}{\rho_a}}.
\end{equation}

Clearly, the shock velocity associated with the bow shock is $v_{bs}=v_{ws}$ and
the shock velocity of the jet shock is:
\begin{equation}
    v_{js}=\frac{v_j}{1+\beta}.
\end{equation}

We now consider the equation of motion $v_{ws}=dx_{ws}/dt$ for the working
surface. Setting $v_j=u_0(\tau)$ (i.e., the velocity of the material entering
the working surface is equal to the ejection velocity at the corresponding
ejection time $\tau$), we then have:
\begin{equation}
    \left(1+\sqrt{\frac{\rho_a}{\rho_j}} \right)\frac{dx_{ws}}{dt} = u_0(\tau).
    \label{motion_rampress}
\end{equation}
The ejection time $\tau$ of the material entering the working surface is obtained
from the free-streaming relation
\begin{equation*}
    x_{ws}=(t-\tau)u_0(\tau),
\end{equation*}
where $t$ is the evolutionary time (see equation  \ref{position}).
Solving for $t$ and differentiating
with respect to $\tau$ we obtain
\begin{equation}
    \frac{dt}{d\tau} = 1 + \frac{1}{u_0}\frac{dx_{ws}}{d\tau} - \frac{x_{ws}}{u_0^2}\frac{du_0}{d\tau},
    \label{diff_times}
\end{equation}
Combining equations (\ref{motion_rampress}) and (\ref{diff_times}), we obtain:
\begin{equation}
    \frac{dx_{ws}}{d\tau}\sqrt{\frac{\rho_a}{\rho_j}} = u_0 -x_{ws}\frac{d\ln{u_0}}{d\tau}.
    \label{diffeqws}
\end{equation}
The density of the jet material entering the working surface is
$\rho_j=\rho(x_{ws},t)$ (see equation \ref{density}), which for our chosen
ejection velocity variability (see equation \ref{ejec_vel}) becomes:
\begin{equation}
    \rho_j= \frac{\rho_0}{1-\frac{x_{ws}}{a\tau^2}}.
    \label{dens_tau_x}
\end{equation}
Finally, from equations (\ref{diffeqws}) and (\ref{dens_tau_x}) we obtain:
\begin{equation}
  \frac{dx_{ws}}{d\tau}= \sqrt{\frac{\rho_0}{\rho_a}}\,a\tau\,
  \sqrt{1-\frac{x_{ws}}{a\tau^2}}.
    \label{diffeqwsrhoc}
\end{equation}

In order to proceed it is now necessary to specify the form of $\rho_0(\tau)$.
We therefore consider the two cases of constant mass loss rate and constant
ejection density.

\begin{enumerate}[a)]
\item { \it Constant mass loss rate }

   for $\dot{m}=\mathrm{const.}$ equation (\ref{diffeqwsrhoc}) takes the form:
    \begin{equation}
        \frac{dx_{ws}}{d\tau} = \sqrt{\frac{\dot{m}}{\rho_a} \left(a\tau - \frac{x_{ws}}{\tau}\right)}.
        \label{ctemasslosseqqdiff}
    \end{equation}
    It is convenient to use the dimensionless variables:
    \begin{equation}
        \eta = \frac{a\rho_a^2}{\dot{m}^2}x_{ws} \;,\;\;\;y=\frac{a\rho_a}{\dot{m}}\tau.
        \label{dimensionless_var}
    \end{equation}
    In terms of these variables, equation (\ref{ctemasslosseqqdiff}) is
    \begin{equation}
        \frac{d\eta}{dy} = \sqrt{y - \frac{\eta}{y}}.
        \label{dimensionless_eqmov}
    \end{equation}
    In Appendix A, it is shown that an approximate analytic solution
    of equation (\ref{dimensionless_eqmov}) can be constructed as a non-linear
    average of a ``near'' and a ``far field'' analytic solution
    (see equations \ref{near_far} and \ref{average}),
    which in terms of the respective dimensional variables is
   \begin{equation}
     x_{ws} = \frac{\dot{m}^2}{a \rho_a^2} \left[\left(\frac{a\rho_a}{\dot{m}}\tau \right)^{-\frac{5}{2}} +\left(\frac{3}{2}\right)^{\frac{5}{4}} \left(\frac{a\rho_a}{\dot{m}}\tau\right)^{-\frac{15}{8}}\right]^{-\frac{4}{5}}.
     \label{xwsolm}
   \end{equation}

    \vspace{0.5cm}
  \item { \it Constant ejection density }

    Setting $\rho_0=\mathrm{const.}$ in (\ref{diffeqwsrhoc}) and taking the
      square of the equation, we obtain:
    \begin{equation}
        \frac{\rho_a}{a\rho_0} \left(\frac{dx_{ws}}{d\tau} \right)^2 + x_{ws} - a\tau^2 = 0.
    \end{equation}
    Proposing a power law solution, we find that
    \begin{equation}
      x_{ws} = \frac{a\rho_0}{8\rho_a} \left[-1 + \sqrt{1+ \frac{16 \rho_a}{\rho_0}} \right] \tau^2.
      \label{xxxws}
    \end{equation}
    We therefore find a quadratic dependency on the ejection time.

    Using equations (\ref{timetau}) and (\ref{xxxws}) we then obtain:
    \begin{equation}
        x_{ws} = \frac{8a\beta_0^2 \left( -1+ \sqrt{1+\frac{16}{\beta_0^2}}\right)}{\left[8 + \beta_0^2 \left(-1 + \sqrt{1+\frac{16}{\beta_0^2}}\right) \right]^2}\,t^2,
        \label{solution_ws_time}
    \end{equation}
    where $\beta_0$ is given by equation(\ref{beta0}).
    This implies a quadratic dependence on the evolutionary time.

    In the $\beta_0 \ll 1$ limit, equation (\ref{solution_ws_time}) becomes
    \begin{equation}
        x_{ws} \approx \frac{a}{2}\,\beta_0\,t^2,
    \end{equation}
    and the $\beta_0 \gg 1$ limit leads to
    \begin{equation}
        x_{ws} \approx \frac{a}{4}\,t^2.
    \end{equation}

\end{enumerate}

\section{The numerical simulations}
\label{sec:n-s}
In order to check our analytical solutions we computed a set of 2D axisymmetric simulations for both the constant mass loss rate and constant ejection density cases.
For the simulations we use a code that solves the ideal gas-dynamic (Euler) equations in a fixed two dimensional grid. The code uses a second order Godunov type method with the HLLC \citep{Toro94} approximate Riemann solver, including a linear reconstruction of the primitive variables with the {\it minmod} slope limiter to avoid spurious oscillations. In order to use cylindrical coordinates, the appropriate geometrical source ($\propto 1/r$) terms are included after each timestep in an operator splitting fashion.
Additionally, we incorporated the
  cooling function dependent on density, metallicity and temperature in the energy equation (also as source term) as proposed by \citet{2014MNRAS.440.3100W}, which we compute assuming
  a solar metallicity.

  In order to stabilize the method an artificial diffusion with a (dimensionless)
  value of $0.01$ was added to all of the equations, and a Courant number
  of $0.2$ was used. The code is written in fortran90
  and is paralelized with the {\it  Message Passing Interface.}
We used a computational grid with a size of 0.05 and 0.2 pc along the $r$ and $x$ directions, respectively. The spatial resolution was $\sim 6.89\,\mathrm{au}$, corresponding to $600\times 6000$ cells along the radial and axial directions.

\subsection{Initial and boundary conditions}

We model jets propagating into a uniform, quiescent environment, and consider
the values $n_{a}=100\,$ cm$^{-3}$ and $n_{a}=5000\,$ cm$^{-3}$ for
the environmental density. For convenience, from now on we will use number density instead of mass density. We consider a mean molecular weight of $\mu =1.3$. We also consider the cases of jets with a constant
mass loss rate and a constant ejection density.

In all simulations, we consider the ejection velocity variability
given by equation (\ref{ejec_vel}) with $a= 100\,\mathrm{km\,s^{-1}}$
per millenium$\;=3.17\times10^{-4}\,\mathrm{cm\,s^{-2}}$.
An initial jet radius $r_j=300\,\mathrm{AU}$ (the cross section
being $\sigma=\pi r_j^2$) and a temperature
$T=100\,\mathrm{K}$ (for both the jet and the environment) was imposed
in all models.

We therefore compute four models, two with constant ejection
density (which we label $n_{100}$ and $n_{5000}$, with the subscript
giving the ambient density in cm$^{-3}$) and two with constant
mass loss rate (which we label ${\dot m}_{100}$ and ${\dot m}_{5000}$).
For the models $\dot{m}_{100}$ and $\dot{m}_{5000}$ we choose a total
mass loss rate $\dot{M}=2.24\times10^{-8}\,\mathrm{M_{\odot}yr^{-1}}$,
which corresponds to a mass loss rate at the jet source per unit
area $\dot{m} = 3.36\times 10^{-13}\,\mathrm{g\,cm^{-2}\,s^{-1}}$.
For the models $n_{100}$ and $n_{5000}$, we choose a density
$n_j=1.66\times10^4\,$ cm$^{-3}$ (for a gas with 90\% H and 10\% He).
These model parameters are summarized in Table~\ref{tab:1}.

In the constant mass flux cases, we find a maximum jet density of $n_j=1.56\times 10^7\,\mathrm{cm^{-3}}$, and a minimum jet density of $n_j=6.24\times 10^3\,\mathrm{cm^{-3}}$, corresponding to an ejection time $\tau=1\,\mathrm{yr}$, and $\tau=2500\,\mathrm{yr}$, respectively.

\begin{table}[htbp]
    \centering
    \caption{\footnotesize{Physical conditions for the numerical simulations.}}
    \label{tab:1}
    \begin{tabular}{c|c|c|c}
    \toprule
    & $\dot{M_j}$ & $n_j$ & $n_a$
    \\
    &\scriptsize$\mathrm{(M_{\odot}\,yr^{-1})}$ &\scriptsize$\mathrm{(cm^{-3})}$& \scriptsize$\mathrm{(cm^{-3})}$ \\
    \hline
    \hline
    $n_{100}$ & \multirow{2}{*}{$1.30\times 10^{-7^{\,*}}$}  & \multirow{2}{*}{$1.66\times 10^{4}$}  & $_{100}$ \\
    $n_{5000} $ & & & $_{5000}$ \\
    \cline{1-4}
    $\dot{m}_{100}$ & \multirow{2}{*}{ $3.35\times 10^{-7}$}  & \multirow{2}{*}{$1.55\times 10^{4^{\,*}}$} & $_{100}$  \\
    $\dot{m}_{5000}$  &  & & $_{5000}$ \\
    \bottomrule
    \end{tabular}
    \\
    \vspace{0.25cm}
    \footnotesize{$^*$ These values correspond to an ejection
      time $\tau=1000\,\mathrm{yr}$.}
\end{table}

We should note that
  the parameters we have chosen are inspired in the characteristics
  of ``classical'' HH jets such as HH~34 and HH~111, which have:
  \begin{itemize}
    \item  radii $\sim 100$~AU (Reipurth et al. 2002) and $\sim 300$~AU
      (Reipurth et al. 1997), for HH~34 and HH~111, respectively.
      We chose the larger of these
  two radii in order to have a better resolution of the jet in
  our numerical simulations. These jet radii are measured
  in the ``jet knot chain'' at distances of $\sim 10''$ (corresponding
  to $\sim 6\times 10^{16}$~cm) from the outflow sources, and
  not in the considerably broader HH~34S and HH~111V ``heads'',
\item full spatial velocities (i.e., combining radial velocities
  and proper motions) between 200 and 300~km~s$^{-1}$
  (Hartigan et al. 2001; Reipurth et al. 2002) and lengths (out to
  HH~34S and HH~111V) of $\approx 8\times 10^{17}$~cm. We have then made
  a choice of $a= 100\,\mathrm{km\,s^{-1}}$ per millenium (see above),
  which after an evolutionary time of $\sim 2500$~yr produces a
  jet head with a position and velocity similar to the ones of
  HH~34S and HH~111V (see the follwing section). The simulations
  are stopped at this time (i.e., at 2500~yr),
\item mass loss rates in the
  $10^{-8}\to 10^{-7}$~M$_\odot$~yr$^{-1}$ range (depending on the emission
  lines used for deriving this value and on position along the jets,
  see Podio et al, 2006). We have therefore chosen densities
  that produce mass loss rates of this order of magnitude,
\item even though the ambient densities in HH 34 and HH 111 are of $\sim 10$~cm$^{-3}$ (Raga et al. 1991), we have chosen higher environmental densities
  (of 100 and 5000~cm$^{-3}$, see above) in order to have a substantial
  braking effect, and to obtain larger differences between
  the centre of mass and ram-pressure balance solutions.
\end{itemize}

\section{Results}

We computed four numerical simulations,  two with a constant ejection density (models $n_{100}$ and $n_{5000}$) and two with a constant mass loss rate (models $\dot{m}_{100}$ and $\dot{m}_{5000}$). The top (purple background) and bottom (black background) subpanels of each subfigure presented in Figure \ref{Figure_2} display the numerical density and temperature stratifications at an evolutionary time of 2500 yr for these four models.

Generally, the highest densities and temperatures are found in the
on-axis regions of the jet heads. The jet heads of the models with
constant ejection density (models $n_{100}$ and $n_{5000}$)
have traveled to larger distances from the outflow source than
the ones of constant mass loss rate ($\dot{m}_{100}$ and $\dot{m}_{5000}$).
This qualitative difference is due to the fact that in the constant
mass loss rate models the injected momentum rate ${\dot M}_jv_j$
scales linearly with $v_j$ (and hence also increases linearly
with time, see equation \ref{ejec_vel}), while this scaling is
quadratic (with ejection velocity and with time) in the constant
ejection density models. This means that the constant ejection density jets have a higher momentum content.

As a result of the lower jet beam density and pressure, the constant
mass loss rate models develop an incident/reflected crossing shock
structure, which (in the incident shock region)
leads to the production of faster off-axis motions
in the jet head (see the two bottom panels of Figure~\ref{Figure_2}).
This kind of structure is always found in jet simulations with appropriate parameters, i.e. simulations where the jet is under-dense and under-pressured with respect to the environment during its evolution (see, e.g., Raga 1988; Downes \& Ray 1999).

\begin{figure*}[htbp]
    \centering
    \includegraphics[scale=1]{ev_allmodels_new-0040.eps}
    \caption{\footnotesize{Numerical density and temperature stratifications obtained from our four models (see the text and Table~1), at an evolutionary time of 2500 years. For each model, we show the number density (top half) and temperature (bottom half) in logarithmic color scales.} The axial and radial axes are given in units of $10^{17}$ cm. In the constant mass loss rate models the jet density, and pressure, drops quite considerably leading to the formation of internal ``crossing shocks''.}
    \label{Figure_2}
\end{figure*}

\begin{figure*}[htbp]
    \centering
    \includegraphics[scale=0.55]{n100mlc_linear_new.eps}
    \includegraphics[scale=0.55]{n100rhoc_linear_new.eps}
    \includegraphics[scale=0.55]{n5000mlc_linear_new.eps}
    \includegraphics[scale=0.55]{n5000rhoc_linear_new.eps}
    \caption{\footnotesize{Position of the jet head as a function of time for the constant mass loss rate models (left panels) and the constant ejection density models (right panels). The red solid line shows the ram-pressure balance solution and the green dashed line shows the center of mass solution. The dash-dotted line shows the bow shock and the dotted line the jet shock positions obtained from the numerical simulations.}}
    \label{posvstime}
\end{figure*}

The upper left panel of Figure~\ref{posvstime} shows the position of the bow shock (dash-dotted blue line) and jet shock (dotted blue line) as a function of time for the ${\dot m}_{100}$ constant mass loss rate model. These positions were obtained by searching for the jumps along the symmetry axis of the pressure stratifications for both the bow shock and jet shock. We compare these shock positions with the analytical solution of the center of mass (dashed green line) and ram-pressure balance (solid red line) analytic models. We find that the position of both the bow shock lies
close to the prediction from the center of mass analytic
model, and the jet shock position lies somewhat below.

The upper right panel of Figure \ref{posvstime} shows the time-dependent position of the jet head for the $n_{100}$ constant ejection density model. For times $<1500$~yr, the jet and bow shock positions (obtained from the numerical simulation) quite closely follow the center of mass analytic solution, At larger times, the jet and bow shock positions lie above the center of mass solution, and begin to approach the ram-pressure balance solution. At all times, the positions of the jet and bow shock lie in between the working surface positions predicted by the two analytic models.

The lower left panel of Figure~\ref{posvstime} shows the evolution of the jet head for the ${\dot m}_{5000}$ constant mass loss rate model. For this model, the positions of the jet and bow shock initially follow the center of mass analytic solution, but for times $t>1500$~yr begin to show large deviations, moving slower than the predictions of the two analytic models. These strong deviations are not surprising given the rather extreme departures from a constant cross section, cylindrical jet beam (assumed in the analytic models) found in the numerical simulation (see the bottom panel of Figure 2).

The strong departures from the simple, cylindrical structure assumed in the analytic models occurs because at increasing times the jet density (and pressure) drop quite considerably in this model, leading to the formation of internal ``crossing shocks''. These shocks then lead to the formation of complex off-axis structures in the jet head (see the two bottom panels of Figure \ref{Figure_2}).

The lower right panel of Figure \ref{posvstime} shows the dynamics of the jet head for the $n_{5000}$ constant ejection density model. This model shows jet and bow shock positions that initially follow the center of mass analytic model, and at times $t>1500$~yr tend to approximate the ram-pressure balance model.

Finally, we calculate the relative position differences $\Delta x_a/x_{cm}=
  x/x_{cm}-1$, with $x_{cm}$ being the position of the ``center of mass jet
  head'' (see equations \ref{sol_xcm_t} and \ref{solution_rhoc}).
  The distance from the source $x_a$ corresponds either to the bow shock or
  jet shock positions (obtained from
  the numerical simulations), or to the position of the ``ram-pressure
  balance jet head'' (see equations \ref{xwsolm} and \ref{solution_ws_time}).
  Figure~\ref{fig:difverl} shows the relative position
  differences (with respect to the center of mass position for the jet head)
  for all of our computed models. The first thing to note (see the
  red lines in the two right panels of Figure 4) is that
  the relative position difference between the center of mass and
  ram pressure balance solutions does not strongly depend on time
  for the constant ejection density models (a result that can be
  seen by comparing equations \ref{solution_rhoc} and \ref{xxxws}).
  This result does not strictly hold for the constant mass loss rate models
  (red lines in the left panels of Figure 4).

  For the constant mass loss rate models we see that for $t>550$~yr
  the bow and jet shock positions remain below the center of mass
  working surface position, with negative $\Delta x$ values (see
  the left panels of Figure 4). For the constant ejection density models (right panels
  of Figure 4) we see an initial convergence to the center of
  mass solution of the bow and jet shock positions, and for
  $t>1300$~yr we see progressively larger, positive values for
  $\Delta x/x_{cm}$ (starting to approach the ram-pressure balance
  working surface position).

\begin{figure*}[htbp]
    \centering
    \includegraphics[scale=0.55]{n100mcl.eps}
    \includegraphics[scale=0.55]{n100rhoc.eps}
    \includegraphics[scale=0.55]{n5000mcl.eps}
    \includegraphics[scale=0.55]{n5000rhoc.eps}
    \caption{The relative differences of the ram-pressure balance
      solution (solid red line), the numerical bow shock (dash-dotted
      blue line) and jet shock positions (dotted blue line) with
      respect to the center of mass working surface solution (the dashed,
      green line showing the $\Delta x/x_{cm}=0$ line) as a function
      of time $t$. The results obtained for our four models (see Table 1)
      are shown (constant mass loss rate models on the left, and constant
      ejection density models on the right).}
 \label{fig:difverl}
\end{figure*}

\section{Discussion}

We have presented a first attempt to compare the analytic
ram-pressure balance (Raga \& Kofman 1992) and center of
mass (Cant\'o et al. 2000) analytic formalisms (for obtaining the motions
of working surfaces in variable ejection jets) with axisymmetric
numerical simulations. To this effect, we choose the relatively
simple problem of the leading head of a jet produced with an ejection velocity that increases linearly with time.
We explore the cases of a time-independent mass loss rate
(in which the ejection density scales as the inverse of the
ejection velocity) and of a constant ejection density. From the simulations we compute
the on-axis positions of the jet shock and bow shock, and compare them with the predictions of the
analytic models (in which the jet and bow shock are assumed
to be spatially coincident).

For three of our four simulations (models $n_{100}$, $n_{5000}$ and
${\dot m}_{100}$), a good agreement is found
between the numerical simulations and the analytic models.
In these three models the jet and bow shock positions initially
follow the center of mass analytic solution, and (for two of these
models) at later
times deviate towards the faster moving, ram-pressure balance
analytic model. This result is completely satisfying, as one
would expect that at early times the shocked jet and environment material
will mostly remain within the working surface (as assumed
in the center of mass formalism), and that a substantial
sideways leakage of material (as assumed in the ram-pressure balance
model) should occur at later times. This is discussed in more
detail below.

It is clear that right after a working surface is formed,
the two associated shocks are very close to each other (with a separation
$d\ll r_j$, where $r_j$ is the local jet beam radius).
Because the surface $S\approx 2\pi r_jl$ through which the
gas exits the working surface is then small ($S\ll \pi r_j^2$), most
of the mass going through the shocks remains within the working surface.
At later times,
  the separation between the two shocks increases, and part of the
  mass processed by the working surface shocks is ejected
  sideways into a jet cocoon. For a jet with
    a time-indepenent ejection velocity, the separation $d$ between the two
    shocks attains a steady, maximum value determined by the balance between
    the incoming and outflowing mass rates. Such a balance between
    inflow and outflow from the working surface is
    also approximately attained in the case of a variable velocity
    jet, but the resulting shock separation (and also mass within
    the working surface) still has a (slower) time-dependence.

  It is possible to obtain
  analytic estimates of the separation between the two shocks
  (and therefore the mass of the gas within the working surface)
  with different degrees of complexity (and accuracy), following
  the approaches developed for determining the shock standoff
  distance in blunt body flows (see, e.g., the book of
  Hayes \& Probstein 2003 and Cant\'o \& Raga 1998). Here,
  we present a very simple estimate of the mass within a
  working surface.

  In the stationary state (i.e., after the initial regime of growing
  sepration between the two working surface shocks), the mass
  $M_{ws}$ within the two working surface shocks is:
  \begin{equation}
    M_{ws}\approx {\dot M}_{in}t_s\,,
    \label{mws}
  \end{equation}
  where ${\dot M}_{in}$ is the mass being fed through
  the two shocks and $t_s$ is the timescale for sideways
  leakage of the material. Assuming that the two shocks are
  highly radiative (with a cooling distance\footnote{The cooling distance is usually defined as the distance from a plane-parallel shock to the point where the temperature has dropped to a certain value (see for example Hartigan et al., 1987).} to $\sim 10^3$~K
  much smaller than the jet radius $r_j$), this leakage will
  occur at the post-cooling sound speed
  $c_0=\sqrt{5kT/3\mu m_H}\approx 3$~km~s$^{-1}$ (for $T=10^3$~K and
  $\mu=1.3$),
  and we then have
  \begin{equation}
    t_s= \frac{r_j}{2c_0}\,.
    \label{ts}
  \end{equation}
  Also, the total mass that has been fed into the working surface
  is
  \begin{equation}
    M_{in}\approx {\dot M}_{in}t_d\,,
    \label{min}
  \end{equation}
  where $t_d=x_{ws}/v_{ws}$ (with $x_{ws}$ being the position and
  $v_{ws}$ the velocity of the working surface) is the dynamical
  timescale.

  Combining equations (\ref{mws}-\ref{min}) we obtain\textbf{:}
  \begin{equation}
    \frac{M_{ws}}{M_{in}}\leq {\rm min}
      \left[1,\frac{r_j}{2x_{ws}}\frac{v_{ws}}{c_0}\right]\,.
      \label{mrat}
  \end{equation}
  In this equation, we have considered that the
  $M_{ws}/M_{in}\leq 1$ condition has to be satisfied
  (with $M_{ws}=M_{in}$ for the initial, mass conserving
  regime, and $M_{ws}<M_{in}$ when the sideways ejection from
  the working surface becomes important), while our analytic
  estimate can give values $>1$ because of its implicit assumption
  of a balance between the mass rates into and out of the working
  surface (which is not satisfied in the early evolution of the
  working surface, as discussed above). Clearly, for the
  case of an accelerating working surface (as the one that we
  obtain from our variable ejection models), the derivation of
  equation (\ref{mrat}) should be done evaluating the appropriate
  integrals (see section 4 of the paper of Lim et al. 2002),
  but the simple derivation shown above still gives
  the correct scaling properties.

  We now consider that the transition between the mass
    conserving and efficient mass loosing working surface
    regimes (which can be described with the ``center of mass''
    and ``ram-pressure balance'' equations of motion, respectively)
    occurs when $M_{ws}/M_{in}\sim 1/2$ (i.e., that the
    working surface has lost half of the processed jet material).
    Therefore, the position $x_t$ of the working surface
    at which the transition between the center of mass and ram
    pressure balance regimes occurs, can be estimated setting
    $M_{ws}/M_{in}=1/2$ in equation (\ref{mrat}), which gives:
  \begin{equation}
    x_t=\frac{v_{ws}}{c_0}r_j\,.
    \label{xtran}
  \end{equation}
  Setting $c_0=3$~km~s$^{-1}$ (see above), $r_j=300$~au (see
  section 3) and a velocity $v_{ws}\approx 200$~km~s$^{-1}$,
  we then obtain $x_t\approx 3\times 10^{17}$~cm, This is the correct
  order of magnitude for the location of the
  transition between the ``mass conserving''
  center of mass and the ``efficient mass loosing'' ram-pressure balance
  regimes for the working surface found in our
    $n_{100}$ and $n_{5000}$ numerical simulations (see Figure 3).

To summarize, we find that numerical simulations of a jet with
a linearly accelerating ejection velocity vs. time produce a
leading head which at early times can be analytically
modeled with the center of mass formalism, and (for appropriate flow
parameters) at later
evolutionary times approaches the ram-pressure balance solutions
(the distance from the source at which the
transition between these two regimes occurs can be estimated
with equation (\ref{xtran})).
This result holds provided that during its time-evolution the
jet does not become under dense/under pressured, in which case
strong departures from the analytic solutions are found.

\section{Conclusions}

This paper presents a first attempt to compare analytic solutions
based on both the ``ram-pressure balance'' and ``center of mass''
formalisms with hydrodynamical (axisymmetric) simulations. For
this comparison, we have studied the case of a linearly increasing
ejection velocity as a function of time (and two different forms
for the ejection density, see \S 2).

This simple ejection variability could apply to the early evolution
of an HH outflow, and is relevant for the problem of entrainment
of environmental molecular material into the jet head (see Lim et al. 2002)
and the potential destruction of molecules by shocks, particularly in the jet head. As a direct application to the observed, evolved HH
jets is not straightforward, the interest of our work is primarily
theoretical.

While full analytic solutions for different forms of the ejection variability
have been found in the past with the ``center of mass'' formalism
(see Cant\'o et al, 2000 and Cant\'o \& Raga 2003), no solutions
(except the trivial one for the constant ejection case) have been
previously found with the ``ram-pressure balance'' equation of
motion for the working surfaces. We were surprised to find that
the ``linearly accelerating'' ejection velocity case does lead
to a full, analytic ram-pressure balance solution for the case of
a constant ejection density (though we managed to find only an
approximate analytic solution for the constant mass loss rate
case). This result \textbf{implies} that ram-pressure balance
analytic solutions might also exist for other ejection velocity/density
variability combinations.

Comparing the center of mass and the ram-pressure balance solutions
with axisymmetric numerical simulations (with the same parameters)
we find that for the case of over-dense jets (with
higher densities in the jet side of the leading working surface)
the numerical results initially follow the center of mass
solution, and at later times approach the ram-pressure balance
analytic solution. This result shows a satisfying consistency
between the analytic approaches and the numerical solutions with
the full (axisymmetric) hydrodynamic equations. Not entirely
unexpectedly, we do not
find a good agreement between the analytic and numerical
results for cases in which the jet density drops to values
lower than the environmental density, as the jet beam then develops
internal shocks (in the numerical solutions) that affect the motion
of the leading working surface.

Clearly, it would be interesting to extend this comparison between
numerical simulations and analytic solutions of variable jets
to the case of jets with ``internal working surfaces'' within
their beams. This problem, which would clearly have much more direct
applications to observed HH jets, is left for a future study.

\vspace{0.25cm}
ARa and JC acknowledge support from the DGAPA-UNAM grant IG100218.
We thank an anonymous referee for helpful comments which lead
(among other things) to many of the issues discussed in section 5.


\subsection{Supporting material}
We begin by recalling the definition and properties of quasiconvex functions on $\R$.
We then recall some basic facts about scalar LQR.

\begin{definition}
	A function $f : \cD \mapsto \R\ $ on the convex domain $\ \cD \subseteq \R^n$
	is \emph{quasiconvex} if its sublevel sets
	\(
		\cD_\alpha = \{ x \in \cD : f(x) \leq \alpha \}
	\)
	are convex for all $\alpha \in \R$.
\end{definition}

\begin{lemma}[\citet{boyd-convex}, \S 3.4]
\label{lem:quasiconvex}
	The following facts hold for quasiconvex functions on a convex $\cD \subseteq \R$:
	\begin{enumerate}[label=(\alph*)]
		\item \label{alternatives}
			If $f: \R \mapsto \R$ is continuous,
			then $f$ is quasiconvex
			if and only if at least one of the following conditions holds on $\cD$:
			\begin{enumerate}[label=\arabic*.]
				\item $f$ is nondecreasing.
				\item $f$ is nonincreasing.
				\item \label{alternative-incdec}
					There exists $c \in \cD$ such that
					for all $t \in \cD$,
					if $t < c$ then $f$ is nonincreasing,
					and if $t \geq c$ then $f$ is nondecreasing.
			\end{enumerate}
		\item \label{secondorder}
			If $f : \R \mapsto \R$ is twice differentiable
			and
			$\partial^2 f / \partial x^2 > 0$
			for all $x \in \cD$ where
			$\partial f / \partial x = 0$,
			then $f$ is quasiconvex on $\cD$.
		\item \label{quotient}
			If $f(x) = p(x) / q(x)$,
			where $p: \R \mapsto \R$ is convex with $p(x) \geq 0$ on $\cD$
			and $q: \R \mapsto \R$ is concave with $q(x) > 0$ on $\cD$,
			then $f$ is quasiconvex on $\cD$.
	\end{enumerate}
\end{lemma}

\noindent The following facts about scalar LQR problems can be derived from
the LQR Riccati equation and some calculus (not shown).
\begin{lemma}
\label{lem:scalar-lqr-facts}
For the scalar LQR problem with $a > 0, b > 0$ and $q = r = 1$,
the optimal linear controller $\kopt{a,b}$ is given by the closed-form expression
	\[
		k^\star_{a,b}
		=
		\min_{k \in \R} J_{a,b}(k)
		=
		\min_{k \in \R} \frac{1 + k^2}{-2(a + bk)}
		=
		-\frac{a + \sqrt{a^2 + b^2}}{b}.
	\]
For fixed $a$, the map from $b$ to $k^\star_{a,b}$
is continuous and strictly increasing on the domain $b \in (0, \infty)$
and has the range $(-\infty, -1)$.
For any $k \in (-\infty, -1)$,
there exists a unique $b_k \in (-\infty, -1)$
for which $k = k^\star_{a,b_k}$, given by
\[
	b_k = \frac{2ak}{1 - k^2}.
\]
\end{lemma}

\newcommand{\cD}{\cD}
\newcommand{\cD^b}{\cD^b}
\newcommand{\cD^k}{\cD^k}

\subsection{Proof of \lemmaref{lem:subopt-convex}}

We first recall the statement of the lemma.
\quasi*

\noindent Instead of a monolithic proof,
we present supporting material in \lemmaref{lem:ratio-quasi,lem:cost-quasi}.
We then show the main result in \lemmaref{lem:increasing-main},
which considers $\alpha$-suboptimal neighborhoods on all of $\R$
instead of restricted to $\mathbf{B}$.
\lemmaref{lem:subopt-convex} will follow as a corollary.

We proceed with more setup.
Recall that the scalar DDF problem\ is defined by
$\mathbf{A} = \{a\}$ and $\mathbf{B} = [\recip \theta, 1]$, where $a > 0$.
For this section,
let
\[
	\cD
	= \{ (b, k) \in (0, \infty) \times \R : a + bk < 0 \}
\]
(note that $J_b(k) < \infty \iff a + bk < 0$).
Denote its projections by
$\cD^b(k) = \{ b: (b, k) \in \cD \}$
and
$\cD^k(b) = \{ k: (b, k) \in \cD \}$.
We compute the suboptimality ratio $r : \cD \mapsto \R$ by
\[
	r(b, k) =
	\frac{J_b(k)}{J^\star_b}
	=
	\nofrac{
		\frac{
			1 + k^2
		}{
			-2(a + bk)
		}
	}{
		\frac{
			a + \sqrt{a^2 + b^2}
		}{
			b^2
		}
	}
	=
	\frac{
		(1 + k^2) b^2
	}{
		-2(a + bk)(a + \sqrt{a^2 + b^2})
	}.
\]
We denote its sublevel sets with respect to $b$ for fixed $k$ by
\[
	\cD^b_\alpha(k) = \{ b \in \cD^b(k) : r(b, k) \leq \alpha \}.
\]
\begin{lemma}
\label{lem:ratio-quasi}
	For fixed $k < 0$, the ratio $r(b, k)$ is quasiconvex on $\cD^b_k$,
	and there is at most one $b \in \cD^b_k$ at which $\partial r / \partial b = 0$.
\end{lemma}
\begin{proof}
	By inspection, $r(b, k)$ is smooth on $\cD^b$.
	We now show that the second-order condition of \lemmaref{lem:quasiconvex}\ref{secondorder} holds.
	\input{ratio_intervals_proof.tex}
	The conclusion follows from \lemmaref{lem:quasiconvex}\ref{secondorder}.
\end{proof}

\begin{lemma}
\label{lem:cost-quasi}
	For fixed $b$, the cost $J_b(k)$ is quasiconvex on $\cD^k(b)$.
	Also, $J_b(k)$ is not monotonic,
	so case \ref{alternative-incdec} of \lemmaref{lem:quasiconvex}\ref{alternatives} applies.
\end{lemma}
\begin{proof}
	We have
	\[
		J_b(k)
		=
		\frac{
			1 + k^2
		}{
			-2(a + bk)
		}.
	\]
	The numerator is nonnegative and convex on $k \in \R$.
	The denominator is linear (hence concave) and positive on $\cD^k(b)$.
	Quasiconvexity follows from \lemmaref{lem:quasiconvex}\ref{quotient}.
	Nonmonotonicity follows from the fact that $J_b(k)$ is smooth on $\cD^k(b)$
	and has a unique optimum at $\kopt{b}$, which is not on the boundary of $\cD^k(b)$.
\end{proof}

\noindent We now combine these into the main result.

\begin{lemma}
\label{lem:increasing-main}
	For a scalar DDF problem,
	if $\alpha > 1$ and $k < -1$,
	then $\cD^b_\alpha(k)$ is either:
	a bounded closed interval $[b_1, b_2]$,
		with $b_1$ and $b_2$ increasing in $k$, or
	a half-bounded closed interval $[b_1, \infty)$,
		with $b_1$ increasing in $k$.
\end{lemma}

\begin{proof}
	By \lemmaref{lem:ratio-quasi}, due to quasiconvexity $\cD^b_\alpha$ is convex.
	The only convex sets on $\R$ are the empty set and all types of intervals:
	open, closed, and half-open.
	We know $\cD^b_\alpha$ is not empty because it contains $b_k$.
	We can further assert that
	$\cD^b_\alpha$ has a closed lower bound
	because $\lim_{b \to (-a/k)} r(b, k) = \infty$
	(see \citet[\S A.3.3]{boyd-convex} for details).
	However, the upper bound may be closed or infinite.
	We handle the two cases separately.

	\paragraph{Bounded case.}
	Fix $k_0 < -1$.
	Suppose $\cD^b_\alpha(k_0) = [b_1, b_2]$ for $0 < b_1 < b_2 < \infty$.
	By the implicit function theorem (IFT), at any $(b_0, k_0)$ satisfying $r(b_0, k_0) = \alpha$,
	if $\partial r / \partial b|_{b_0, k_0} \neq 0$
	then there exists an open neighborhood around $(b_0, k_0)$
	for which the solution to $r(b, k) = \alpha$ can be expressed as
	$(g(k), k)$, where $g$ is a continuous function of $k$ and
	\[
		\left. \frac{\partial g(k)}{\partial k} \right|_{k_0} =
		\left. - \left( \frac{\partial r}{\partial b} \right)^{-1}
		\frac{\partial r}{\partial k} \right|_{b_0, k_0}.
	\]
	By the continuity and quasiconvexity of $r$,
	and the fact that $\partial r / \partial b = 0$ only at $b_k$
	(\lemmaref{lem:ratio-quasi})
	we know that $r(b_1, k_0) = r(b_2, k_0) = \alpha$
	and
	\[
		\left. \frac{\partial r}{\partial b} \right|_{b_1, k_0} < 0
		\quad \text{and} \quad
		\left. \frac{\partial r}{\partial b} \right|_{b_2, k_0} > 0.
	\]
	By \lemmaref{lem:scalar-lqr-facts},
	since $k < -1$
	there exists $b_k > 0$ satisfying
	$k = k^\star_{b_k}$.
	Since $r(b_k, k) = 1$ and $\alpha > 1$, we know $b_{k_0} \in (b_1, b_2)$.
	Again by \lemmaref{lem:scalar-lqr-facts},
	the map from $b$ to $k^\star_b$ is increasing in $b$.
	Therefore,
	\(
		{k^\star_{b_1} < k_0 < k^\star_{b_2}}.
	\)
	By the quasiconvexity and nonmonotonicity of $J_b(k)$ from \lemmaref{lem:cost-quasi},
	via \lemmaref{lem:quasiconvex}\ref{alternatives}
	we have
	\[
		\left. \frac{\partial r}{\partial k} \right|_{b_1, k_0} \geq 0
		\quad \text{and} \quad
		\left. \frac{\partial r}{\partial k} \right|_{b_2, k_0} \leq 0.
	\]
	Therefore, the functions $g_1, g_2$ satisfying the conclusion of the IFT
	in the neighborhoods around $(b_1, k_0)$ and $(b_2, k_0)$ respectively
	also satisfy
	\[
		\left. \frac{ \partial g_1(k)}{ \partial k } \right|_{b_1, k_0} \geq 0
		\quad \text{and} \quad
		\left. \frac{ \partial g_2(k)}{ \partial k } \right|_{b_2, k_0} \geq 0.
	\]
	Therefore, $b_1$ and $b_2$ are locally nondecreasing in $k$.

	\paragraph{Unbounded case.}
	Suppose $\cD^b_\alpha(k) = [b_1, \infty)$ for $b_1 < \infty$.
	By the same IFT argument as in the bounded case, $b_1$ is increasing in $k$.
	By the quasiconvexity of $r$ in $b$,
	the value of $r$ is increasing for $b  > b_k$,
	but the definition of $\cD^b_\alpha(k)$ implies that
	$r(b, k) \leq \alpha$ for all $b > b_k$.
	Therefore, $\lim_{b \to \infty} r(b, k)$ exists
	and is bounded by $\alpha$.
	In particular,
	\begin{equation*}
		\begin{split}
			\lim_{b \to \infty} r(b, k)
			&=
			\lim_{b \to \infty}
			\frac{
				(1 + k^2) b^2
			}{
				-2(a + bk)(a + \sqrt{a^2 + b^2})
			}
			\\ &=
			\lim_{b \to \infty}
			\frac{
				(1 + k^2) b^2 / b^2
			}{
				-2(a + bk)(a + \sqrt{a^2 + b^2}) / b^2
			}
			\\ &=
			-
			\frac{
				1 + k^2
			}{
				2k
			}.
		\end{split}
	\end{equation*}
	Taking the derivative shows that this value is decreasing in $k$ for $k < 0$.
	Therefore, if $k < k' < 0$ then
	\[
		\lim_{b \to \infty} r(b, k') \leq \lim_{b \to \infty} r(b, k) \leq \alpha.
	\]
	The property that $r(b, k')$ is increasing in $b$ for $b > b_k$
	further ensures that
	$r(b, k') \leq \alpha$ for all $b > b_k$.
	Therefore, $\cD^b_\alpha(k')$ is also unbounded.
\end{proof}

\noindent For completeness, we prove \lemmaref{lem:subopt-convex}.

\quad

\begin{proof} (of \lemmaref{lem:subopt-convex}).
	By \lemmaref{lem:increasing-main},
	$\cD^b_\alpha(k)$ is either
	a bounded closed interval $[b_1, b_2]$,
		with $b_1$ and $b_2$ increasing in $k$, or
	a half-bounded closed interval $[b_1, \infty)$,
		with $b_1$ increasing in $k$.
	Recall that $\suboptneighb{\alpha}{k} = \cD^b_\alpha(k) \cap \mathbf{B}$
	with $\mathbf{B} = [\recip \theta, 1]$.
	Therefore, the half-bounded case can be reduced to the bounded case with $b_2 = 1$.
	The intersection can be expressed as
	\[
		\suboptneighb{\alpha}{k} = [\max\{b_1, \textstyle \recip \theta\}, \min\{b_2, 1\}],
	\]
	where the interval $[a, b]$ is defined as the empty set if $a > b$.
	Taking the maximum or minimum of a nonstrict monotonic function and a constant
	preserves the monotonicity, so we are done.
\end{proof}

\section{Introduction}
\input{introduction.tex}

\section{Problem setting}
\input{setting.tex}

\paragraph{Standard LQR problem}
\input{lqr.tex}

\section{Theoretical results}
\input{results_scalar.tex}

\section{Empirical results}
\input{empirical_upperbound.tex}

\subsection{Suboptimal neighborhood visualizations}
\input{empirical_lowerbound.tex}

\section{Related work}
\input{related.tex}

\section{Conclusion and future work}
\input{conclusion.tex}

\newpage

\subsection{Scalar upper bound}

We are now ready to bound the covering number for scalar systems.
The first lemma bounding $J^\star_{a,b}$ will be useful for the lower bound also.
We then construct a cover inductively.

\begin{lemma}
\label{lem:scalar-cost-ub}
	In a scalar LQR problem,
	if $a > 0$ and $0 < b \leq 1$,
	then the optimal scalar LQR cost satisfies the bounds
	\(
		\nofrac{2a}{b^2}
		<
		J^\star_{a,b}
		<
		\nofrac{(2a + 1)}{b^2}.
	\)
\end{lemma}
\begin{proof}
	The lower bound is visible from
	the closed-form solution for the scalar Riccati equation, which is
	\(
		J^\star_{a,b} = \frac{a + \sqrt{a^2 + b^2}}{b^2}.
	\)
	The upper bound is obtained by substituting $a^2 + b^2 \leq (a + 1)^2$.
\end{proof}

\begin{lemma}
\label{lem:scalar-cover-recursion}
	If $p, \tau, k = {\normalfont \texttt{GCC}}(a, b_1, b_2, q)$,
	then for any $\beta \in (0, 1)$,
	there exists $k' \in \R$
	such that
	$\displaystyle p', \tau, k' = {\normalfont \texttt{GCC}} \left( a, \beta b_1, \beta b_2, \beta^{-2} q \right)$,
	where $\ p' = \beta^{-2} p$.%
\end{lemma}
\begin{proof}
	In the scalar system,
	the GCC matrix Riccati equation \eqref{eq:petersen-riccati}
	reduces to the quadratic equation
	\begin{equation}
	\label{eq:scalar-ric}
		\left(
			\textstyle \recip \tau b_1^2
			- \recip{1 + \tau} b_2^2
		\right)
		p^2
		+ 2a p
		+ q
		= 0.
	\end{equation}
	Substituting $p' = \beta^{-2} p$ into \eqref{eq:scalar-ric}
	and multiplying by $\beta^{-2}$ yields
	a new instance of~\eqref{eq:scalar-ric}
	with the parameters
		$b_1' = \beta b_1$, \
		$b_2' = \beta b_2$, \
		$q' = \beta^{-2} q$,
	for which $p'$ is a solution
	with $\tau$ unchanged.
\end{proof}

\begin{theorem}
\label{thm:covering-scalar}
	For the scalar DDF problem\ defined by
		$\mathbf{A} = \{ a \}$,
		where $a > 0$,
		and
		$\mathbf{B} = \left[\textstyle \recip \theta, 1 \right]$,
	if $\alpha \geq \frac{2a+1}{2a}$,
	then $\coveringnum{\alpha}{\mathbf{B}} = O(\log \theta)$.
\end{theorem}
\begin{proof}
	We construct a cover from the upper end of $\mathbf{B}$.
	By \lemmaref{lem:scalar-cost-ub},
	the condition ${\alpha \geq \frac{2a+1}{2a}}$ implies that
	\(
		J^\star_{b=1} < \alpha 2 a < \alpha J^\star_{b=1}.
	\)
	Therefore, by
	\lemmaref{lem:petersen-gcc,lem:petersen-existence},
	there exists $\beta \in (0, 1)$
	and
	$p, \tau, k$
	such that
	$p, \tau, k = {\normalfont \texttt{GCC}}(a, (1-\beta)/2, (1+\beta)/2, 1)$
	and
	$p \leq \alpha 2a$.

	Proceeding inductively,
	suppose that for $N \geq 1$,
	we have covered $[\beta^N, 1]$
	by the intervals $\mathbf{B}_n = [\beta^{n+1}, \beta^{n}]$
	for $n \in \{0, \dots, N-1\}$,
	and each $\mathbf{B}_n$ has a controller $k_n$ such that
	\[
		\beta^{-2n}p, \tau, k_n = {\normalfont \texttt{GCC}} \left(
			a,
			\nofrac{(\beta^{n} - \beta^{n+1})}{2},
			\nofrac{(\beta^{n} + \beta^{n+1})}{2},
			\beta^{-2n}
		\right).
	\]
	Then the existence of the desired $\mathbf{B}_N, k_N$
	follows immediately from
	\lemmaref{lem:scalar-cover-recursion}.

	By \lemmaref{lem:lancaster-ARE-domination},
	for each $\mathbf{B}_n$
	the GCC state cost $q_n = \beta^{-2n} \geq 1$
	is an upper bound on the cost if we replace $q_n$ with $1$
	to match the DDF problem.
	Therefore,
	for each interval $\mathbf{B}_n$,
	for all $b \in \mathbf{B}_n$,
	\[
		\alpha J^\star_{b}
		\geq
		\alpha J^\star_{\beta^n}
		>
		\beta^{-2n} \alpha 2 a
		\geq
		\beta^{-2n} p
		\geq
		J_b(k_n),
	\]
	where first inequality is due to \lemmaref{lem:lancaster-ARE-domination},
	the second is due to \lemmaref{lem:scalar-cost-ub},
	the third is by construction of $p$,
	and last is due to the GCC guarantee of $k_n$.
	Hence, $\mathbf{B}_n \subseteq \suboptneighb{\alpha}{k_n}$.
	We cover the full $\mathbf{B}$ when
	$\beta^{N} \leq \recip \theta$,
	which is satisfied by $N \geq - \log \theta / \log \beta$.
\end{proof}

\subsection{Scalar lower bound}

For the matching lower bound,
we begin by
deriving a simplified overestimate of $\suboptneighb{\alpha}{k}$.
We then show that the true $\suboptneighb{\alpha}{k}$
is still a closed interval moving monotonically with $k$.
Finally, we argue that
the gaps between consecutive elements of a cover grow at most geometrically,
while the range of $k$ values in a cover must grow linearly with $\theta$.

\begin{lemma}
\label{lem:scalar-neighborhood-optimistic}
	For a scalar DDF problem\ with $a \geq 1,\ \mathbf{B} = [\recip \theta, 1]$,
	for any $k < 0$,
	if $\alpha \geq 3/2$, then
	$\suboptneighb{\alpha}{k} \subseteq \recip{|k|} [c_1 \!- \!c_2,\ c_1\! + \!c_2]$,
	where $c_1$ and $c_2$ are constants depending on $\alpha$ and $a$.
\end{lemma}
\begin{proof}
	Beginning with the closed-form solution for $J_b(k)$,
	which can be derived from \eqref{eq:cost-suboptimal},
	we define
	\begin{equation}
	\label{eq:Jbk}
		J_b(k) = \frac{1 + k^2}{-2(a + bk)}
		\geq \frac{k^2}{-2(a + bk)}
		\defeq \underline{J_b}(k)
		.
	\end{equation}
	By \lemmaref{lem:scalar-cost-ub}, we have
	\(
		J^\star_b < \nofrac{3a}{b^2} \defeq \overline{J^\star_b}
	\),
	so $\tilde r = \nofrac{\underline{J_b}(k)}{\,\overline{J^\star_b}}$
	is a lower bound on the suboptimality of $k$.
	Computing $\partial^2 \tilde r / \partial b^2$
	shows that $\tilde r$ is strictly convex in $b$ on the domain $a + bk < 0$,
	so the $\alpha$-sublevel set of $\tilde r$ is the closed interval
	with boundaries where $\tilde r = \alpha$.
	This equation is quadratic in $b$ with the solutions
	\(
		b = -a(3 \alpha \pm \sqrt{9 \alpha^2 - 6 \alpha})/k
	\).
	The resulting interval contains $\suboptneighb{\alpha}{k}$.
\end{proof}

\begin{restatable}{lemma2}{quasi}
\label{lem:subopt-convex}
	For a scalar DDF problem, if $\alpha > 1$
	and $k < -1$,
	then $\suboptneighb{\alpha}{k}$ is either empty or a closed interval
	$[b_1, b_2]$,
	with $b_1$ and $b_2$ positive and nondecreasing in $k$.
\end{restatable}
\begin{proof}
The result follows from quasiconvexity of
both the suboptimality ratio $J_b(k) / J^\star_b$
and the cost $J_b(k)$.
Showing these requires some tedious calculations.
\ifextended{
	For details, see \Cref{sec:appendix-theory}.
}{
	For details, see the extended version of this paper.
}
\end{proof}

\begin{theorem}
\label{thm:scalar-lb}
	For a scalar DDF problem\ with
	$a = 1,\ \mathbf{B} = [\recip \theta, 1]$,
	if $\alpha \geq \frac 3 2$, %
	then
	$\coveringnum{\alpha}{\mathbf{B}} = \Omega(\log \theta)$.
\end{theorem}
\begin{proof}
	From the closed-form solution
	\(
		k^\star_{a,b} = -(a + \sqrt{a^2 + b^2})/b
	\),
	we observe that $k^\star_b < -1$ for all $b \in \mathbf{B}$.
	This, along with the quasiconvexity of $J_b(k)$ in $k$,
	implies that there exists a minimal $\alpha$-suboptimal cover $\cC$ for which all $k_i < -1$.
	Suppose $\cC = k_1, \dots, k_N$ is such a cover, ordered such that $k_i < k_{i+1}$.
	Then
	by \lemmaref{lem:subopt-convex},
	$\suboptneighb{\alpha}{k_i}$
	and
	$\suboptneighb{\alpha}{k_{i+1}}$
	must intersect,
	so their overestimates according to
	\lemmaref{lem:scalar-neighborhood-optimistic}
	certainly intersect,
	therefore satisfying
	\[
		\frac{c_1 + c_2}{-k_{i+1}}
		\geq
		\frac{c_1 - c_2}{-k_i}
		\implies
		\frac{k_{i+1}}{k_i} \leq \frac{c_1 + c_2}{c_1 - c_2}
		\implies
		\frac{k_N}{k_1} \leq \left(\frac{c_1 + c_2}{c_1 - c_2}\right)^{N-1}.
	\]
	By \lemmaref{lem:scalar-neighborhood-optimistic},
	to cover $b = 1$
	controller $k_1$ must satisfy
	$k_1 \geq -(c_1 + c_2)$,
	and to cover $b = \recip \theta$,
	controller $k_N$ must satisfy
	$k_N \leq -\theta(c_1 - c_2)$.
	Along with the previous result, this implies
	\[
		\left(\frac{c_1 + c_2}{c_1 - c_2}\right)^{N-1}
		\geq
		\theta
		\frac{c_1 - c_2}{c_1 + c_2}
		\implies
		N \geq \frac{\log \theta}{\log \frac{c_1 + c_2}{c_1 - c_2}}.
	\]
	Recalling that $c_1$ and $c_2$ only depend on $a$ and $\alpha$,
	the $\Omega(\log \theta)$ dependence on $\theta$ is established.
\end{proof}

\paragraph{Remarks}
\begin{itemize}
	\item For the upper bound, it may be possible to compute or bound $\beta$ in the scalar case
		as a function of $a$ and $\alpha$,
		but the analogous result will likely be much more complicated in the matrix case.
	\item \theoremref{thm:covering-scalar} imposes a lower bound on $\alpha$ greater than $1$.
		We believe this is a mild condition in practice:
		if the application demands a suboptimality ratio very close to 1,
		then the size of the suboptimal cover is likely to become impractical for storage.
		However, further theoretical results building upon suboptimal coverings may require eliminating the bound.
\end{itemize}


\section{Introduction}

A \textbf{translation surface} is a closed Riemann surface with an additional structure which can be described by certain gluing data. We construct such a surface as finitely many polygons in the Euclidean plane with edge identifications along pairs of parallel edges. If all polygons are unit squares one obtains an \textbf{origami} (also known as \textbf{square-tiled surface}). Each origami naturally defines a torus cover sending each square in the tiling to the torus. We are interested in the case where this cover is normal. Then the symmetry group of the origami is maximal and we call the origami \textbf{regular}. A regular origami is completely determined by its deck transformation group $G$ and two deck transformations $x$ and $y$ mapping a fixed square to its right and upper neighbor, respectively (see e.g. \cite{flake2020strata}). We denote such an origami by the tuple $(G,x,y)$.

The matrix group $\operatorname{SL}(2,\mathbb{R})$ acts on translation surfaces by sheering the polygons in the Euclidean plane. Sometimes the orbit defines an algebraic curve in the moduli space of complex algebraic curves called \textbf{Teichmüller curve}. The stabilizer under this action captures whether this happens and - assuming a positive answer - much of the geometry of the Teichmüller curve. For a translation surface $X$, the stabilizer is called the \textbf{Veech group} of $X$ and is denoted by $\operatorname{SL}(X)$. Origamis define always Teichmüller curves. The Veech groups of \textbf{reduced origamis}, i.e., the geodesic segments between singularities span $\mathbb{Z}^2$, are finite index subgroups of $\operatorname{SL}(2,\mathbb{Z})$. Here, a singularity means a vertex of a square in the tiling with cone angle larger than $2\pi$. Since non-trivial regular orgiamis are reduced, we restrict ourselves to studying the $\operatorname{SL}(2,\mathbb{Z})$-action. On a regular origami $\mathcal{O}=(G,x,y)$ this action is defined by $S\cdot \mathcal{O}=(G,y^{-1},x)$ and $T\cdot \mathcal{O}=(G,x,yx^{-1})$, where $S=\left(\begin{smallmatrix}0&-1\\1&0\end{smallmatrix}\right),~ T=\left(\begin{smallmatrix}1&1\\0&1\end{smallmatrix}\right).$ For more details, see e.g. \cite{schmithuesen2004algorithm} and \cite{WS-deficiency}.

We are interested in the following open questions: for which origamis are the Veech groups congruence subgroups and for which are they far away from being a congruence subgroup? Weitze-Schmithüsen showed in \cite{schmithuesen2004algorithm} that almost all congruence groups occur as Veech groups. However, Hubert and Leliévre proved that in the stratum $\mathcal{H}(2)$ all but one of the occurring Veech groups are not congruence groups (see \cite{HubertLelievre06}). In \cite{WS-deficiency}, Weitze-Schmithüsen introduced the deficiency of finite index subgroups of $\operatorname{SL}(2,\mathbb{Z})$. It measures how far the group is from being a congruence subgroup. She also established the notion of totally non-congruence groups. Such a group projects surjectively onto $\operatorname{SL}(2,\mathbb{Z}/n\mathbb{Z})$ for each $n\in \mathbb{Z}_+$, i.e., no information about the group itself can be recovered from the images under these natural projections. In \cite{sp-ws-2018}, an infinite family of origamis with totally non-congruence subgroups as Veech groups are constructed for each stratum. These origamis had only few symmetries. In this article, we present sufficient conditions for regular origamis to have totally non-congruence subgroups as Veech groups and introduce a class of regular origamis satisfying this condition.

A version of this paper will be published in the proceedings of the 3rd BYMAT Conference 2020.

\section{Prerequisites and preliminary results}\label{sec prerequisites}

In this section, we introduce basic concepts and preliminary results, which are used in \Cref{sec application to regular origamis}. Note that the Euclidean metric on $\mathbb{R}^2$ lifts to a metric on a translation surface. Therefore, notions as directions and geodesics are well-defined on translation surfaces. A \textbf{cylinder} on a translation surface is a maximal collection of parallel closed geodesics. Given a cylinder on a translation surface there exist $w,h>0$ such that the cylinder is isometric to a Euclidean cylinder $\mathbb{R}/w\mathbb{Z}\times (0,h)$. One calls $w$ the \textbf{circumference}, $h$ the \textbf{height}, and the quotient $\frac{h}{w}$ the \textbf{modulus} of the cylinder. If the genus of the translation surface is larger than one, a cylinder is bounded by geodesics between singularities. We call such a geodesic a \textbf{saddle connection}. The direction of a saddle connection bounding a cylinder is called the \textbf{direction of the cylinder}. A \textbf{cylinder decomposition} is a collection of pairwise disjoint cylinders such that the union of their closures covers the whole surface.

The cylinder decomposition of an origami leads to a parabolic element in its Veech group given the following situation. Let $\mathcal{O}$ be an origami, $v\in \mathbb{Z}^2$ be a rational direction, and $A\in \operatorname{SL}(2,\mathbb{Z})$ be a matrix mapping $e_1 = \tbinom{1}{0}$ to $v$. If $\mathcal{O}$ decomposes into cylinders $C_1,\dots,C_k$ with inverse moduli $m_1,\dots,m_k$ and $m$ is the smallest common integer multiple of all the $m_i$, then the matrix
$A\cdot
\bigl(\begin{smallmatrix}
    1 & m\\
    0 & 1
\end{smallmatrix}\bigr)
\cdot A^{-1}$ is contained in the Veech group $\operatorname{SL}(\mathcal{O})$ (see e.g. \cite[Section 2.1]{WS-deficiency}).

\begin{center}
\begin{figure}
  \begin{minipage}[c]{0.27\textwidth}
    \hspace{0.45cm}\input{ex-cylinderdecomp}
  \end{minipage}\hfill
  \begin{minipage}[c]{0.72\textwidth}\vspace{0.35cm}
    \caption{Let $D_8$ denote the dihedral group $\langle r,s ~|~ r^{4}=s^{2}=1,~s\inv rs=r\inv  \rangle.$ The clyinder decomposition of the $2$-origami $\mathcal{O}=(D_8,r,s)$ in horizontal direction consists of two cylinders shaded in green and blue, respectively. The inverse modulus of both cylinders is $4$. Choosing $A$ as the identity element, we obtain that the parabolic matrix $\bigl(\begin{smallmatrix}
    1 & 4\\
    0 & 1
\end{smallmatrix}\bigr)$ lies in the Veech group $\operatorname{SL}(\mathcal{O})$.
    } \label{fig::example-cylinderdecomp}
  \end{minipage}
\end{figure}
\end{center}

We conclude this section giving a sufficient condition when finite index subgroups of $\operatorname{SL}(2, \mathbb{Z})$ are totally non-congruence groups. It is used in the proofs of \Cref{cor. cond. totally NCGp} and \Cref{(abc)-gps totally NCG}.

\begin{thm}[{\cite[Theorem 1]{sp-ws-2018}}]\label{thm sp-ws tot NCG}
 Let $\Gamma$ be a finite index subgroup of $\operatorname{SL}(2, \mathbb{Z})$. Suppose that for each prime $p$ there exist matrices $A_1 , A_2 \in \operatorname{SL}(2, \mathbb{Z})$ with the following properties:
 \begin{enumerate}
     \item $\forall j \in \mathbb{Z}_+ : A_1 e_1 \neq j \cdot A_2 e_1$ modulo $p$.
     \item There exist $m_1 , m_2 \in \mathbb{Z}_+$ with $A_1T^{m_1}A_1^{-1},~A_2T^{m_2}A_2^{-1}\in \Gamma$ such that $p$ divides neither $m_1$ nor $m_2$.
 \end{enumerate}
Then $\Gamma$ is a totally non-congruence group.
\end{thm}

\section{Application to regular origamis}\label{sec application to regular origamis}

In this section, we use cylinder decompositions in different directions to construct the parabolic matrices occurring in \Cref{thm sp-ws tot NCG}. The moduli of the considered cylinders coincide with the order of certain deck transformations. We apply \Cref{thm sp-ws tot NCG} to regular origamis with deck transformation groups that admit for each prime number a cylinder decomposition with a suitable modulus. The following lemma computes the inverse moduli of the cylinders in the directions of interest.

\begin{lemma}\label{lem moduli of cylinders}
    Let $\mathcal{O}=(G,x,y)$ be a regular origami. For $m\in \mathbb{Z}_{\ge0}$, the inverse modulus of all cylinders in direction $\tbinom{1}{-m}$ coincides with the order of $xy^m$.
\end{lemma}

\begin{proof}
    Denote $\tbinom{1}{-m}$ by $v$. Acting with the matrix
    $A=\bigl(
    \begin{smallmatrix}
    1 & 0\\
    -m & 1
    \end{smallmatrix}
    \bigr)=\left( S^{3}TS\right)^m\in \operatorname{SL}(2,\mathbb{Z})$ maps the horizontal direction to the direction $v$, i.e., $A\cdot e_1=v$. The inverse modulus of all horizontal cylinders of the origami $A\cdot\mathcal{O}=(G,xy^m,y)$ coincides with the order of $xy^m$. Note that acting by matrices in $\operatorname{SL}(2,\mathbb{Z})$ does not change the modulus of a cylinder. Hence, the inverse modulus of the cylinder in direction $v$ of the origami $\mathcal{O}$ equals the order of $xy^m$.
\end{proof}

Using \Cref{thm sp-ws tot NCG} and \Cref{lem moduli of cylinders}, we deduce a sufficient condition for regular origamis to have a totally non-congruence group as Veech group.

\begin{proposition}\label{cor. cond. totally NCGp}
    Let $\mathcal{O}=(G,x,y)$ be a regular origami. If for each prime $p$ one of the following holds
    \begin{enumerate}
        \item $\gcd(p,\operatorname{ord}(y)\cdot\operatorname{ord}(yx))=1$ or
        \item there exist $m_1,m_2\in \mathbb{Z}_{\ge0}$ with $\overline{m}_1\not\equiv \overline{m}_2 \mod p$ and $\gcd(p,\operatorname{ord}(xy^{-m_1})\cdot\operatorname{ord}(xy^{-m_2}))=1$.
    \end{enumerate}
    Then the Veech group $\operatorname{SL}(\mathcal{O})$ is a totally non-congruence group.
\end{proposition}

\begin{proof}
    Fix a prime $p$. If condition (1) holds, consider the matrices $S^{-1}T^{-1}=\bigl(
    \begin{smallmatrix}
    0 & 1\\
    -1 & 1
    \end{smallmatrix}
    \bigr)$ and $TS^{-1}=\bigl(
    \begin{smallmatrix}
    -1 & 1\\
    -1 & 0
    \end{smallmatrix}
    \bigr).$  We obtain $S^{-1}T^{-1}\cdot\mathcal{O}=(G,yx,x^{-1})$ and $TS^{-1}\cdot\mathcal{O}=(G,y,x^{-1}y^{-1})$. The moduli of the horizontal cylinders of the regular origamis $(G,yx,x^{-1})$ and $(G,y,(yx)^{-1})$ are $\operatorname{ord}(yx)=\vcentcolon a$ and $\operatorname{ord}(y)=\vcentcolon b$, respectively. Hence, $S^{-1}T^{-1} \cdot T^a\cdot TS$ and $TS^{-1}\cdot T^b \cdot ST^{-1}$ lie in the Veech group $\operatorname{SL}(\mathcal{O})$. Moreover, we obtain for each $j\in\mathbb{Z}$ the inequality $S^{-1}T^{-1}\cdot e_1=\tbinom{0}{-1}\neq j\cdot \tbinom{-1}{-1}= j\cdot TS^{-1}\cdot e_1$ modulo $p$.

    If condition (2) holds, then let $m_1,m_2$ be natural numbers satisfying condition (2). Define the matrices
    $A_1=\bigl(
    \begin{smallmatrix}
    1 & 0\\
    m_1 & 1
    \end{smallmatrix}
    \bigr),~
    A_2=\bigl(
    \begin{smallmatrix}
    1 & 0\\
    m_2 & 1
    \end{smallmatrix}
    \bigr).$ Since ${m}_1\not\equiv {m}_2 \mod p$, we have $A_1 e_1 \neq j \cdot A_2 e_1$ modulo $p$ for each $j\in\mathbb{Z}_+$.

    As $\gcd(p,\operatorname{ord}(xy^{-m_1})\cdot\operatorname{ord}(xy^{-m_2}))=1$, set $k_1=\operatorname{ord}(xy^{-m_1})$ and $k_2=\operatorname{ord}(xy^{-m_2})$. Using \Cref{lem moduli of cylinders}, we conclude that the matrices $A_iT^{k_i}A_i^{-1}$ are contained in the Veech group of the origami $\mathcal{O}$. By \Cref{thm sp-ws tot NCG}, the claim follows.
\end{proof}

In the following corollary, we construct generating sets $\{x,y\}$ of alternating groups $A_n$ satisfying the conditions given in \Cref{cor. cond. totally NCGp}. Consequently, the infinite family of regular origamis $(A_n,x,y)$ have totally non-congruence groups as Veech groups.

\begin{corollary}\label{A_n totally NCG}
    For each prime $n\ge 5$, the regular origami $(A_n, (1,2,3), (1,2,3,\dots,n))$ has a totally non-congruence group as Veech group.
\end{corollary}

\begin{proof}
    Set $x\vcentcolon=(1,2,3)$ and $y\vcentcolon=(1,2,3,\dots,n)$. For each prime $p\neq n$, we consider the group elements $y$ and $yx$. Since the orders of $y$ and $yx$ are equal to $n$, the prime $p$ does not divide $\operatorname{ord}(y)\cdot\operatorname{ord}(yx).$

    For the prime $n$, we consider the group elements $xy^{n-1}$ and $x$, i.e., $m_1=1-n$ and $m_2=0$. Note that ${1-n}\not\equiv 0\mod n$. The permutation $xy^{n-1}$ has the fixed point $2$ and thus $n$ does not divide the order of $xy^{n-1}$. Since $\operatorname{ord}(x)=3<n$, the prime $n$ does not divide the order of $x$ either. By \Cref{cor. cond. totally NCGp}, the claim follows.
\end{proof}

\begin{example}
    We consider the regular origami $\mathcal{O}\vcentcolon=(A_5, (1,2,3), (1,2,3,4,5))$ given in \Cref{A_n totally NCG} for the prime $n=5$.
    \begin{center}
    \captionsetup{type=figure}
    \input{A5example}
    \caption{Each unlabeled edge is glued to the opposite edge. For clarity, these edges are not labeled.}
    \label{fig::example-A5}
\end{center}
Using the \texttt{GAP} package \cite{OrigamiPackage}, we compute that the Veech group $\operatorname{SL}(\mathcal{O})$ is generated by the matrices
$$S^2, TST^{-1}, T^3,T^{-1}STS^{-1},STST^{-3}S^{-1}$$
and has index $9$ in $\operatorname{SL}(2,\mathbb{Z})$. Representatives of the $\operatorname{SL}(2,\mathbb{Z})$-orbit are given by the following origamis
\begin{align*}
    &\mathcal{O} = (A_5, (1,2,3), (1,2,3,4,5)),\hphantom{\vcentcolon 4_3.,5} \quad\mathcal{O}_2\vcentcolon=(A_5, (2,4)(3,5), (1,2,3,4,5)),\\
    &\mathcal{O}_3\vcentcolon=(A_5, (1,2,4,5,3), (1,2,3,5,4)),\quad~\mathcal{O}_4\vcentcolon=(A_5, (3,5,4), (1,2,3,4,5)),\\
    &\mathcal{O}_5\vcentcolon=(A_5, (1,3,2,5,4), (1,2)(3,4)),\hphantom{5}\quad ~\mathcal{O}_6\vcentcolon=(A_5, (1,2,3,4,5), (1,2,3)),\\
    &\mathcal{O}_7\vcentcolon=(A_5,(1,3,5,4,2), (1,2,3)), \hphantom{4),5}\quad~\mathcal{O}_8\vcentcolon=(A_5, (1,2,3,4,5), (1,2,3,5,4)), \\
    &\mathcal{O}_9\vcentcolon=(A_5, (3,4,5), (1,2,3)).
\end{align*}
\end{example}

\Cref{A_n totally NCG} motivates to examine finite simple groups more generally. Simple groups form an interesting class of $2$-generated groups. The natural question, how the orders of generators for a fixed group can be chosen, has been studied intensively (see e.g. \cite{JLM} for further information). This question suggests to consider $(a,b,c)$-groups.

\begin{definition}
    A finite group generated by two elements $x,y$ with $\operatorname{ord}(x)=a, \operatorname{ord}(y)=b,$ and $\operatorname{ord}(xy)=c$ is called an \textbf{$(a,b,c)$-group}. We call such generators \textbf{$(a,b,c)$-generators}.
\end{definition}

Each $(a,b,c)$-group is a finite quotient of the \textbf{triangle group} $$T_{(a,b,c)}=\langle x,y,z~|~x^a=y^b=z^c=xyz=1\rangle.$$ The following theorem shows that $(a,b,c)$-groups where $a,b,$ and $c$ are chosen pairwise coprime produce regular origamis with a totally non-congruence group as Veech group.

\begin{theorem}\label{(abc)-gps totally NCG}
    Let $a,b,c\in\mathbb{Z}_{\ge0}$ be pairwise coprime and $G$ be an $(a,b,c)$-group with $(a,b,c)$-generators $x,y$. The Veech group of the regular origami $(G,y,x)$  is a totally non-congruence group.
\end{theorem}

\begin{proof}
    We prove that the assumptions of \Cref{thm sp-ws tot NCG} are satisfied for the Veech group of the regular origami $\mathcal{O}=(G,y,x)$. Let $p$ be a prime. Since $a,b,$ and $c$ are pairwise coprime, $p$ divides at most one of the numbers $a,b,$ and $c$. We consider each of the three cases separartely.

    If $p$ is coprime to $b\cdot c$, then consider the matrices $I$ and $S^{-1}T^{-1}=\bigl(
    \begin{smallmatrix}
    0 & 1\\
    -1 & 1
    \end{smallmatrix}
    \bigr).$
    We obtain $I\cdot\mathcal{O}=\mathcal{O}$ and $S^{-1}T^{-1}\cdot\mathcal{O}=(G,xy,y^{-1})$. The inverse moduli of the horizontal cylinders of the regular origamis $\mathcal{O}$ and $(G,xy,y^{-1})$ are $\operatorname{ord}(y)=b$ and $\operatorname{ord}(xy)=c$, respectively. Hence, $T^b$ and $S^{-1}T^{-1}\cdot T^c \cdot TS$ lie in the Veech group $\operatorname{SL}(\mathcal{O})$. Moreover, we obtain $S^{-1}T^{-1}\cdot e_1=\tbinom{0}{-1}\neq j\cdot\tbinom{1}{0}$ modulo $p$ for each $j\in\mathbb{Z}$.

    If $p$ is coprime to $a\cdot c$, then consider the matrices $S^{-1}T^{-1}=\bigl(
    \begin{smallmatrix}
    0 & 1\\
    -1 & 1
    \end{smallmatrix}
    \bigr)$ and $TS^{-1}=\bigl(
    \begin{smallmatrix}
    -1 & 1\\
    -1 & 0
    \end{smallmatrix}
    \bigr).$ We obtain the regular origamis $S^{-1}T^{-1}\cdot\mathcal{O}=(G,xy,y^{-1})$ and $TS^{-1}\cdot\mathcal{O}=(G,x,y^{-1}x^{-1})$. The moduli of the horizontal cylinders of the regular origamis $(G,xy,y^{-1})$ and $(G,x,(xy)^{-1})$ are $\operatorname{ord}(xy)=c$ and $\operatorname{ord}(x)=a$, respectively. Hence, $S^{-1}T^{-1} \cdot T^c\cdot TS$ and $TS^{-1}\cdot T^a \cdot ST^{-1}$ lie in the Veech group $\operatorname{SL}(\mathcal{O})$. Moreover, we obtain for each $j\in\mathbb{Z}$ the inequality $S^{-1}T^{-1}\cdot e_1=\tbinom{0}{-1}\neq j\cdot \tbinom{-1}{-1}= j\cdot TS^{-1}\cdot e_1$ modulo $p$.

    If $p$ is coprime to $a\cdot b$, then consider the matrices $I$ and $S^{-1}=\bigl(
    \begin{smallmatrix}
    0 & 1\\
    -1 & 0
    \end{smallmatrix}
    \bigr)$. We obtain the regular origamis $I\cdot\mathcal{O}=\mathcal{O}$ and $S^{-1}\cdot\mathcal{O}=(G,x,y^{-1})$. The moduli of the horizontal cylinders of the regular origamis $\mathcal{O}$ and $(G,x,y^{-1})$ are $\operatorname{ord}(y)=b$ and $\operatorname{ord}(x)=a$, respectively. Hence, $T^b$ and $S^{-1}\cdot T^a \cdot S$ lie in the Veech group $\operatorname{SL}(\mathcal{O})$. Moreover, we have $S^{-1}\cdot e_1=\tbinom{0}{-1}\neq j\cdot\tbinom{1}{0}$ modulo $p$ for each $j\in\mathbb{Z}$.
\end{proof}

\begin{example}\label{ex hurwitz groups}
    A well-studied family of groups satisfying the assumption in \Cref{(abc)-gps totally NCG} are $(2,3,7)$-groups, which are also called Hurwitz groups. Hurwitz groups are of interest from a geometric point of view because they arise as automorphism groups of compact Riemann surfaces of genus $g>1$ with minimal order $84(g-1)$. The smallest Hurwitz group is the projective linear group $\operatorname{PSL}(2,7)$ and has order $168$. For further information about Hurwitz, see e.g. \cite{hurwitz1} and \cite{hurwitz2}.
\end{example}

\printbibliography[heading=bibintoc]

\end{document}