text:
\section{Introduction \label{sec:Int}}

{The Hwa-Kardar differential stochastic equation\cite{HK,HK1} is a continuous model of self-organized criticality. Unlike equilibrium systems that require ``fine-tuning'' of control parameters (e.g., critical temperature in second order phase transitions) to become critical\cite{Book3}, systems with self-organized criticality arrive at critical state in the natural course of their evolution\cite{Bak3}. Such systems are wide-spread\cite{Col0,Col1,Col2,Col3} with examples ranging from neural systems to online social networks with many other systems in between.}

{The Hwa-Kardar model is expected to describe
the infrared (IR, large-scale, long-time) asymptotic behaviour of the system with self-organized criticality, i.e. to yield the power laws that characterize the behaviour.}

{The system is an anisotropic ``running'' sandpile with an average flat profile that has a constant tilt along the preferred direction. The new sand keeps entering the system resulting in avalanches that cause other sand to exit the system. Eventually, the system arrives at a steady state with a self-similar profile that is anisotropic (due to the tilt).}

{The height of the profile is described by the scalar field $h(x)=h(t,{\bf x})$ where $t$ and a $d$-dimensional ${\bf x}$ stand for the time-space coordinates. Note that $h(x)$ is measured from the average tilted profile, i.e. it is a deviation of the profile height from its average value. The tilt corresponds to the preferred direction defined by the unit constant vector $ {\bf n}$. Any vector ${\bf q}$ can be decomposed as ${\bf q} = {\bf q}_{\perp} + {\bf n}\, q_{\parallel}$ where
$({\bf q}_{\perp} \cdot \, {\bf n}) =0$. The spatial derivative ${\bf \partial}=\partial/ \partial {x_{i}}$, $i=1,\dots, d$, is then ``split'' into two derivatives: a derivative ${\bf \partial_{\perp}}=\partial/ \partial {x_{i}}$, $i=1,\dots, (d-1)$, in the subspace orthogonal to the vector $ {\bf n}$, and a one-dimensional derivative
$\partial_{\parallel} = ({\bf n} \cdot {\bf \partial})$.}

The stochastic differential equation that describes the evolution of the system is
\begin{equation}
\partial_{t} h= \nu_{\perp 0}\, {\bf \partial}_{\perp}^{2} h + \nu_{\parallel 0}\,
\partial_{\parallel}^{2} h -
\partial_{\parallel} h^{2}/2 + f.
\label{eq1}
\end{equation}
Here it is assumed that $\partial_{t} = \partial/ \partial {t}$ and ${\bf \partial}_{\perp}^{2}=({\bf \partial}_{\perp}\cdot{\bf \partial}_{\perp})$. The factors $\nu_{\parallel 0}$ and
$\nu_{\perp 0}$ stand for the diffusivity coefficients, and $f(x)$ is a Gaussian random noise with a zero mean and the correlation function
\begin{equation}
\langle f(x)f(x') \rangle = D_{0}\,
\delta(t-t')\, \delta^{(d)}({\bf x}-{\bf x}'), \quad D_0> 0.
\label{forceD}
\end{equation}

{Critical behaviour of a system is sensitive to various external disturbances including, in particular, motion of the environment\cite{Onuki,Nelson,Satten}. Thus, it is desirable to consider the effect of the motion on a system with self-organized criticality. }

{The coupling with the velocity field ${\bf v}(x)$ is achieved through a ``minimal substitution''
\mbox{$\partial_{t} \to \nabla_{t} = \partial_{t} + ({\bf v}\cdot {\bf \partial})$}.}

{Let us ascribe the following statistic to the velocity field ${\bf v}(x)$:
\begin{equation}
\begin{aligned}
\langle v_{i} (t, {\bf x}) v_{j}(t',{\bf x}')\rangle =  \delta(t-t')
D_{ij}({\bf x}-{\bf x}'),
\\
D_{ij}({\bf r}) = B_{0}
\int_{k>m} \frac{d{\bf k}}{(2\pi)^{d}} \frac{1}{k^{d+\xi}}
 P_{ij}({\bf k})
\exp ({\rm i} {\bf k}\cdot{\bf r}), \quad B_{0}>0.
\label{white}
\end{aligned}
\end{equation}
{This is the Kazantsev-Kraichnan ``rapid-change'' velocity ensemble\cite{FGV}.
Here $k\equiv |{\bf k}|$ stands for the wave number while $P_{ij}({\bf k})$ is
the transverse
projector defined in a standard way: $P_{ij}({\bf k}) = \delta_{ij} - k_i k_j / k^2$}. The velocity field is, thus, incompressible, i.e.
$({\bf  \partial}\cdot {\bf v}) = 0$. The IR regularization is provided by the sharp cutoff $k>m$. The parameter $\xi$ is non-negative and smaller than $2$.}

{In this paper, we apply field-theoretic renormalization group (RG) analysis to problem~(\ref{eq1})~-- (\ref{white}) and find out that it is the coupling of the anisotropic equation~(\ref{eq1}) and the isotropic velocity ensemble~(\ref{white}) that causes emergence of ``nonconventional'' scaling behaviour. In RG analysis, stochastic problem~(\ref{eq1})~-- (\ref{white}) is substituted with a field theory. Critical exponents from the power laws that describe the IR asymptotic behaviour of the system can be found from RG equations for the field theory. In particular, each IR attractive fixed point of
the RG equations is connected to a regime of asymptotic behaviour, i.e. to a universality class. It was shown in [\refcite{WeU}] that the system's critical behaviour is divided into four universality classes. Two of them correspond to ordinary isotropic scaling in which infrared irrelevant parameters are kept fixed and the coordinates $x_{\parallel}$ and ${\bf x_{\bot}}$ are not scaled independently. The other two regimes, however, involve nonconventional or anisotropic scaling. The aim of this paper is to explore these two regimes.}

\section{Renormalization and scaling of pure Hwa-Kardar model\label{sec:QFT1}}

{Let us begin by deriving scaling of the Hwa-Kardar equation without the turbulent motion of the environment. This ``pure'' scaling was first analysed by RG in its Wilsonion form in [\refcite{HK,HK1}]. We reproduce those results here (expressed in terms of the field theoretic RG) so that we can compare them with the scaling regimes of the model with turbulent advection.}

{Instead of the
stochastic problem~(\ref{eq1})~-- (\ref{forceD}), one can consider the equivalent field theory\cite{Book3} with the doubled set of fields $\Phi=\{h,h'\}$ and the action}
\begin{equation}
    S(\Phi) = \frac{1}{2}h'D_0\,h' + h'(-\partial_th + \nu_{\parallel0}\,\partial^2_{\parallel}h + \nu_{\perp0}\,{\bf \partial}^2_{\perp}h - \frac{1}{2}\,\partial_{\parallel}h^2).
    \label{action_pure}
\end{equation}
{The integrations over the space-time coordinates $x=\{t,{\bf x}\}$ and the summations over the vector indices are implied throughout.}

{As model (\ref{eq1}), (\ref{forceD}) is anisotropic, it involves two independent spatial scales $L_{\parallel}$ and $L_{\perp}$ instead of a single scale $L$. Generally, a quantity's canonical dimension in dynamical models is described by the momentum dimension $d^{k}_{F}$ and the frequency dimension $d^{\omega}_{F}$ related to the spatial scale $L$ and the temporal scale $T$. Since there are two independent spatial scales, an arbitrary quantity $F$ is described by three canonical dimensions:}
\begin{equation}
    [F] \sim [T]^{-d^{\omega}_{F}}[L_{\parallel}]^{-d^{\parallel}_{F}}[L_{\perp}]^{-d^{\perp}_{F}}.
\end{equation}
{Table~\ref{canon_dim_pure} contains the canonical dimensions of all fields and parameters for the theory~(\ref{action_pure}) obtained from the condition that each term of the action~\eqref{action_pure} be dimensionless. We denoted a total canonical dimension $d_{F}$ as $d_{F} = d_{F}^{k} + 2d_{F}^{\omega}$ and a total momentum dimension $d_{F}^{k}$ as $d_{F}^{k} = d_F^{\parallel} + d_F^{\perp}$. The coupling constant $g_0$ is defined as} $D_0 = g_0\nu_{\perp0}^{3/2}\nu_{\parallel0}^{3/2}$.

{We can see from Table~\ref{canon_dim_pure} that  $g_0\sim\Lambda^{\varepsilon}$ where $\Lambda$ is a characteristic ultraviolet momentum scale and $\varepsilon= 4-d$, therefore, the logarithmic dimension is $d=4$. Let us also define the parameter $u_0 = \nu_{\parallel0}/\nu_{\perp0}$ that we will need later. Note that it possesses nontrivial momentum canonical dimensions $d_u^\parallel$ and $d_u^\perp$ while its total canonical dimension $d_{u}$ is zero.}

\begin{table}[t]
\tbl{Canonical dimensions for the theory~(\ref{action_pure}); $\varepsilon=4-d$.}
{\begin{tabular}{cccccccccc} \toprule
$\pmb F$ & $\pmb h'$ & $\pmb h$ & $\pmb D_0$ & $\pmb \nu_{\pmb \parallel\pmb 0}$ &$ \pmb \nu_{\pmb \perp\pmb 0}$& \boldmath {$ u_0$}&\boldmath {$g_0$}&$\pmb g$&$\pmb \mu$\pmb ,$\pmb m$\pmb ,$\pmb \Lambda$ \\  \colrule
$d_F^{\omega}$&$-1$&$1$&$3$&$1$&$1$&$0$&$0$&$0$&$0$\\
$d_F^{\parallel}$&$2$&$-1$&$-3$&$-2$&$0$&$-2$&$0$&$0$&$0$\\
$d_F^{\perp}$&$d-1$&$0$&$1-d$&$0$&$-2$&$2$&$\varepsilon$&$0$&$1$\\
$d_F^{k}$&$d+1$&$-1$&$-d-2$&$-2$&$-2$&$0$&$\varepsilon$&$0$&$1$\\
$d_F$&$d-1$&$1$&$4-d$&$0$&$0$&$0$&$\varepsilon$&$0$&$1$\\ \botrule
\end{tabular} \label{canon_dim_pure}}
\end{table}

{Canonical dimensions analysis augmented with symmetry analysis reveals that the model~(\ref{action_pure}) is multiplicatively renormalizable. It involves the sole non-trivial renormalization constant $Z_{\nu_{\parallel}} $. The only parameters that require renormalization are $\nu_{\parallel0}$ and $g_0$. They are related to their renormalized~counterparts in the following way ($\mu$ is the renormalization mass):}
\begin{equation}
\nu_{\parallel0} = \nu_{\parallel} Z_{\nu_{\parallel}} , \quad
g_0= g\mu^{\varepsilon} Z_g.
\label{REN1}
\end{equation}

{The canonical differential equations for a renormalized Green function $G^R=\langle \Phi \dots \Phi \rangle$ read}
\begin{equation}
\begin{aligned}
    \left(\sum_i d_i^{\omega} {\cal D}_i - d_{G}^{\omega} \right)G^R &= 0, \\
    \left(\sum_i d_i^{\perp} {\cal D}_i - d_{G}^{\perp} \right)G^R &= 0,\quad
    \left(\sum_i d_i^{\parallel} {\cal D}_i - d_{G}^{\parallel} \right)G^R &= 0.
\label{ScInvX}
\end{aligned}
\end{equation}
{The index $i$ enumerates the arguments of $G^R$ which are $\omega$, $k_{\perp}$, $k_{\parallel}$, $\mu$, $\nu_{\perp}$, and $\nu_{\parallel}$ (or $u$ depending on our choice of parameters). The operator ${\cal D}_a$ is ${\cal D}_a=a\partial_a$ for any $a$.}

{To calculate the critical exponents (or dimensions) that stand in scaling power laws, we need to supplement equations (\ref{ScInvX}) with the differential RG equation:}
\begin{equation}
{\left(
{\cal D}_{\mu} + \beta_g\partial_g - \gamma_{\nu_\parallel}{\cal D}_{\nu_\parallel}
- \gamma_{\nu_\perp}{\cal D}_{\nu_\perp} - \gamma_G \right) G^{R} = 0.}
\label{RGeqX}
\end{equation}
{Here $\gamma$ and $\beta$ are RG functions
(anomalous dimensions and a  $\beta$-function respectively)}
{for the parameters of the system and coupling constant $g$. Note that quantities that are not renormalized have vanishing anomalous dimensions, i.e. $\gamma_{\nu_{\perp}}=0$ in the model (\ref{action_pure}), but we keep this term for the future.}

{Universality classes of asymptotic (critical) behaviour are determined by fixed points of the RG equations, i.e. by the zeroes of the $\beta$ functions, $\beta(g^*)=0$. Real parts of all the eigenvalues $\lambda_i$ of the matrix $\Omega_{ij}=\partial_{g_i}\beta_{g_j}$ (here $g =\{g_i\}$ is a full set of the charges) must be positive for a fixed point to be IR attractive.}

The substitution $g\to g^*$ and, hence, $\gamma_F \to\gamma_F ^*$ turns equation~\eqref{RGeqX} into equation with constant coefficients, i.e. into equation of the same type as equations~\eqref{ScInvX}.
Each of these equations describes a certain independent scaling behaviour of the function $G^R$, in which some of its variables are scaled and some are kept fixed. A~parameter is scaled if the corresponding derivative enters the differential operator; otherwise the parameter is fixed.

We are interested in the critical scaling behaviour where the frequencies and momenta (or, equivalently, times and coordinates) are scaled, while the IR irrelevant parameters \mbox{(namely, $\mu$, $\nu_{\perp}$
and $\nu_{\parallel}$)} are kept fixed.
Thus, we combine all these equations to exclude
the derivatives with respect to all the IR irrelevant parameters
and arrive at the critical scaling equation for a given fixed point:
\begin{equation}
    \left({\cal D}_{k_\perp} +{\cal D}_{k_\parallel} \Delta_\parallel+ \Delta_\omega {\cal D}_\omega
        - \Delta_G  \right)G^{R} = 0
    \label{HK-scalingXZ}
\end{equation}
with $\Delta_{\parallel} = 1 + \gamma_{\nu_{\parallel}}^{*}/2$ and $\Delta_\omega=2-\gamma_{\nu_\perp}^*$.

As we will see below, such an exclusion is not always possible: it requires some balance between the numbers of IR relevant and IR irrelevant parameters and the number of independent scaling equations.

The model has two fixed points. The Gaussian (free) fixed point with $g^{*}=0$ is IR attractive for $\varepsilon<0$; the corresponding
critical dimensions coincide with canonical ones.
The nontrivial fixed point with $g^*=32\varepsilon/9+O(\varepsilon^2)$ is IR attractive for $\varepsilon>0$. The corresponding canonical dimensions read
\begin{equation}
        \Delta_h = 1-\varepsilon/3, \quad \Delta_{h'} = 3 - \varepsilon/3, \quad
        \Delta_{\omega}=2, \quad   \Delta_{\parallel} = 1+ \varepsilon/3.
        \label{HKcrit}
\end{equation}

\section{Renormalization and fixed points of the model with turbulent advection\label{sec:QFT2}}

As was stated above, inclusion of turbulent advection in the Hwa-Kardar model is achieved by replacing $\partial_t$ with $\nabla_t$ in the action (\ref{action_pure}) and adding Gaussian averaging over the velocity ensemble with the correlation function~(\ref{white}). This leads to the following action functional:
\begin{equation}
\begin{aligned}
    S(\Phi) = \frac{1}{2}h'D_0\,h' + h'(-\nabla_th + \nu_{\parallel0}\partial^2_{\parallel}h + \nu_{\perp0}\partial^2_{\perp}h - \frac{1}{2}\partial_{\parallel}h^2) + S_{{\bf v}},
    \label{action_ve}
\end{aligned}
\end{equation}
\begin{equation}
    S_{{\bf v}} = -\frac{1}{2}\int dt \int d{\bf x} \int d{\bf x}' v_i(t,{\bf x})D^{-1}_{ij}({\bf x}-{\bf x}')v_j(t,{\bf x}').
\label{action_Sv}
\end{equation}
Here $D^{-1}_{ij}({\bf x} - {\bf x'})$ is the kernel of the inverse operator $D^{-1}_{ij}$ for the integral operator $D_{ij}$ from~(\ref{white}) as $S_{{\bf v}}$ describes  Gaussian averaging over the field ${\bf v}$ with the correlation function~(\ref{white}).

Since that velocity ensemble is isotropic, it is no longer possible to define two independent spatial scales in this model. Thus,
\begin{equation}
    [F] \sim [T]^{-d^{\omega}_{F}}[L]^{-d^{k}_{F}}
    \label{t5}
\end{equation}
and $d_{F} = d_{F}^{k} + 2d_{F}^{\omega}$.
Canonical dimensions of the fields and parameters of the model are presented in Table~\ref{canon_dim_ve}.
We note that for the quantities that
are present in the pure Hwa-Kardar model with the action
(\ref{action_pure}) these dimensions coincide with their counterparts in Table~\ref{canon_dim_pure}.

\begin{table}[t]
\tbl{Canonical dimensions for the theory (\ref{action_Sv}); $\varepsilon=4-d$.}
{\begin{tabular}{ccccccccccccc}
\toprule
\boldmath {$F$} & \boldmath {$h'$} & \boldmath {$h$} & \boldmath {$D_0$} & \boldmath {$\nu_{\parallel0}$} &\boldmath {$ \nu_{\perp0}$}&\boldmath {$g_0$}&\boldmath {$v$}&\boldmath {$B_0$}&\boldmath {$x_0$}&\boldmath {$u_0$},\boldmath {$u$}&\boldmath {$g$},\boldmath {$x$}&\boldmath {$\mu$}, \boldmath {$m$}, \boldmath {$\Lambda$} \\ \colrule
$d_F^{\omega}$&$-1$&$1$&$3$&$1$&$1$&$0$&$1$&$0$&$0$&$0$&$0$&$0$\\
$d_F^{k}$&$d+1$&$-1$&$-2-d$&$-2$&$-2$&$\varepsilon$&$-1$&$\xi-2$&$\xi$&$0$&$0$&$1$\\
$d_F$&$d-1$&$1$&$4-d$&$0$&$0$&$\varepsilon$&$1$&$\xi$&
$\xi$&$0$&$0$&$1$\\\botrule
\end{tabular}
\label{canon_dim_ve}}
\end{table}

In contrast to the pure Hwa-Kardar model and its modifications with anisotropic velocity ensembles\cite{AK1,Serov,SerovEP} (where the diffusivity coefficients $\nu_{\parallel0}$ and $\nu_{\perp0}$ have different momentum dimensions $d^\parallel$ and $d^\perp$), only total momentum dimension $d^k$ can be defined for the theory~\eqref{action_ve}.
From Table~\ref{canon_dim_ve} it follows that
the ratio $u_0=\nu_{\parallel0}/\nu_{\perp0}$ is completely dimensionless, that is, dimensionless with respect to the frequency and momentum dimensions separately.
Therefore, according to the general rules, $u_0$~should be treated as an additional charge.
The third coupling constant is $x_0\sim\Lambda^{\xi}$ related
 to the amplitude $B_0$ of the correlation function~(\ref{white}) as
$B_0 = x_0\nu_{\perp0}$. The theory is logarithmic when $\varepsilon=0$ and $\xi=0$.

The theory is renormalized by introducing the renormalization constants $Z_i$:
\begin{equation}
\begin{aligned}
    \nu_{\parallel0} = \nu_{\parallel} Z_{\nu_{\parallel}} , \quad &\nu_{\perp0} = \nu_{\perp} Z_{\nu_{\perp}}, \quad
    g_0=Z_g g\mu^{\varepsilon},\quad x_0=&Z_x x\mu^{\xi},\quad  u_0 = Z_u u.
    \label{renorm_param}
\end{aligned}
\end{equation}
The anomalous dimensions calculated from the Feynman graphs have the following one-loop expressions (we omit details of calculations\cite{AKL} for brevity):
\begin{equation}
    \gamma_{\nu_{\parallel}} = \frac{3}{8}\frac{x}{u} + \frac{3}{16}g, \quad \gamma_{\nu_{\perp}} = \frac{3}{8}x.
    \label{gamma1_gamma2}
\end{equation}
The anomalous dimensions for the coupling constants read
\begin{equation}
    \gamma_{g} = - \frac{3}{2}\left(\gamma_{\nu_{\parallel}} + \gamma_{\nu_{\perp}}\right), \quad \gamma_x = -\gamma_{\nu_{\perp}}, \quad \gamma_u = \gamma_{\nu_{\parallel}} - \gamma_{\nu_{\perp}}.
    \label{an-2}
\end{equation}
Similarly to the case of the pure Hwa-Kardar model (see Section \ref{sec:QFT1}), the~anomalous dimensions for the fields $h$, $h'$, and ${\bf v}$ vanish due to the absence of their renormalization: \mbox{$\gamma_h=\gamma_{h'}=\gamma_v=0$}.

From explicit results for the anomalous dimensions~(\ref{gamma1_gamma2}), it follows that the  one-loop expressions for the $\beta$ functions have the form
\begin{equation}
\begin{aligned}
\beta_g &= g\left(-\varepsilon +\frac{9}{32}g+\frac{9}{16}\frac{x}{u}+\frac{9}{16}x\right), \quad
\beta_x =x\left(-\xi+\frac{3}{8}x\right), \\
\beta_u &= u\left(-\frac{3}{16}g-\frac{3}{8}\frac{x}{u}+\frac{3}{8}x\right).
\end{aligned}
\label{beta-expr1}
\end{equation}
Analysis of the system~\eqref{beta-expr1}
reveals two possible IR attractive fixed points: the Gaussian point FP1 with the coordinates $g^*=0$, $x^*=0$ and arbitrary
$u^*$, and the fixed point FP2 with the coordinates $g^*=0$, $x^*=8\xi/3$, $u^*=1$ which corresponds to the regime of simple turbulent advection
(the nonlinearity of the Hwa-Kardar equation is irrelevant in the sense of Wilson).

The point FP1 is IR attractive for $\varepsilon<0$, $\xi<0$, whilst
the point FP2 is IR attractive for $\xi>\varepsilon/3$, $\xi>0$. This is the full set of fixed points with finite and non-zero value of $u^*$; there are no fixed points with $g^*\neq 0$ in the set, i.e. the Hwa-Kardar universality class is not realized for such values of $u$.

However,  the system~(\ref{beta-expr1}) may lose possible solutions with $u^*=0$ or $1/u^*=0$.
In order to explore those exceptional values, we have to pass to new variables: $w=x/u$ (to study the case
$u^*=0$) and to $\alpha=1/u$ (to study the case
$u^*\to\infty$).

For the first case we obtain no IR attractive fixed points.
In the second case we have
two more fixed points:
the point FP3 with the coordinates $g^*=32\varepsilon/9$, $x^*=0$, $\alpha^*=0$ and the point FP4 with  $g^*=32\varepsilon/9-16\xi/3$, $x^*=8\xi/3$, $\alpha^*=0$. The point FP3 corresponds to the regime of critical behaviour where only the nonlinearity of the Hwa-Kardar equation is relevant, and the point FP4 corresponds to the regime where both the nonlinearity
and the turbulent advection are relevant. The point FP3 is IR attractive if $\varepsilon>0$, $\xi<0$, and the point FP4 is IR attractive if $\xi<\varepsilon/3$, $\xi>0$.

\vspace{-.3em}
\section{Scaling Regimes and Critical Dimensions in the Model with Turbulent Advection\label{sec:scal}}

Unlike model (\ref{action_Sv}), the original Hwa-Kardar model\cite{HK,HK1} allows to introduce two independent momentum dimensions.
Thus, it involves three equations~(\ref{ScInvX}) related to the canonical scale invariance. On the other hand, the model with turbulent advection contains only two such equations.
Combining them with differential RG equation
\begin{equation}
\left(
{\cal D}_{\mu} + \beta_g\partial_g + \beta_x\partial_x+ \beta_u\partial_u - \gamma_{\nu_\perp}{\cal D}_{\nu_\perp} -\gamma_G\right) G^{R} = 0
\label{RGeqXX}
\end{equation}
taken at a fixed point, we arrive at
\begin{equation}
    \left({\cal D}_{k_\parallel}+{\cal D}_{k_\perp} +\Delta_\omega{\cal D}_\omega-d^k_G -\Delta_\omega d^\omega_G-\gamma^*_G\right)G^{R} = 0,
    \label{G-scaling}
\end{equation}
where $\Delta_\omega=2-\gamma_{\nu_\perp}^*$.
This is the critical scaling equation; note that it does not involve derivatives over $\mu$ and $\nu_\perp$.
The critical dimension $\Delta_F$ of a quantity $F$ now reads
\vspace{-.9em}
\begin{equation}
    \Delta_F=d^k_F +\Delta_\omega d^\omega_F+\gamma^*_F.
    \label{crete}
\end{equation}
For the fixed points FP2 critical dimensions are
\begin{equation}
    \Delta_h=1-\xi, \quad \Delta_v=1-\xi, \quad \Delta_{h'}=3-\varepsilon+\xi, \quad \Delta_{\omega}=2-\xi.
\end{equation}
These results are perturbatively exact due to the fact that $\gamma_{\nu_\perp}^*$ is known exactly while $\gamma_{h}=\gamma_{h'}=\gamma_{v}=0$.

The fixed point FP1 is Gaussian so the corresponding critical dimensions coincide with the canonical ones.

Surprisingly, the fixed points FP3 and FP4 are related to scaling that is majorly different from the scaling described by the points FP1 and FP2. The difference arises when one sets $\beta_{\left\{g_i^*\right\}}=0$ in the RG equation (\ref{RGeqXX}). The Green functions have well-defined finite limits at $g\to0$ and $x\to 0$, thus, the $\beta$-functions at the points FP1 and FP2 can be set to zero without the need for any further analysis. However, such a straightforward substitution cannot be performed for $\beta_\alpha$ at the fixed points with $\alpha^*=0$ (or alternatively at the fixed points with $u\to\infty$). Instead, $\beta_\alpha$ should be replaced with the first nontrivial order of its expansion~(\ref{RGeqXX}).
Thus, the critical scaling equation at the fixed points FP3 and FP4 reads
\begin{equation}
    \left({\cal D}_{k_\parallel}+{\cal D}_{k_\perp} +\Delta_\omega{\cal D}_\omega-
    \lambda^*{\cal D}_\alpha -d^k_G -\Delta_\omega d^\omega_G-\gamma^*_G\right)G^{R} = 0,
    \label{34-scaling}
\end{equation}
where $\lambda=\partial \beta_\alpha/\partial\alpha$ taken at $\alpha=0$, $\Delta_\omega=2-\gamma_{\nu_\perp}^*$, and $\lambda^*=\lambda(g^*,x^*)$.

The point FP3 corresponds to the regime where only the nonlinearity of the Hwa-Kardar equation is relevant; its scaling equation takes the form
\begin{equation}
    \left({\cal D}_{k_\parallel}+{\cal D}_{k_\perp} +2{\cal D}_\omega-
    \frac{2\varepsilon}{3}
    \, {\cal D}_\alpha -d^k_G -2d^\omega_G-\gamma^*_G\right)G^{R} = 0.
    \label{3-scaling}
\end{equation}
This equation corresponds to scaling where the momenta $k_\parallel$, ${\bf k}_{\bot}$, the~frequency~$\omega$, and the ratio $\alpha$ are scaled. However, in this special case where only the nonlinearity is relevant, the model~\eqref{action_ve} coincides with the pure Hwa-Kardar model~\eqref{action_pure}. It means that additional canonical symmetry arises and turns $d^\parallel$ and~$d^\perp$ into independent dimensions, see~Equations~\eqref{ScInvX}.

The general solution of the system of equations that includes two canonical invariance equations, the
RG equation~(\ref{RGeqXX}) taken at a fixed point, and~homogeneous counterpart of equation~\eqref{3-scaling} is an arbitrary function of three independent variables chosen here for definiteness as
\begin{equation}
    z_1=\frac{\omega}{\nu_\perp k_\perp^2}, \quad
    z_2=\frac{k_\parallel}{k_\perp}, \quad \text{and}
    \quad z_3=  \alpha \left(\frac{k_\perp}{\mu}\right) ^{2\varepsilon/3}.
    \label{sol1}
\end{equation}

Additional symmetry requires the variables $z_1$, $z_2$, and $z_3$ to be dimensionless with respect to Table~\ref{canon_dim_pure}.
While the variables $z_2$ and $z_3$ do not satisfy this requirement, it is possible to introduce a new variable $z_0=z_2z_3^{-1/2}$ with needed canonical dimensions that serves as the second solution (along with $z_1$) of the homogeneous part of the equation~\eqref{HK-scalingXZ}:
\begin{equation}
    z_0=
    \frac{k_\parallel }{k_\perp^{\Delta_\parallel} \alpha^{\wp}}\, \mu^{\varepsilon/3}
    \quad \text{with}
     \quad \wp=\frac{3}{2\varepsilon}\left(\Delta_\parallel-1\right),
    \label{sol2}
\end{equation}
where $\Delta_\parallel=1+\varepsilon/3$ is in agreement with expressions~\eqref{HKcrit}. This implies that the fixed point FP3 allows for the scaling where the coordinates $x_\parallel$ and ${\bf x}_{\bot}$ (or momenta $k_\parallel$ and ${\bf k}_{\bot}$)
are scaled simultaneously with nontrivial $\Delta_\parallel\neq1$, while all the IR irrelevant parameters (including $\alpha$) are kept fixed. Thus, the results derived for the pure Hwa-Kardar model~\eqref{action_pure} are reproduced.

The fixed point FP4 corresponds to the regime where both the nonlinearity of the Hwa-Kardar equation and the turbulent advection are relevant; its equation of critical scaling~\eqref{34-scaling} reads
\begin{equation}
    \left({\cal D}_{k_\parallel}+{\cal D}_{k_\perp} +\Delta_\omega{\cal D}_\omega-
    \left(\frac{2}{3}\varepsilon-2\xi\right){\cal D}_\alpha -d^k_G -\Delta_\omega d^\omega_G-\gamma^*_G\right)G^{R} = 0,
    \label{4-scaling}
\end{equation}
where $\Delta_\omega=2-\xi$. Three combinations $\hat{z}_1$, $\hat{z}_2$, and $\hat{z}_3$ are possible solutions of its homogeneous part:
\begin{equation}
    \hat{z}_1= \frac{\omega}{\nu_\perp k_\perp^{2}}\left(\frac{k_\perp}{\mu}\right)^\xi, \quad
    \hat{z}_2= \frac{k_\parallel}{k_\perp}, \quad
    \hat{z}_3= \alpha \left(\frac{k_\perp}{\mu}\right) ^{2\varepsilon/3-2\xi}.
    \label{FP4-an}
\end{equation}
The variables~\eqref{FP4-an} describe generalized critical scaling where
$\alpha$ (i.e. the ratio of $\nu_\perp$ and~$\nu_\parallel$) is also scaled.

This resembles various modified similarity hypotheses for systems with  different characteristic scales or different scaling laws\cite{Stell,Stell2,Stell3,Stell4}.

Of course, it is not impossible that the Green functions  have some peculiar dependence on the parameters (\ref{FP4-an}) so that the ordinary critical scaling at fixed $u$, with some nontrivial scaling dimensions, takes place. However, as far as we can see, there is no anticipated reason for this situation to occur.

\section{Conclusion \label{sect:Con}}

We studied nonconventional scaling behavior of the self-organized critical system in an insotropic turbulent environment. The system was modelled by the anisotropic Hwa-Kardar equation~(\ref{eq1})~-- (\ref{forceD}) while the turbulent advection was described by Kazantsev-Kraichnan ``rapid-change'' ensemble~(\ref{white}).

We analyzed a field theory equivalent to the stochastic problem and established that the theory is multiplicatively renormalizable. The coordinates of the four fixed points of RG equations were calculated in the one-loop approximation or found exactly.

Among four universality classes of critical behaviour, two classes turned out to correspond to nonconventional scaling behavior. It was shown that a kind of dimensional transmutation takes place in one of those regimes, precisely, in the one where the nonlinearity of the Hwa-Kardar equation is relevant. The transmutation results in the ratio $u$ of the two diffusivity coefficients $\nu_\parallel$ and $\nu_\perp$ acquiring a nontrivial canonical dimension. Due to this, the new canonical symmetry arises in the model~\eqref{action_ve} and allows canonical dimensions $d^\parallel$ and $d^\perp$ to be scaled independently. Thus, the scaling behaviour involves simultaneous scaling of the coordinates $x_\parallel$ and ${\bf x}_{\bot}$ with a nontrivial relative exponent $\Delta_\parallel\neq1$. The~regime is in agreement with the predictions of the pure Hwa-Kardar model analysis~\eqref{action_pure}.

In the regime where both the advection and nonlinearity are relevant, the self-similar behaviour necessarily involves some dilation of the ratio $u$.
The ordinary critical behaviour with fixed definite critical dimensions can be established only if the correlation functions of the model have some special dependence on $u$. Otherwise, the IR behaviour will remind
of various generalized self-similarity hypotheses, like the
weak scaling due to George Stell\cite{Stell,Stell2,Stell3}
or the parametric scaling in the spirit of Michael Fisher\cite{Stell4}
for systems with different characteristic scales and modified scaling laws.

\section*{Acknowledgments}
The reported study was funded by the Russian Foundation for Basic Research, project number~20-32-70139.
The work by N.V.A.~and P.I.K. was also supported by the Theoretical Physics and Mathematics Advancement Foundation ``BASIS.''

\pagebreak


text:
\section{Introduction}

In this paper we introduce push-out space  for an immersion
 $f:M^m\longrightarrow \Bbb R^{m+k}$, where $M$ is a smooth connected
 $m$-dimensional manifold without boundary. To do this, we give  some examples
in 3-dimensional Euclidean space, $\Bbb R^3$, infact, in these
examples we calculate normal holonomy group and push-out space
geometrically. We consider the case when $\mathcal{H}ol(f)$ is
non-trivial. This extends the work of Carter and Senturk \cite{2},
who obtained results about the case when $\mathcal{H}ol(f)$ is
trivial. In these examples we show that some of the properties of
push-out space which they obtained is not true for the case when
$\mathcal{H}ol(f)$ is non-trivial.

\section{Basic definitions}

\begin{defn}[\cite{2}]\label{de1} Let $f:M^m\longrightarrow \Bbb R^{m+k}$ be a
smooth immersion where $M$ is a smooth connected $m$-dimensional
manifold without boundary. The total space of the normal bundle of
$f$ is defined by
  $$N(f)=\{(p,x)\in M\times\Bbb R^{m+k}:<x,v>=0\quad  \forall v\in
  f_{*}T_{p}(M)\}$$
The endpoint map $\eta:N(f)\longrightarrow\Bbb R^{m+k}$ is defined
by $\eta(p,x)=f(p)+x$ and, the set of singular points of $\eta$ is
subset $\Sigma(f)\subset N(f)$ called the set of critical normals of
$f$ and the set of focal points of $\eta$ is a subset
$\eta({\Sigma(f)})\subset \Bbb R^{m+k}$.\\
For $p\in M$, we put $N_p(f)=\{x:(p,x)\in N(f)\}$ and
${\Sigma}_{p}(f)=\{x:(p,x)\in \Sigma(f)\}$ respectively, normal
space at $p$ and the set can be thought of as focal points with base
$p$.
\end{defn}
\begin{defn}[\cite{1}]\label{de2}
For $p_0\in M$ and  $p\in M$ and path $\gamma:[0,1]\longrightarrow
M$ from $p_0$ to $p$ define
$\varphi_{p,\gamma}:N_{p_0}(f)\longrightarrow N_p(f)$ by parallel
transport along $\gamma$. The $\varphi_{p,\gamma}$'s are isometries.
The normal holonomy group on $N_{p_0}(f)$, is
$$\mathcal{H}ol(f)=\{\varphi_{p_0,\gamma}:\gamma:[0,1]\longrightarrow M,\quad \gamma(0)=\gamma(1)=p_0\}$$
If the closed path $\gamma$ at $p_0$ is homotopically trivial then
$\varphi_{p_0,\gamma}$ is an element of the restricted normal
holonomy group ${\mathcal{H}ol}_0(f)$.
\end{defn}
\begin{defn}[\cite{3}]\label{de3}
For a fix $p_0\in M$ the push-out space for an immersion
$f:M^m\longrightarrow \Bbb R^{m+k}$ is defined by
$$\Omega(f)=\{x\in N_{p_0}(f):\forall p\in M,\forall \gamma\mbox{ s.t. }\gamma(0)=p_0,\gamma(1)=p \mbox{ then } \varphi_{p,\gamma}x\notin\Sigma_p(f) \}$$
(i.e. $\forall p\in M$, $f(p)+\varphi_{p,\gamma}(x)$ is not a focal
point with base $p$ when $x$ belongs to $\Omega(f)$). Therefore
$\Omega(f)$ is the set of normals at $p_0$, where transported
parallely along all curves, do not meet focal points. So $\Omega(f)$
is invariant under the action of $\mathcal{H}ol(f)$.
\end{defn}

\begin{defn}[\cite{4}]\label{de4}
Let $B\subset N(f)$ be a smooth subbundle with type fiber S where\\
1) S is a smooth submanifold of $\Bbb R^k$\\
2) $B\cap\Sigma(f)=\emptyset$\\
3) B is invariant under parallel transport (along any curve in M).
Then B is a smooth manifold and $g\equiv\eta|_B:B\longrightarrow\Bbb
R^{m+k}$ is a smooth immersion called a partial tube about f.
\end{defn}

\begin{thm}[\cite{2}]\label{2.5}
Let $ \mathcal{H}ol(f) $ is trivial and $ M $ be a compact manifold. Then each path-connected component of
$\Omega(f)$ is open in $\Bbb R^k$.
\end{thm}

\begin{thm}[\cite{2}]\label{2.6}
Let $ \mathcal{H}ol(f) $ is trivial then
Each path-connected component of $\Omega(f)$ is convex.
\end{thm}

\begin{rem}
In Example \ref{ex:3.2}, if $\frac{\alpha}{\pi}$ is irrational then
$\Omega(\bar{f})$ is not open in $\Bbb R^2$ but $M=\Bbb S^1$ is
compact. Also, in Example \ref{ex:3.5}, $\Omega({f})=\{O\}$ hence
$\Omega({f})$ is closed in $\Bbb R^2$ but $\mathcal{H}ol(f)$ is
trivial. This shows that Theorem \ref{2.5} is false when M is not
compact or $\mathcal{H}ol(f)$ is non-trivial.
\end{rem}

\begin{rem}
In Example \ref{ex:3.6} one of path-connected components of
$\Omega(f)$, which is the complement space of cone and two other
components in $\Bbb R^3$, is not convex. This shows that Theorem
\ref{2.6} is false when $\mathcal{H}ol(f)$ is non-trivial.
\end{rem}
we conclude that the properties of push-out space that are proved in
the case $\mathcal{H}ol(f)$ is trivial, is not true in general.

\section{Examples of normal holonomy groups and push-out spaces}

\begin{ex}\label{ex.31}
We start with a curve as below
 \begin{center}
  \includegraphics[width=4cm]{1.jpeg}
\end{center}
 suppose this curve is given by $s\mapsto
(\xi(s),\eta(s))$ where $s\in [0,1]$ and at $(1,0,0)$:$s=0$
,${\frac{\partial\xi} {\partial s}}=1,({\frac{\partial}{\partial
s}})^{r}\eta=0$ for all $r\geq 0$ and at
$(0,1,0)$:$s=1$,${\frac{\partial\eta} {\partial
s}}=1,({\frac{\partial} {\partial s}})^{r}\xi=0$ for all $r\geq 0$.
Now,we take this curve in $\Bbb R^3$ and consider the same curves in
$yz$-plane and $xz$-plane and fit together to make a smooth closed
curve in $\Bbb R^3$

Now by identifying $\Bbb S^1$ with $\frac {\Bbb R} {3\Bbb Z}$, the
curve in $\Bbb R^3$  can be redefined as $f:\Bbb
S^1\longrightarrow\Bbb R^3$ where:

$f(s)= \left\{
  \begin{array}{cc}
     (\xi(s),\eta(s),0) \qquad\quad \quad 0\leq s\leq 1 \\
     (0,\xi(s-1),\eta(s-1))  \quad 1\leq s\leq 2\\
     (\eta(s-2),0,\xi(s-2)) \quad 2\leq s\leq3\\
  \end{array}
\right.$

\begin{center}
  \includegraphics[width=5cm]{5}
\end{center}
\end{ex}
To find the normal holonomy group of the above curve, we will
consider normal vector to the curve under parallel transport. As
each part of the curve lies in a 2-plane, the normal plane at a
point of the curve is spanned by the perpendicular direction to the
2-planes.\\
{\bf Step1.}Start with the normal vector at (1,0,0), in the diagram,
it stays in the xy-plane under parallel transport.
\begin{center}

 \includegraphics[width=6cm]{6}
\end{center}
The normal vector (0,1,0) at (1,0,0) goes to normal vector (-1,0,0)at (0,1,0).\\
 {\bf Step2.}At (0,1,0)the the normal vector (-1,0,0) is
 perpendicular to the yz-plane, it stays perpendicular to the
 yz-plane under parallel transport form (0,1,0) to (0,0,1).
 \begin{center}

 \includegraphics[width=6cm]{7}
\end{center}
The normal vector (-1,0,0)at(0,1,0) goes to normal vector(-1,0,0) at (0,0,1).\\
{\bf Step3.}The the normal vector (-1,0,0) is in the xz-plane at
(0,0,1) and stays in the xz-plane from (0,0,1) to (1,0,0).
\begin{center}

\includegraphics[width=6cm]{8}
\end{center}
The normal vector (-1,0,0,) at (1,0,0)  by going once around the
curve the normal vector will turn about $\frac {\pi}{2}$.
\begin{center}

 \includegraphics[width=6cm]{9}
\end{center}
Going around of curve again, the normal vector moves through another
$\frac {\pi}{2}$ and after four times around the curve back to its
original position. this shows that $\mathcal{H}ol(f)$ is generated
by a rotation through $\frac {\pi}{2}$.
\begin{center}

\includegraphics[width=6cm]{10}
\end{center}
Now we find the push-out space of $f$. Except at end-points of three
areas, locally the curve lies in a 2-plane so the focal points with
base $s$,$f(s)+\Sigma_s(f)$, consists of a straight line through the
center of curvature, $c(s)$, of the curve at $s$, perpendicular to
the line joining $c(s)$ and $f(s)$.
\begin{center}

\includegraphics[width=6cm]{12}
\end{center}
At end-points of three areas, and possibly some other points, the
focal set is empty as the center of curvature "at infinity".\\
 so
$\Sigma_s(f)$ is a line in $N_s(f)$. The image of $\Sigma_s(f)$
under normal holonomy group is obtained by rotating it through
$\frac {\pi}{2}$ until it returns to the original position.
\begin{center}

 \includegraphics[width=6cm]{13}
\end{center}
Now, fix the normal plane $N_{s_0}(f)$ at $f(S_0)=(1,0,0)$ where
$s_0=0$ and use parallel transport to identify all the normal planes
with the normal plane  $N_{s_0}(f)$. The push-out space is
complement of all the $\Sigma_s(f)$ and their images under normal
holonomy group.
\begin{center}

 \includegraphics[width=6cm]{14}
\end{center}
Therefore the push-out space of $f:\Bbb S^1\longrightarrow\Bbb R^3$
is an open square $Q$ with sides of length $2\rho$ where $\rho$ is
the minimum absolute value of the radius of curvature of the
original curve in the xy-plane. (i.e. $\Omega(f)$)is the interior of
the smallest square on $N_{s_0}(f)$.)

\begin{ex}\label{ex:3.2}
We consider the immersion $\bar{f}$ as in Example \ref{ex.31} except that
the xz-plane is tilted through an angle $\alpha$.

In other words, $\bar{f}=Lof$ where f is the immersion in example \ref{ex.31} and $L$ is the linear
transformation given by\\
$L=\left(
    \begin{array}{ccc}
    1 & 0 & 0 \\
    0 & 1 & \tan\alpha \\
    0 & 0 & 1 \\
    \end{array}
  \right)$ where
$0<\alpha<\frac {\pi} {2}$.
So, we have
 $$\bar{f}(s)=\begin{cases}
     (\xi(s),\eta(s),0)& 0\leq s\leq 1 \\
     (0,\xi(s-1)+\eta(s-1)\tan{\alpha},\eta(s-1))  &1\leq s\leq 2\\
     (\eta(s-2),\xi(s-2)\tan{\alpha},\xi(s-2)) & 2\leq s\leq3
\end{cases}$$

 The end-points of three areas of this immersion are
(1,0,0),(0,1,0) and (0,$\tan\alpha$,1). Note that at these points
the tangent to the curve is the radial line form (0,0,0) so unit
tangent at (1,0,0) is (1,0,0)and unit tangent at (0,$\tan\alpha$,1)
is
$\frac{(0,\tan\alpha,1)}{\sqrt{1+\tan^{2}\alpha}}$ etc.\\
As in Example \ref{ex.31}, under parallel transport, the normal vector
(0,1,0) at (1,0,0) goes to the normal vector (-1,0,0) at (0,1,0),
which goes to the normal vector (-1,0,0) at (0,$\tan\alpha$,1),
which goes to the normal
vector$\frac{(0,\tan\alpha,1)}{\sqrt{1+\tan^{2}\alpha}}$ at (1,0,0).
So going once around the curve the normal vector has moved through
$\frac {\pi}{2}-\alpha$.
\begin{center}

  \includegraphics[width=6cm]{11}
\end{center}

 This shows that  $\mathcal{H}ol(\bar{f})$
is generated by a rotation through $\frac
  {\pi}{2}-\alpha$
As in Example 3.1, the image of $\Sigma_s(\bar {f})$ under normal
holonomy group is obtained by rotating the line $\Sigma_s(\bar {f})$
through $\frac {\pi}{2}-\alpha$. It depends on $\alpha$ and is
obtained as:
\begin{center}

\includegraphics[width=6cm]{20}
\end{center}
If $R$ is the rotating through an angle $\frac {\pi}{2}-\alpha$ then
 $\Omega({\bar{f}})=\bigcap\{(R)^nQ: n\in\Bbb Z\}$ where Q is a
square as in Example \ref{ex.31}, thus if $\frac {\alpha}{\pi}$ is rational
then $(R)^n(\Sigma_s(\bar {f}))=\Sigma_s(\bar {f})$ for some
$n\in\Bbb Z $ and so,$\Omega({\bar{f}})$  is the interior of the
smallest polygon. If $\frac {\alpha}{\pi}$ is irrational then
$(R)^n(\Sigma_s(\bar {f}))\neq \Sigma_s(\bar {f})$ for any $n\in\Bbb
Z $ and so,$\Omega({\bar{f}})$  is an open disk of radius $\rho$
together with a dense set of points on the boundary circle where
$\rho$ is minimum absolute value of the radius of curvature of the
immersed curve by $\bar{f}$.
\end{ex}

\begin{ex}\label{ex:3.3} In Example \ref{ex:3.2}, we replace the immersion $\bar{f}$ with the
immersion $\bar{f}oh$ where $h:\Bbb R\longrightarrow\Bbb
S^1\equiv\frac {\Bbb R} {3\Bbb Z}$ is covering projection. Since
${\Bbb R}$ is simply connected, for any arbitrary point $s\in {\Bbb
R}$, any closed path at $s$ is nulhomotopic with constant path at
$s$, hence definition \ref{de2} shows that, the normal holonomy
group of $\bar{f}Oh$ is trivial
(i.e.$\mathcal{H}ol(\bar{f}oh)=\mathcal{H}ol_0(\bar{f}oh)$ ). To
calculate $\Omega(\bar{f}oh)$, we prove the theorem \ref{th:3.4}, in
general. It will show that $\Omega(\bar{f}oh)=\Omega(\bar{f})$.
\end{ex}

\begin{thm}\label{th:3.4}
Let $f:M^m\longrightarrow \Bbb R^{m+k}$ be an immersion and
$\hat{M}$ be any covering space with covering projection
$h:\hat{M}\longrightarrow M$. If $\hat{f}=foh$, then
$\Omega(\hat{f})=\Omega({f})$.
\end{thm}
\begin{proof}
 Let $x\in\Omega({f}) $ and fix $p_0\in M$. Then
definition \ref{de3} implies that, $\forall p\in M,\forall \gamma$ s.t.
$\gamma(0)=p_0, \gamma(1)=p$; $\varphi_{p,\gamma}x\neq \Sigma_p(f).$\\
we define the total space of the normal bundle of $\hat{f}$ by
$$N(\hat{f})=\{(\hat{p},x)\in \hat{M}\times\Bbb R^{m+k}:<x,v>=0\quad  \forall v\in
 \hat{ f}_{*}T_{\hat{p}}(\hat{M})\}$$
  Also, for any $\hat{p}\in h^{-1}(p)$
  we have
  \begin{align*}
\hat{f}_{*}T_{\hat{p}}(\hat{M})&=(foh)_{*}T_{\hat{p}}(\hat{M})\\
&=(f_{*}oh_{*})T_{\hat{p}}(\hat{M})\\
&=f_{*}T_p(M)
  \end{align*}
  this shows that, for any $\hat{p}\in h^{-1}(p)$ we have $N_{\hat{p}}(\hat{f})=N_p(f)$
  and so$\Sigma_{\hat{p}}(\hat{f})=\Sigma_p(f)$. Further, we fix$\hat{p_0}\in
  h^{-1}(p_0)$ then
  $\hat{\varphi}_{\hat{p},\hat{\gamma}}=\varphi_{p,\gamma}$ where
  $\hat{\gamma}:[0,1]\longrightarrow\hat{M}$
   s.t. $\hat{\gamma(0)}=\hat{p_0}, \hat\gamma(1)=\hat{p}.$
   Therefore,$\forall \hat{p}\in \hat {M}, \forall \hat{\gamma}$;
   $\hat{\varphi}_{\hat{p},\hat{\gamma}}x \neq \Sigma_{\hat{p}}(\hat{f}).$
Now using definition \ref{de3} again, follows that, $x\in
\Omega(\hat{f})$. By the same way proves that
$\Omega(\hat{f})\subseteq \Omega(f)$.
\end{proof}

\begin{ex}\label{ex:3.5}
If $\frac {\pi}{2}-\alpha=\frac {2\pi}{n}$, then Example \ref{ex:3.3} can be
modified by replacing h by the n-fold covering $\bar{h}:\Bbb
S^1\longrightarrow \Bbb S^1$.

Going once around the first $\Bbb S^1$ in $\bar{h}:\Bbb
S^1\longrightarrow \Bbb S^1$ corresponds to moving n times around
the second $\Bbb S^1$ so parallely transporting a normal n times
around the second $\Bbb S^1$ which gives a rotation of
\begin{align*}
n(\frac {\pi}{2}- \alpha)&=n(\frac{2\pi}{n})\\
&=2\pi
  \end{align*}
  i.e. the identity, so
$\mathcal{H}ol(\bar{f}o\bar{h})=\mathcal{H}ol_0(\bar{f}o\bar{h})$.
Since, the immersed curve by $\bar{f}o\bar{h}$ and the immersed
curve by $\bar{f}$ have same figure in $\Bbb R^3$ and
$\mathcal{H}ol(\bar{f}o\bar{h})$ is trivial so the singular sets of
them also the same (i.e. $\Sigma(\bar{f}o\bar{h})=\Sigma(\bar{f})$).
This implies that $\Omega(\bar{f}o\bar{h})=\Omega(\bar{f})$.
\end{ex}

\begin{ex}\label{ex:3.6}
We consider a sequence of curves $f_n$ in $\Bbb R^3$ defined as in
Example \ref{ex.31} except that $||f_n(s)||$ and the curvature tends to
infinity with n when $s=\frac {1}{2},\frac {3}{2}$ or $\frac {5}{2}$
but is bounded otherwise.\\
\begin{center}

 \includegraphics[width=6cm]{17}
\end{center}
Now, we define the immersion $f:\Bbb R\longrightarrow\Bbb R^3$ by
$f(s\pm3n)=f_n(s)$. When n tends to infinity, the immersion $f:\Bbb
R\longrightarrow\Bbb R^3$ has a sequence of points where the
curvature tends to infinity and the radius of curvature at these
points can be arbitrary small; in other words, $\exists s$ where
$\Sigma_s(f)$ is arbitrary close to "O" in $N_s(f)$. So $\{O\}$ is
the only point not in the image of $\Sigma_s(f)$ under normal
holonomy group for all $s\in \Bbb R.$ Then $\Omega(f)=\{O\}$. In
this case because $\Bbb R$ is simply connected then
$\mathcal{H}ol(f)$ is trivial.
\end{ex}

The following results  have been proved in \cite{2}, when
$\mathcal{H}ol(f)$ is trivial.
\begin{thm}\label{th:3.7}
Let
M be a compact manifold, then each path-connected component of
$\Omega(f)$ is open in $\Bbb R^k$.
\end{thm}

\begin{thm}\label{th:3.8}
Each path-connected component of $\Omega(f)$ is convex.
\end{thm}

\begin{rem}
In Example \ref{ex:3.2}, if $\frac{\alpha}{\pi}$ is irrational then
$\Omega(\bar{f})$ is not open in $\Bbb R^2$ but $M=\Bbb S^1$ is
compact. Also, in Example \ref{ex:3.5} $\Omega({f})=\{O\}$ so
$\Omega({f})$ is closed in $\Bbb R^2$ but $\mathcal{H}ol(f)$ is
trivial. This shows that Theorem \ref{th:3.7} is false when M is not
compact or $\mathcal{H}ol(f)$ is non-trivial.
\end{rem}

\begin{rem}
In Example \ref{ex:3.6} one of the path-connected components of
$\Omega(f)$, which is the complement space of cone and two other
components in $\Bbb R^3$, is not convex. This shows that Theorem
\ref{th:3.8} is false when $\mathcal{H}ol(f)$ is non-trivial.
\end{rem}
Thus the properties of push-out space that are proved in the case
$\mathcal{H}ol(f)$ is trivial, is not true in general.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

On September 11$^{th}$, 2001, World Trade Centers (WTC) 1 and 2 were each hit by a
commercial airplane.
Approximately 56 minutes after being hit WTC 2 collapsed to the ground in what appears
to be a single-front avalanche, which started in the upper half of the building and reached
the ground floor in $_2T_c \simeq 11$ s.
WTC 1 followed the similar destiny some 50 minutes later, where its collapse lasted
$_1T_c \simeq 13$ s.
In 2002, Federal Emergency Management Agency (FEMA) published a WTC performance
study~\cite{fema:2002} that concluded that the buildings were built properly.
National Institute for Standard and Technology (NIST) followed with
its own report on collapse in 2005~\cite{nist:2005}.
The NIST report agreed with the FEMA findings that the buildings did not have design flaws.
It dealt mainly with the conditions throughout the impacted floors and how the
damage from the impact and subsequent fires could have compromised the load carrying
capabilities of the central core of the building in the impact zone.

In this report we use mathematical models of progressive collapse
with the structural parameters of the buildings to examine the
consistency between static and dynamics features of collapse.
Here, under the static features we assume the conditions in the building
on the verge of collapse, while under the dynamic we assume the
observed duration of collapse as calculated using the mathematical model of
progressive collapse.

The report is organized as follows.
In Sec. II we present the mathematical models
of progressive collapse, where we introduce two classes (discrete and continuous)
and two types (elastic and plastic) of models.
We briefly discuss the underlying physics and how the structural
parameters enter each class of the models.
In Sec. III we discuss the structural parameters of WTC 1 and 2 and introduce
their two main load-bearing structures, the weaker perimeter columns (PCs)
and the stronger core columns (CCs).
In Sec. IV we discuss the NIST scenario of collapse, and show how it is inconsistent
with the static features of collapse (collapse initiation lines) unless
one neglects the contribution of CCs.
In Sec. V we present our conclusions, while in Sec. VI we discuss some less obvious
assumptions.
One group assumptions is hidden behind the one-dimensional formulation of progressive
collapse, while the other behind, so-called, rigidity assumption, by which the range
of elastic/plastic deformation of a column does not spread beyond the stretch of a
single storey.
In Appendix A we present mathematical model of collapse of WTC 7 and analyze it
in the same way as it was done for WTC 1 and 2.

\section{Mathematical Models of Progressive Collapse}

The World Trade Centers 1 and 2 were complex three-dimensional structures, yet it appears
that the progressive collapse in each of them was one-dimensional:
an avalanche formed at the top of the impact zone and propagated to
the ground.
Based on that fact, we describe two classes of the mathematical models
that were proposed for description of the progressive collapse of the buildings:
continuous and discrete.
An assumption, both models critically depend on, is that the avalanche front -
a fictional plane that moves with the velocity of the avalanche and through
which the building underneath joins the avalanche - maintained its
shape (flatness, being parallel to the ground) for the whole duration of collapse.

Less important assumptions are those regarding the loss of mass and momentum
and the rigidity.
As for the losses, in the case of a massive building, such as the WTCs, one could expect
lateral (material leaking on the sides) and longitudinal
(significant variation in resistive force along the avalanche front prevents formation
of a sufficiently compacted layer that would be able to crush strong vertical columns
in its path).
The rigidity assumption, on the other hand, requires that the interaction between
the avalanche and the building is localized only to the point of contact, and not
any further into the building beneath.
In our exposition of the mathematical models we
thus assume no losses of any kind, and a reasonable rigidity in the building,
all of which contributed to the apparent one-dimensionality of the avalanche propagation.

\subsection{Continuous Models}

For the purpose of our models we assume that the avalanche propagates
along the $z$-axis, where the positive direction is along the direction
of Earth's gravity.
The building in which the avalanche propagates thus initially stretches
from 0 (top) to H (bottom), and its longitudinal mass density
is described by $\mu=\mu(Z)$, where $M=\int_0^H\,dZ'\ \mu(Z')$ is the
total mass of the building.

In the continuous class of mathematical models of progressive collapse, the avalanche
at position $Z$ consists of two objects moving together: the top section,
and the avalanche front.
The top section is the part of the building initially stretching from 0 to $Z_0$,
where $Z_0$ is the point at which the avalanche started.
The avalanche front is a zero-size object that contains all the mass the avalanche
has collected during its fall.
The total mass of the avalanche is thus the mass of the top section
and the mass of the avalanche front, and is given by $m(Z) = \int_0^Z\,dZ'\ \mu(Z')$.
The destruction of the building occurs at the interface between the
avalanche front and the stationary part of the building, where all the
transfer of mass and momentum between the two occurs.

We consider the perfectly inelastic (plastic) collisions between the avalanche
front and the building first.
Here, the rate of change of momentum of the avalanche follows from elementary
Newtonian dynamics, and is given by
\begin{equation}
  \label{eq:nc}
  \frac{d}{dt} \left( m(Z) \, \dot Z \right) = m(Z) \, g + {R^*}(Z),
\end{equation}
where $\dot Z = dZ/dt$ is the velocity of the avalanche.
In \eqn{eq:nc} the term ${R^*} = {R^*}(Z)$ comes from the average resistive force
produced by the building.
It measures how does the building oppose its own destruction by the avalanche,
and is intimately related to, so called, the yield strength of the
structural members that provide the vertical support to the building.
The term $m(Z) \, g$ is the usual gravitational force, where $g$ is Earth's gravity.
In this formulation the resulting equation resembles that of a falling
water drop~\cite{ews2001,krane1981} or a falling chain.~\cite{keiffer2001,thornton:2004}

The collisions between the avalanche front and the building can also be
considered as perfectly elastic.
This assumption leads to the Lagrangian formalism, which requires
that all relevant energies be expressible in terms of a small set of
configuration coordinates and their time derivatives (generalized velocities).
For this particular problem position of the avalanche, $Z$, and
its velocity, $\dot Z$, are the most obvious choice.
The relevant energies necessary for Lagrangian are as follows:
$K(Z,\dot Z) = \frac 1 2 \, m(Z) \, \dot Z^2$ is the kinetic energy, while
$U(Z) = - m(Z) \, g \,Z - \int_Z^H dX \, \mu(X) \, g \, X$ is its
potential energy.
Additionally, we need energy associated with the structure of the building,
$L = - \int^Z\,dX \ R^*(X)$, as a result of which the resistive
force appears.
In Lagrangian formulation the equation of motion follows from
$\frac {d}{dt} \, \partial{\cal L}/\partial \dot Z = \partial{\cal L}/\partial Z$,
where ${\cal L} = K - U - L$ is the Lagrangian:
\begin{equation}
  \label{eq:c}
  \frac{d}{dt} \left( m(Z) \, \dot Z \right) = m(Z) \, g
  + \frac 1 2 \, \mu(Z) \, \dot Z^2 + R^*(Z).
\end{equation}
A practical physical realization of such system is a standard
textbook problem~\cite{wy2005,thornton:2004} of
a falling chain.
Please note, \eqn{eq:c} can also be derived from \eqn{eq:nc} by considering
not only a mass transfer but also a transfer of momentum.~\cite{thorpe1962,tiersten1969}

\eqns{eq:nc} and \noeqn{eq:c} can be jointly written as
\begin{equation}
  \label{eq:cnc}
  \frac{d}{dt} \left( m(Z) \, \dot Z \right) =
  m(Z) \, g  + \frac 1 2 \, \epsilon \, \mu(Z) \, \dot Z^2 + R^*(Z),
\end{equation}
where we introduce an additional parameter, elasticity, $\epsilon\in[0,1]$,
which describes the type of collisions between the building and the avalanche front,
of which plastic and elastic collisions are just the extremae,
$\epsilon=0$ and $\epsilon=1$, respectively.
One can envision obtaining the value of $\epsilon$ by fitting the
actual progression of the collapse in time to that predicted by the
model~\eqn{eq:cnc}.

To simplify the analysis of progressive collapse further we assume that
the building is homogeneous height-wise, i.e., $\mu = M / H$.
Finally, we convert \eqn{eq:cnc} to dimensionless form where we use
as the unit of length the height of the building $H$, and as the
unit of time the free fall time from height $H$, $T=\sqrt{2\, H / g}$.
The position of the avalanche $z=Z/H$ in time $\tau=t/T$ is described as a solution
of the following ordinary differential equation:
\begin{equation}
  \label{eq:cnc:ode}
  z''(\tau) = 2 - \left(1 - \frac{\epsilon}{2}\right) \,
  \frac{z'^2(\tau)}{z(\tau)} + \frac{2}{z(\tau)}\,\frac{R^*(z)}{M\,g},
\end{equation}
where prime denotes the differentiation with respect to dimensionless time $\tau$.
We elaborate on the initial conditions $z(0)$ and $z'(0)$ later, while we
discuss the resistive force $R(z)$ next in the context of the discrete models.

\subsection{Discrete Models}

In the discrete models the building is represented by a periodic structure
of stacked stories,~\cite{bazant2006}
where each storey consists of a floor and a number of vertical
columns which keep the stories at their positions
$Z_F$, where $F=1,..F_T$ is the floor number, while $F_T$ is the total number
of floors in the building.
The height of each storey is $\Delta H = H/F_T$.

Just like in one-dimensional continuous model \eqn{eq:cnc}, the avalanche
in discrete model propagates only in vertical direction.
Let us assume that the consuming of a storey by the avalanche comprises of
three phases, {\em resistive}, {\em free fall} and
{\em adsorption},
where each phase takes a definite amount
of fractional floor height, $\lambda_1$ for resistive,
$\lambda_2$ for free fall, and $\lambda_{\infty}$ for adsorption.
Obviously, $\lambda_1 + \lambda_2 + \lambda_\infty = 1$.
The fractional height $\lambda_\infty$ is sometimes called the compaction
ratio.~\cite{bazant}
During first two phases the avalanche moves without acquiring the mass,
where the resistance of the vertical columns is $R^1$ and $R^2\equiv0$, respectively.
After dropping by the fractional height $(\lambda_1+\lambda_2)$ the avalanche
makes a discontinuous jump by $\lambda_\infty$ and adsorbs
the compacted storey in the process.
The progressive collapse of the $j$-th storey of
mass $\Delta m = M/F_T$ is described by the following
set of dimensionless equations:
\begin{subequations}
  \label{eq:d}
  \begin{equation}
    \label{eq:d:1}
    \left.
      \begin{array}{lcl}
        a_k &=& 2 + \frac{2}{j-1} \, \frac{R^k_{j-1}}{M\,g}\\
        u_k^2 &=& u_{k-1}^2 + 2 \, \frac{\lambda_k}{F_T} \, a_k\\
        \Delta \tau^k   &=& \frac{u_k - u_{k-1}}{a_k}
      \end{array}
    \right\} \mbox{ for } k=1,2,
  \end{equation}
  \begin{equation}
    \label{eq:d:2}
    v_j = \frac {j-1 + \frac{\epsilon}{2}}{j} \, u_2,
  \end{equation}
\end{subequations}
where $j = F_T - F$, with $F$ being the floor number of the current avalanche position.
Here, $u_0 = v_{j-1}$ is the velocity of the avalanche at the beginning of the
storey, while $\tau_j -\tau_{j-1} = \Delta \tau^1  + \Delta \tau^2$ is how long it took
the avalanche to transverse the storey.
In \eqn{eq:d:2} we introduce a parameter $\epsilon$ in the same
way it appears in the continuous model \noeqn{eq:cnc},
\begin{equation}
  \label{eq:d:e}
  \epsilon = 2 \, \frac{u'}{u_2},
\end{equation}
where $u'$ is the average velocity of the compacted storey at the moment it
is adsorbed by the avalanche, the velocity of which is $u_2$, see \eqn{eq:d:1}.

In \eqns{eq:d:2} and \noeqn{eq:d:e} one can see that the rigidity assumption
put forth by Bazant {\em et al.}\cite{bazant2000} is expected to hold, at best, only for
the height scales greater than the storey height, but not for smaller.
The compaction of the storey by the avalanche begins the moment the avalanche enters
the storey, and by the time the avalanche reaches the position $\lambda_1+\lambda_2$
some of the storey's mass is already moving with the velocity of
the avalanche $u_2$, so in fact $\epsilon>0$.

\subsection{Resistive Force: discrete vs. continuous}

The resistive force $R$ is the force the building opposes its own destruction
by the avalanche.
Bazant {\em et al.}\cite{bazant2006} call it the ``crushing force,'' and for the most part of their
calculation they assume that its magnitude is a constant with respect to height.
For our investigation we relax that assumption by allowing
$R$ to vary with height in a simple linear fashion,
\begin{equation}
  \label{eq:R}
  f(z) = -\frac{R(z)}{M\,g} = r + s \cdot z,
\end{equation}
where $M$ is the total mass of the building, and $z=Z/H$ is the scaled height.
As we will see later, when we discuss the load-bearing capacity of the structure
of WTC 1 and 2, this will turn out to be more accurate description.

Let us now discuss the resistive force within the discrete (floor) model.
We posit that the building is rigid, i.e., that the floor that avalanche
is currently penetrating sits on an infinitely strong structure.
With that assumption the interaction between the building and the avalanche
is localized to the floor of contact between the two and not beyond.
Now, as the avalanche enters the storey it first encounters the resistive force from
the vertical columns, call it $R^1$.
Over the fractional length $\lambda_1$ the vertical columns
maintain ultimate yield force under compression.
In other words, for the first phase of floor destruction
the fractional distance $\lambda_1$ corresponds to the yield strain $\lambda_Y$,
while $R^1$ corresponds to the ultimate yield force.

Assuming that avalanche has passed the fractional distance $\lambda_1$,
the columns fail.
The failure mode of the column is debatable:
if the columns  are compressed by the avalanche than they fail by buckling. However, if
the avalanche front consists of crushed material, so that the ends of the vertical
columns facing the avalanche are free, then the failure mode may be bending, as well.
Either way, we assume that during this phase the column does not offer any resistance.
In other words, for the fractional distance of $\lambda_2$ the avalanche falls freely,
$R^2 \equiv 0$.

The last phase of propagation occurs when the avalanche reaches the fractional distance
of $\lambda_1+\lambda_2 = 1 - \lambda_\infty$, that is, it compacts the storey.
We assume that then the position of the avalanche makes a discontinuous jump
from $1 - \lambda_\infty$ to 1 (bottom of the current storey) and adsorbs the mass
of the storey as discussed earlier, see \eqn{eq:d:2}.

The resistive force $R^*=R(r^*,s^*)$ that enters the continuous models
is related to $R=R(r,s)$ in the discrete model via the energy argument.
By that argument the energy $\Delta L$ over the fractional length
$\lambda_1+\lambda_2$ is the same in both models, yielding
\begin{subequations}
  \label{eq:rc:all}
  \begin{equation}
    \label{eq:rc:1}
    r^* = \frac{\lambda_1}{1-\lambda_\infty} \, r,
  \end{equation}
  \begin{equation}
    \label{eq:rc:2}
    s^* \simeq \frac{\lambda_1}{1-\lambda_\infty} \, s.
  \end{equation}
\end{subequations}
As we have two descriptions of the resistive force, we choose the discrete
as the primary.
Thus, whenever we talk about the parameters of the building we always assume
those entering the discrete model, thus $(r,s)$.
If we use continuous models for what ever reason then we imply that
the nominal, or discrete, values $(r,s)$ were converted to $(r^*,s^*)$
via~\eqn{eq:rc:all} before commencing the computations.

A question may arise if the two models for the same nominal values of $r$ and $s$
give the same or, at least, similar results.
The question is not so much relevant from the computational point of view:
recursive relations~\noeqn{eq:d} comprising the discrete model are actually simpler
to numerically implement than the ODE comprising the continuous model~\noeqn{eq:cnc:ode}.
On the other hand, the continuous model allows us to use a number of theoretical
tools (Lagrangian formalism, energy analysis, \& c.).
The main difference between the models is in the distance the avalanche
propagates.
In the continuous model the avalanche covers the full height of each storey, $\Delta H$.
In the discrete model, on the other hand, in each storey a section of
fractional length $\lambda_\infty$ is excluded from propagation,
so the effective height per story is $(1-\lambda_\infty) \, \Delta H$.
For the quantity we are interested in, which is the duration of collapse
or {\em collapse time}, we expect both models to be similar for the following
reason.
The collapse time $t$ depends on the velocity $u$ and the distance $d$ as
$t \sim d /u$.
When we change from discrete to continuous model the path increases but so does the velocity
as there is more energy available.
In the lowest order of approximation the two effects are thus expected to
counteract each other.
In fact, we use the spread between the two models as a measure of uncertainty
when determining $r$ and $s$ that yield desired (observed) collapse time.

On the other hand, by comparing plastic ($\epsilon\equiv0$) to
elastic case ($\epsilon\equiv1$) the following is obvious: given $(r,s)$
the collapse time is shorter in elastic than in plastic case.
This difference, as is common knowledge, comes from the energy dissipation inherent
in plastic model.
In spirit of ``no-energy-and-momentum-left-behind,'' by which we try to neglect as many
losses as possible to obtain the lower estimates on collapse times,
we assume $\epsilon\equiv1$ for most of the computations that follow.

\subsection{Bazant's model}

Bazant {\em et al.}\cite{bazant2006} proposed the following mathematical
model to describe the progressive collapse in World Trade Center 2,
\begin{equation}
  \label{eq:bazant:1}
  \left(1 - \lambda_\infty \right) \,
  \frac{d}{dt} \left( m(Z) \, \dot Z \right) = m(Z) \, g + {R}(Z),
\end{equation}
where, as before, $Z$ is the position of the avalanche front and $R(Z)$ is the
average floor ``crushing'' force.
Additionally, the resistive force that enters the computations is
to be halved along the whole length of the building due to the ``heat'' effects
(see p.15 on-line edition, top paragraph).
Finally, for the most part of their presentation the authors assume that the
computational value for the parameters $(r,s)$ of the resistive force
are $s=0$ and $r=0.15$
(for crushing energy of $W_f = 2.4$ GNm, $\Delta H=3.7$ m, where the total mass
of the building is $M=450\cdot10^6$ kg).
Based on their modeling the authors concluded ({\em sic!}) that there was nothing
unusual in the way the building collapsed.

The rest of this paper is to show that the alternative interpretation
is possible.
As already discussed, the scaling of kinetic term on the right-hand-side
in \eqn{eq:bazant:1} is unnecessary to equate the discrete
to the continuous model.
In fact, this procedure is equivalent to boosting the force that accelerate
the avalanche by a factor of $(1 - \lambda_\infty)^{-1}$.
Similarly, the model \eqn{eq:bazant:1} is non-conservative to embody rigidity
assumption against which, the reader might recall, we argued.
Finally, no attempts were made on authors' behalf to correlate the information about
the status of the building to the values of the parameters that enter the model.
Heat argument that the authors used to halve the magnitude of the resistive force in
the whole building is in direct contradiction with the NIST report, which explicitely
states that no temperatures higher than 200$^o$C were measured anywhere in the building
below the impact zone.

Overall, it appears that the authors tried to compensate for slowness of their
non-conservative model by boosting the energy $U+L$ and arbitrary halving
the strength of the building.
For these reasons, we exclude Bazant's model from the analyses that take the
rest of the paper.

\section{Structural parameters of the World Trade Centers}

World Trade Centers 1 and 2 were designed as the external tubular frame around
the strong central core.~\cite{web:wtc}
Tubular frame consisted of 236 perimeter columns (PC), while
the central core consisted of 51 core columns (CC), all from structural
steel of nominal strength 36-100 KSI.
These provided load-bearing support to the whole structure, while the
combination of steel trusses and concrete panels stretching between the
two groups of columns provided actual floors for occupants to dwell on.
We recall that the destruction is determined by the ultimate, rather than the
nominal strength of the vertical columns.

PCs, the external dimensions of which were 14''-by-14'', were made of structural
steel the yield strength of which varied from nominal 36 KSI
(ultimate 58 KSI) and thickness $1\over4$'' at
the top of the building to nominal 100 KSI (ultimate 110 KSI)
at the bottom.\cite{nist:2005}
The thickness of the steel plates at the bottom, to the best of author's
knowledge, is not yet publicly available.
CCs were made of structural steel which varied from nominal 36 KSI (ultimate 58 KSI)
at the top to nominal 42 KSI (ultimate 60 KSI)  at the bottom.
Neither the dimensions or the thickness of the steel plates used for CCs,
to the best of the author's knowledge, are yet publicly available.
The NIST report~\cite{nist:2005} claims that there were two types of
CCs: ``standard'' and ``massive,'' where the four ``massive'' columns,
one at each core's corner, ``together provided 20\% of load-bearing
capacity of the core.''
This was recently proven to be misleading:\cite{911:core}
information released to the public shows that of 51 CCs at least 16 (along two longer
sides) were of external dimensions 22''-by-55'' and of thicknesses
up to 5''.
To obtain the values for our models we proceed as follows.
For each structural element, the PCs and the CCs, we find the scaled ultimate
yield force $f=f(z)$, at the top ($z=0$) and at the bottom ($z=1$).
The parameters $r$ and $s$ are then obtained using
\begin{subequations}
  \label{eq:rs:def}
  \begin{equation}
    \label{eq:rs:def:r}
    r = f(0),
  \end{equation}
  \begin{equation}
    \label{eq:rs:def:s}
    s = f(1)-f(0).
  \end{equation}
\end{subequations}

For the PCs we assume that the thickness of the plates at the bottom was $1\over2$'',
yielding
\begin{equation}
  \label{eq:pc:yield}
  f_{PC}(z) = \frac{R_{PC}}{M\,g} \simeq 0.2 + 1.2 \cdot z,
\end{equation}
where $M = 4.5\cdot10^8$ kg is the estimated mass of the building, while
$g$ is gravity.

As for the CCs, we do base our expectations
on the following: all 51 column were of dimensions 22''-by-55'', while
their thickness and strength varied from 1$1\over4$'' and nominal 36 KSI
(ultimate 58 KSI) at the top, to 5'' and nominal 42 KSI (ultimate 60 KSI)
at the bottom.
This yields,
\begin{equation}
  \label{eq:cc:yield}
  f_{CC}(z) = \frac{R_{CC}}{M\,g} \simeq 0.6 + 1.5 \cdot z.
\end{equation}
Please note, the assumption of all 51 columns being the same appears to
be corroborated by the floor plans:~\cite{nist:2005} all CCs appear to
have the same footprint.
Also, while this estimate is purely hypothetical, it turns out that,
it is also irrelevant for the analysis of collapse that follows.

We observe that from \eqns{eq:pc:yield} and \noeqn{eq:cc:yield}
the ultimate safety factor of the WTCs is $f_{CC}(1)+f_{PC}(1) \sim 3.5$.
This is a reasonable estimate considering the dimensions of the buildings
and the other safety requirements that entered the structural calculation
(ability to withstand hurricane winds and an airplane impact).

From the properties of structural steel\cite{materials} it is known that the yield
strain under tension and compression are fairly similar, and is $\sim 21-25\%$.
In our model this is represented by $\lambda_1$, which we take to be $\lambda_1 = 0.2$.
The value of compaction limit we take from Bazant,\cite{bazant2006}
$\lambda_\infty=0.2$, which leaves $\lambda_2 = 1-\lambda_1-\lambda_\infty = 0.6$.
From there, $(r^*,s^*)$ in the continuous model are related to
 $(r,s)$ in the discrete values as,
\begin{equation}
  \label{eq:rs:wtc}
  \begin{array}{lcl}
    r^* &=& 0.25 \cdot r,\\
    s^* &=& 0.25 \cdot s.
  \end{array}
\end{equation}

\section{Failure of the NIST scenario}

As is known, WTCs 1 and 2 were subjected to two damaging events, as a result of
which the buildings eventually collapsed in ${_1T}_c=13$ s (WTC 1),
and ${_2T}_c=11$ s (WTC 2).
The first event was an airplane collision in which some of the vertical columns
were cut or otherwise compromised.
We label that event with a letter $M$, for {\em Mechanical}.
The second event was a gradually developing heat damage from ambient and
jet-fuel fires.
We label that event with a letter $H$, for {\em Heat}.
Both types of damage were delivered to each building only in the
primary or impact zone, while the secondary zone (below the primary)
was left intact.
We note that the primary zone in WTC 1 spread over floors $F_0=99$ to $F_1=93$,
while in WTC 2 it spread between the floors $F_0=85$ and $F_1=77$.
Both buildings had the
same number of floors $F_T=110$.
\begin{figure}[h]
  \singlespacing
  \centering
  \includegraphics[scale=0.4,clip]{fig1.eps}
  \caption{%
    \label{fig:wtc}
    Schematic of airplane collisions with WTC 1 and 2 on September 11,
    2001~\cite{web:wtc:scheme}.
  }
\end{figure}

Let us now formulate, what we call, the NIST scenario for the collapse of
WTC 1 and 2.
For that we need to define first {\em the primary zone damaging event} or NIST($\mu$):
\begin{quote}
  {\bf Definition:} \em A primary zone damaging event, NIST($\mu$), is an incident
  as a result of which the nominal strength of the resistive force in the primary zone
  of the building, $\tilde R$, is reduced by a factor $\mu\ge1$ with respect to
  the value in intact building, $R=R(z)$.
\end{quote}
This is written as,
\begin{equation}
  \label{eq:nist:mu}
  \tilde R(z) = \left\{
  \begin{array}{ll}
    \frac{1}{\mu}\, R(z), & \mbox{ for } z_0 < z < z_1,\\
    R(z), & \mbox{ for } z_1 > z,
  \end{array}
  \right.
\end{equation}
where $z_i = 1 - F_i/F_T$, for $i=0,1$.
The NIST scenario regarding the collapse of WTC 1 and 2 is thus formally expressed
as,
\begin{quote}
  {\bf Definition:} \em
  A NIST scenario of collapse, or NIST($\mu \cdot \nu$), is a sequence
  of two primary zone damaging events to which a building is exposed:
  $M$, with $\mu\ge1$, and $H$, with $\nu\ge1$.
  As a result of the first event, NIST($\mu$), the building does not collapse.
  It collapses only after the second primary zone damaging event,
  NIST($\mu \cdot \nu$), where the duration
  of collapse is known to be $T_c$.
\end{quote}

Let us estimate factors $\mu$ and $\nu$, as they are essential for the NIST scenario.
The NIST report~\cite{nist:2005} documented the damage to the primary zone
fairly well.
WTC 1 suffered frontal hit, and after the aircraft penetrated the building it
exploded inside.
We assume that for WTC 1 we have
\begin{equation}
  \label{eq:mu:wtc1}
  \frac{1}{_1\mu} = 0.5,
\end{equation}
that is, we assume that the aircraft destroyed 50\% of the PCs and of the CCs.
WTC 2, on the other hand, was hit in the corner.
Parts of the airplane flew through the building and most of the fuel
exploded outside the building.
We thus assume that for WTC 2 we have
\begin{equation}
  \label{eq:mu:wtc2}
  \frac{1}{_2\mu} = 0.75,
\end{equation}
that is, we assume that the aircraft destroyed 25\% of the PCs and the CCs.
As for the heat damage, it is reported\cite{nist:2005} that the temperatures of up
to 600$^o$C were measured at some of the locations in the primary zone.
As is known, at that temperature the structural steel loses approximately
one-half of its strength.
We thus have
\begin{equation}
  \label{eq:nu:wtcs}
  \frac{1}{\nu} = 0.5,
\end{equation}
for both buildings.

Now that we have $\mu$ and $\nu$ for each building, we recall that the
collapse initiation occurred at the top of the primary zone,
at position $z_0 = 1 - F_0/F_T$, where for
WTC 1, ${_1z}_0 =0.1$, and for WTC 2, ${_2z}_0 = 0.23$.
This allows us to construct the collapse initiation lines.
The collapse starts at point $z_0$ because the yield force of
the compromised building, $\tilde  R(z_0)$, is not sufficient to resist the
weight of the building above,
\begin{equation}
  \label{eq:damage}
  -\tilde R(z_0) < m(z_0) \cdot g.
\end{equation}
Please observe, the collapse initiation lines are derived from static properties of the
buildings, thus they do not depend on the mathematical model
used to describe the dynamics of collapse.

The collapse initiation lines according to the NIST scenario of collapse
are as follows:
\begin{itemize}
\item
  {\em 0-damage line}:
  \begin{equation}
    \label{eq:0-damage}
    \frac r {{_j z}_0} + s = 1.
  \end{equation}
\item
  {\em M-damage line}:
  \begin{equation}
    \label{eq:M-damage}
    \frac r {_j\mu \cdot {_jz}_0} + \frac{s}{_j\mu} = 1.
  \end{equation}
\item
  {\em M+H-damage line}:
  \begin{equation}
    \label{eq:M+H-damage}
    \frac r {{_j\mu} \cdot \nu \cdot {_jz}_0} + \frac{s}{{_j\mu} \cdot \nu} = 1.
  \end{equation}
\end{itemize}
where $j=1,2$, for WTC 1 and 2, respectively.
Please note, in \eqns{eq:0-damage}-\noeqn{eq:M+H-damage}, $r$ and $s$ represent
the intact building parameters.
Obviously, these lines limit $r$ and $s$ the buildings could have had.
E.g., the buildings do not collapse originally, thus their $r$ and $s$ are
somewhere above the 0-damage line,  \eqn{eq:0-damage}.
Similarly, the building does not collapse after an $M$-event,
thus $r$ and $s$ are somewhere above the $M$-damage line,
\eqn{eq:M-damage}.
On the other hand, the buildings do collapse after a combined $M+H$-event,
thus $r$ and $s$ are somewhere between the $M$-damage line and
the $M+H$ damage line, \eqn{eq:M+H-damage}.

\begin{figure}[htp]
  \singlespacing
  \centering
  \includegraphics[scale=0.7,clip]{fig2a.eps}
  \includegraphics[scale=0.7,clip]{fig2b.eps}
  \includegraphics[scale=0.7,clip]{fig2c.eps}
  \includegraphics[scale=0.7,clip]{fig2d.eps}
  \caption{%
    \label{fig:2}
    Collapse time in the $NIST(\mu\cdot\nu)$ scenario,
    in the discrete model (top row) and in the continuous model (bottom row),
    in WTC 1 (left column) and WTC 2 (right column).
    The base contains the contours $T(r,s) = {_j\tau}_c, 2$ for each building.
  }
\end{figure}
\begin{figure}[htp]
  \singlespacing
  \centering
  \includegraphics[scale=0.5,clip]{fig3a.eps}
  \includegraphics[scale=0.5,clip]{fig3b.eps}
  \caption{%
    \label{fig:3}
    Collapse initiation lines (black from \eqn{eq:0-damage},
    blue from \eqn{eq:M-damage}, purple from \eqn{eq:M+H-damage}),
    with the collapse time $\tau(r,s)=\tau_c$ as calculated
    using the continuous (solid red and orange) and the
    discrete model (dashed red and orange), for $NIST(\mu)$ (red) and
    $NIST(\mu\cdot\nu)$ scenario (orange).
    Both collapse duration lines, $T^{M}=\tau_c$ (red), and
    $T^{M+H}=\tau_c$ (orange), are consistent with the collapse
    initiation lines, $M$-damage line (blue) and $M+H$-damage line (violet),
    respectively.
    The collapse duration lines do not prefer either scenario, and do not
    reveal new information except that the heat damage was below the
    assumed (maximal) heat damage.
    The points representing the structural parameters of the building,
    known ultimate yield force of PCs (green), guessed for CCs (green, with
    question mark) and guessed for their sum (pink, with question mark),
    show that the collapse is consistent with the resistive force coming
    from the PCs only.
  }
\end{figure}
Now we use the mathematical models of progressive collapse and calculate the
collapse duration lines
${_j\tau}^M(r,s) = {_j\tau}_c$ and ${_j\tau}^{M+H}(r,s) = {_j\tau}_c$, for $j=1,2$,
where we assume that the avalanche started at the top of the primary zone, $z(0) = {_jz}_0$,
with zero velocity, $z'(0) = 0$.
The total collapse times are scaled with the free fall time $T = \sqrt{2\,H/g} = 9.22$ s,
from the height of the building, $H=417$ m,
so that ${_2\tau}_c = 1.41$ (${_2T}_c \simeq 13$ s, for WTC 2) and ${_1\tau}_c = 1.19$
(${_1T}_c \simeq 11$ s, for WTC 1).
In our mathematical models we use $\epsilon\equiv1$, as this provides the shorter
collapse times than $\epsilon\equiv0$.
We show in \fig{fig:2} collapse time as a function of building
parameters $r$ and $s$ in each building for the $NIST(\mu\cdot\nu)$ scenario.

In \fig{fig:3} we plot the collapse duration lines (from discrete and continuous model)
for $NIST(\mu)$ and $NIST(\mu\cdot\nu)$ scenarios, the collapse initiation lines
and the points that correspond to the known ultimate
yield strength of the perimeter columns (PCs), the guessed ultimate yield
strength of the core columns (CCs, with question mark), and the
guessed ultimate yield strength of both (CC+PC, with question mark).
We observe that the position of the of the collapse duration lines,
${_j\tau}_c = \tau(r,s)$, varies somewhat with respect to the class of the model
(discrete or continuous), but the difference is not critical.
In fact, we use this spread to estimate the approximate range
of $(r,s)$ that yield the observed collapse time $\tau_c$.

For the $M$-event (mechanical damage from the airplane impact),
we observe that the respective collapse initiation line overlap with the
collapse duration lines.
On its own, the overlap implies that the spontaneous collapse of the building
as a result of mechanical damage only is not impossible.
However, because $r$ and $s$ corresponding to the $M$-event are well below the
estimated strength of the PCs
The PCs in each building's secondary zone appear to have
been non-compromised prior to collapse, so this must be a minimum resistive
force the building could have provided.
Thus, we conclude that the buildings survived the $M$-event, which in fact they did.

For the $M+H$ event (heat damage that followed mechanical damage), in the expected
range of parameters ($r$ and $s$ above that of PCs) the collapse initiation line
does not overlap with the collapse duration lines.
This can be interpreted in two ways:
first one is that the collapse was not spontaneous, while
the second is that the reduction of strength due to heat damage
was smaller than presumed, that is, $\nu < 2$.
We observe that, if the latter were the case, by decreasing $\nu \rightarrow 1$,
the lines for the $M+H$ event (collapse initiation and duration) move toward the same
lines for the $M$ event.
However, the collapse duration lines move ``slower''
than its $M+H$-damage line, so for some $\nu < 2$ the two can overlap.

More importantly, we find that the $r$ and $s$ of the PCs are located
in the region bounded by the collapse duration lines.
This means that in each building the collapse initiation and duration are consistent
with the NIST($\mu\cdot\nu$) scenario being applied to the PCs only,
while the stronger core columns (CCs) are not present at all.
This in turn implies that the NIST scenario is incomplete:
the collapse of the buildings to the ground requires yet another damaging event,
the sole purpose of which is a destruction of the CCs in the secondary zone.
We label this damaging event the ``wave of massive destruction'' (WMD), because of
its catastrophic nature.
Interestingly, the avalanche we have discussed so far can only appear in its wake,
and is thus a result of the WMD rather then the other way around.

\section{Conclusion}

We have determined the static and the dynamic features of a progressive collapse
in the WTCs using the structural properties of the building and the
mathematical models of the avalanche propagation.
We have formally expressed the destruction scenarios proposed by NIST
as a sequence of damaging events in the primary (or impact) zone of each building,
which leave the secondary zone (below) intact.
We have shown that the static and dynamic features of collapse are mutually consistent.
On the other hand, we have demonstrated that the NIST scenarios are inconsistent
with the structural parameters of the building.
More precisely, the features of the avalanche propagation (initiation and duration)
indicate that in their final moments the buildings did not have the core columns (CCs).
We conclude that the buildings did not perish because of combined
mechanical and heat damage to their primary zones,
but because of yet another catastrophic event:
a wave of massive destruction (WMD) that destroyed the CCs, following which
the buildings collapsed to the ground.

\section{Discussion}

A physical situation where a building is being pushed to the ground by the impact of its top
section, where its structural strength comes from a few vertical elements,
is implausible.
Consider the following:
\begin{itemize}
\item[1.] {\bf Destruction of vertical columns:}
Floor plans of the WTCs show floors shaped as flat rectangular doughnuts.
The perimeter columns (PCs) are on the outside, while the core columns (CCs)
are on the inside.
Each floor consists of a concrete surface supported by trusses stretching between
the vertical columns.
If we thus think of the avalanche front as being made of the floor material,
we see that the avalanche front does not stretch horizontally far enough
to reach the PCs and CCs.
In the proposed model of collapse, however, it is implied that $(i)$ the
avalanche front is wide enough to reach both columns, and $(ii)$ provides
sufficient pressure at the edges so that the vertical columns fail.
These assumptions are well hidden in the one-dimensional formulation of our
mathematical models:
all the pressures are integrated over the perpendicular cross section
of the building, and only then the equations of motion are derived.

Consider the average pressure created by the avalanche at the bottom of
the primary zone, $z_1 = 1-F_1/F_T$, which is given by
\begin{equation}
  \label{eq:comment:1}
  p = \frac {M\, g}{a^2} \cdot  z_1.
\end{equation}
Here, $a = 206'=2472''$ is an approximate length of the side of the building.
We get ${_1p} = 0.025$ KSI for WTC 1, and ${_2p} = 0.05$ KSI for WTC 2,
which is three orders of magnitude smaller than nominal 36-100 KSI (ultimate 58-110 KSI)
the vertical columns were able maintain while yielding in plastic deformation.
Bazant {\em et al.}\cite{bazant2006} argued that an avalanche propagating through
the primary zone would get sufficiently compacted so that it
could provide necessary pressure.
We see two insurmountable problems with this suggestion.
First, the avalanche front can only ``grow'' thicker - it cannot
expand laterally in such a fashion that would allow its edges to be strong enough to
crush the vertical columns.
Second, for compaction to happen the floor material has to be compressed between two
solid surfaces, and we see that there are no such surfaces on either end of the avalanche
front.
In fact, the strength of the vertical columns will redirect the avalanche (which now
consists only of destroyed floor material) to the region in-between the columns.
The formation of such avalanche is promoted by the relative
weakness of the floors, the resistive force $f$ of which is
$f \sim 0.02$,\cite{nist:Lewetal2003} per each floor,
as compared to the resistive force of the intact
vertical columns, $f_{CC}+f_{PC} \simeq 0.8 + 2.7\cdot z$.
\item[2.] {\bf Rigidity assumption:}
The NIST report claims that the collapse started because the vertical columns
could not absorb the energy of the falling top section of the building.~\cite{nist:2005}
By design, all vertical columns were continuous structures that stretched
from the ground floor to the top of the building.
Lateral support was added to them to prevent them from buckling under
load, so that they would behave as ``short columns.''
For our models we assumed that the vertical columns are indeed short columns:
under compression they maintain their ultimate strength until the yield strain
is reached.
The rigidity assumption enters here as the location where the fracture occurs -
at (according to Bazant{\em et al.}\cite{bazant2006}) or near
the interface between the avalanche and the vertical column.
However, this is a slow compression of the column (the velocity of the source
of compression is much smaller than the sound velocity in the steel) so
the stress has time to propagate throughout the whole column causing
the strain to do the same.
As a result, the fractional distance $\lambda_1$ should be
applied to the full length of the column ($\sim H$, the height of the building)
and not to the storey height $\Delta H = H/F_T$.
Actual yield strain $\lambda_y$ can be estimated from the
yield force, $f(z) = r + s \, z$, where a local contribution to the
yield strain, $\delta\lambda_y$, is reversely proportional
to the local yield force, $\delta\lambda_y = C/f(z)$,
in accord with the uniform distribution of deformation energy.
This gives
\begin{equation}
  \label{eq:comment:2}
  \lambda_y = \sum \delta\lambda_y \simeq
  \int_{z_1}^1 \, dz \ \frac{C}{f(z)} =
  \lambda_1 \cdot \frac {r}{s} \cdot \log\frac{r+s}{r+s\,z_1}.
\end{equation}
Here, the limiting case of $s \rightarrow 0$ gives $C = \lambda_1 \, r$,
in which case $\lambda_y = \lambda_1 \cdot (1-z_1)$.
Calculation with intact parameters, $r=0.8$ and $s=2.7$, gives
${_1\lambda}_y = 0.063$ ($z_1 = 1 - F_1/F_T = 0.15$, for WTC 1 ) and
${_2\lambda}_y = 0.046$ ($z_1 = 0.3$, for WTC 2).
We leave as an exercise to the reader to show that these distances are sufficient
to stop the fall of the top section even if one makes a radical assumption that the
avalanche propagated through the primary zone without resistance ($r=s=0$).

Solution:
In the elastic case $(\epsilon\equiv1$) the velocities $v_1$ at the exit of the primary
zone are found from the energy conservation, ${_1v}_1 = 0.42$ and ${_2v}_1=0.5$
(for $\epsilon\equiv0$, one would use $P^2$ as an integral of motion, where $P$
is the momentum of the avalanche).
Let $\lambda\,H$ be the length needed to stop the fall, and let $r=0.8$ and
$s=2.7$ for the intact secondary zone, as before. Then, $\lambda$ is
given by
\begin{equation}
  \label{eq:comment:3}
  \lambda = \frac{z_1 \, v_1^2}{2 \,( f(z_1) - z_1)}.
\end{equation}
We find ${_1\lambda} = 0.013$ for WTC~1, and ${_2\lambda} = 0.029$ for WTC~2,
which are considerably smaller then their yield strains $\lambda_y$'s.
Thus, contrary to the NIST claim, the total plastic deformation of the
intact vertical columns in the secondary zone was more than sufficient to
arrest the fall of the top section.
\end{itemize}
From our discussion so far it follows that $(i)$, the secondary zones of both WTC 1 and
WTC 2 were not intact - in agreement with the hypothetical WMD destroying the
CCs discussed earlier, and $(ii)$, the destruction of the remaining vertical
columns (PCs) was either not through compression, or there
had to be a mechanism present that would pull the PCs inwards and into the path of the
avalanche.

This said, let us propose a consistent hypothetical model of an avalanche.
The avalanche is created by severing the core columns (CCs) at some distance
from the primary zone.
This makes the avalanche consist of the intact top section, the intact CCs of
which penetrate the secondary zone, and so give it an overall wedge-like appearance.
As a result of weight redistribution, the avalanche now interacts only with the
perimeter columns (PCs) from the top of the primary zone down to the level
at which the CCs were severed from the secondary zone CCs.
The avalanche's CCs pull the secondary zone PCs inwards, and so compromise them,
while the intact top section finishes the PCs as it goes down.
In this way, the avalanche never encounters the rigidity of the whole building,
just of its small section, as discussed earlier.
Furthermore, the pulling action is realized with the intact floor structure in the secondary
zone, through the tension of the floor trusses.
As is known, the tension yield of the floor trusses is much greater then
their shear yield force.
From the outside, it appears as if the avalanche starts at the weakest point of
the remaining structure: the compromised PCs in the primary zone.
By propagating so, the avalanche sees mostly the resistive force of the PCs in the
secondary zone, and some friction from the penetration of avalanche CCs into the
floor structure of the secondary zone.
The compromising of the secondary zone CCs continues so that the next severing point
is always ahead of the avalanche: otherwise, the avalanche's CCs might interfere
with the severing, which if prevented would result in a slowing down
of the avalanche.
The process continues until the avalanche reaches the ground floor.
We show the schematic of such collapse in \fig{fig:4}.
\begin{figure}[h]
  \singlespacing
  \centering
  \includegraphics[scale=0.5,clip]{fig4.eps}
  \caption{%
    \label{fig:4}
    Hypothetical model of progressive collapse in WTC 1 and 2, caused
    by compromising the central core columns (CCs) at predetermined
    locations along the height of the building and times, labeled
    ``the wave of massive destruction'' (WMD).
    As the core columns (CCs) are disabled, the collapse is opposed by the
    perimeter columns (PCs) mainly,  as discussed in text.
  }
\end{figure}

\newpage