\section{Introduction}
The higher derivative field theories are notorious for the stability
problems at interacting level, and also in quantum theory, see
\cite{PU, Pavsic, Woodard, Smilga} for discussions and further
references. The most frequently discussed higher derivative field
theories are the modified models of gravity where the stability is
an issue \cite{Salvio, Tomboulis, Belenchia}. Among the modified
gravity theories, $f(R)$ models (for review and further references
see \cite{Sotiriou, Felice}) provide the best known example of
stable nonlinear higher derivative field theory. Some other models
with similar properties can be found in the references
\cite{Smilga1, Bergshoeff2, Nitta}. The stability of this
exceptional class of higher derivative theories is related to the
fact that the canonical energy is bounded because of strong second
class constraints. For discussion of stability in various non-linear
higher derivative mechanical models we refer to \cite{Andrzejewski,
Smilga2, Chen, Pavsic1, Avendao-Camacho} and references therein.

In this paper, we consider a certain class of higher derivative
models which we call \textit{derived theories}. At free level, the
field equations of the derived theory are defined by the higher
derivative wave operator $M$ being a polynomial of another
differential operator $W$. The latter is supposed to be of the first
or second order. As we demonstrate, this class of systems admits,
under certain conditions, inclusion of stable interactions, and the
stability persists at quantum level.

Many well-known higher derivative models fall into the class of
derived theories. For example, the higher derivative scalar field of
the Pais-Uhlenbeck type \cite{PU} is a derived system, where $W$ is
the d'Alembert operator.  Podolsky electrodynamics \cite{Pod} is a
derived system, with the wave operator being a second-order
polynomial of Maxwell operator.  The extended Chern-Simons
\cite{Deser} is a derived theory of the vector field in $3d$
Minkowski space, with the wave operator being a third-order
polynomial in the Chern-Simons operator $*d$, where $*$ is the Hodge
star operator and $d$ is the de Rham differential. The Podolsky and
Chern-Simons electrodynamics have been discussed for many years from
various viewpoints, see \cite{Buffalo1, Buffalo3, Nogueira} and
references therein. In the conformal gravities in 4 and 6 dimensions
\cite{Pope1, Pope2} the linearized equations of motion for spin 2
fields belong to the derived type. The stability  is studied in
these works for the small fluctuations in the vicinity of constant
curvature backgrounds further extending earlier observation of the
work \cite{Maldacena}. With a special choice of boundary conditions,
it is observed that the theory might be stable. In the present work,
we follow a different idea being unrelated to the special boundary
conditions.

The simplest derived-type field equations correspond to the wave
operator being the second-order polynomial in another operator. This
case has been studied in \cite{KLS14}. In this class of models, it
turns out that at free level one can connect two different conserved
quantities to the time-shift symmetry. One of these quantities is
the canonical energy, while another one is a different independent
integral of motion. If the second quantity is bounded, the theory
will be stable at classical level, and the stability can be promoted
to quantum level \cite{KLS14}. In the paper \cite{KKL}, a more
general setup of the derived systems is studied when the wave
operator is a polynomial of arbitrary finite order $n$ in another
operator (see relation (\ref{derived}) in the present work). Once
the primary operator $W$ admits some symmetry, it can be connected
to $n$ conserved quantities of the derived system. In particular, if
one of the symmetries of $W$ is the homogeneity of time, then one
gets $n$ independent conserved quantities, which includes the
canonical energy. In this paper, we adopt slightly more general
setup for the derived systems aimed at relaxing restrictions on
inclusion of stable interactions. The fields are divided into
several subsets such that the wave operator in each subset is a
polynomial in a certain operator. Unlike a more simple case studied
in \cite{KKL}, the primary operators can be different for different
subsets of fields, and the wave operators are defined by different
polynomials for different fields. This more general setup at free
level provides more flexibility for inclusion of stable
interactions.

The stability of interactions in various particular types of higher
derivative theories of derived type has been previously studied in
\cite{KLS14, KKL, AKL18, KKL18}. The proper deformation method has
been suggested in \cite{KLS16} to systematically include stable
interactions in the derived theory which is stable at free level.
The crucial ingredient of the proper deformation method is the
Lagrange anchor. The Lagrange anchor has been first introduced in
the work \cite{KazLS} to BRST embed and quantize not necessarily
Lagrangian dynamics. Later, it has been found \cite{KLS10} that the
Lagrange anchor connects conserved quantities to symmetries for any
system of field equations, be they Lagrangian or not. For the
Lagrangian system the unit operator serves as the Lagrange anchor
that establishes one-tow-one correspondence between symmetries and
conserved currents. In principle, the Lagrange anchor is not
necessarily unique for given system of field equations. Once the
field equations admit multiple Lagrange anchors, the same symmetry
can be connected to different conserved quantities. As is noticed in
the paper \cite{KLS14}, even the simplest derived system, with the
wave operator $M$ being the second-order polynomial of another
differential operator $W$, admits two different Lagrange anchors.
This explains the existence of one more conserved quantity connected
to time-independence besides the canonical energy. Given the
Lagrange anchor, the general method of \cite{KLS16} allows one to
consistently include the interactions into field equations of motion
deforming the conserved quantities connected with the symmetries by
the anchor. If the anchor connects the symmetry with the bounded
quantity, the system remains stable upon inclusion of interaction by
the proper deformation scheme of \cite{KLS16}.

In this paper, we provide a more simple scheme for inclusion stable
interactions in a wider class of derived theories, skipping to
explicitly employ the Lagrange anchor. The basic idea is that we
have two subsets of fields such that the wave operator in each
subset is a polynomial in certain primary operator. The primary
operators are assumed to be Poincar\'e invariant. One of the primary
wave operators is supposed to be gauge invariant. We also assume
that the primary theories admit appropriate covariant interaction
vertex consistent with the gauge symmetry. To construct the stable
interactions between two derived theories, we identify the series of
conserved tensors at free level such that every item is connected to
the space-time translation symmetry. Some of the conserved tensors
can have bounded $00$-component, while the other ones are unbounded.
Then, we seek for the interactions in the derived theories such that
generalize couplings between primary models and keep the appropriate
bounded quantity conserved at interacting level. If the theory
admits bounded conserved quantity, the stability will persist at
interacting level once it is constructed by this method. The
proposed construction can reproduce all the perviously known stable
interactions in higher derivative systems of derived type
\cite{KLS14,KKL,AKL18} and we also add a new example in this paper,
to illustrate the method.

As the illustration of the general method described in the paper, we
construct the consistent and stable interaction between the order
$N$ extension of Chern-Simons and order $2n$ Pais-Uhlenbeck-type
higher derivative complex scalar field. These models have been
discussed earlier at free level in the papers \cite{KLS14, KKL,
AKL18, KKL18} for some specific orders $n$ and $N$. In this paper,
the $N+n$-parameter series of conserved second-rank tensors is
constructed for the free model being connected to the translation
invariance. Also $n$-parameter series of conserved currents is
identified to be connected to a single $U(1)$-symmetry of the
scalar. There are bounded quantities among the conserved observables
that make the theory stable at free level, while the canonical
energy is unbounded. Non-minimal gauge invariant couplings are
identified such that the higher derivative theory remains stable at
interacting level. The stable interaction vertices are inevitably
non-Lagrangian. That does not necessarily obstruct quantization as
we explain in the conclusion.

The article is organized as follows. In the next section, we define
the free derived theories and elaborate on the series of symmetries
and conserved quantities connected to each symmetry of the primary
wave operator $W$. In Section 3, we identify the gauge invariant
couplings such that provide conservation of any fixed representative
of the series of conserved quantities of the derived-type free
system. In Section 4, we consider interactions between the higher
derivative complex scalar field and higher derivative extension of
Chern-Simons. Following general pattern of the previous section, we
explicitly construct the conserved quantities and keep track of the
stability making use of the bounded integrals of motion. In
conclusion, we discuss the results.

\section{Derived-type theories, higher symmetries, and conservation laws}
Consider a set of fields $\phi^a$ in $d$-dimensional Minkowski space
with local coordinates $x^\mu,\mu=0,\ldots,d-1$. We systematically
use the DeWitt condensed notation \cite{DeWitt} in this paper. The
fields are labeled by the condensed indices $a,b, \ldots$. Every
condensed index accommodates all the vector, tensor, spinor,
isotopic, etc. indices and the space-time coordinates. Summation in
the condensed indices implies integration in $x^\mu$. In this
notation, the linear differential operators are represented by
matrices with condensed indices. We assume that the theory admits a
constant metric which is used to raise and lower multi-indices. In
this way, every linear operator has the quadric form. These
quadratic forms correspond to the local functionals which are
bilinear in the fields. We also imply that the fields vanish at the
infinity, so if the quadratic form vanishes in the condensed
notation, this means that the corresponding local functional is the
integral of total divergence.

In the condensed notation, any system of linear field equations
reads:
\begin{equation}\label{linear}
    M_{ab}\phi^b=0\,,
\end{equation}
where $M_{ab}$ is the integral kernel of matrix differential
operator. For the sake of simplicity, we assume that the matrix
$M_{ab}$ is square, so that the number of equations in each
space-time point coincides with the number of fields.  In this class
of theories, $M_{ab}$ is usually called the \emph{wave operator}.
The formal adjoint of the wave operator is defined by
\begin{equation}\label{}
    M^\dagger{}_{ab}=M_{ba}\,.
\end{equation}
The field equations are variational whenever $M^\dagger=M$, in which
case the action functional $S[\phi]$ reads
\begin{equation}\label{L-action}
    S[\phi]=\frac{1}{2}\langle\phi, M\phi\rangle\,,
\end{equation}
where the brackets $\langle\,\,,\,\rangle$ denote the natural pairing between
the fields,
\begin{equation}\label{}
    \langle\phi, M\phi\rangle\equiv\phi^a(M\phi)_a\,,\qquad (M\phi)_a\equiv M_{ab}\phi^b\,,
\end{equation}
and summation  is implied over the repeated multi-indices $a$ when
they stand at different levels. If $M^\dagger_\ast=-M_\ast$ for some
other operator $M_\ast$, its diagonal elements are total
divergencies,
\begin{equation}\label{div-j}
    \int d^dx \partial_\mu j^\mu(\phi)=\langle\phi,
    M_\ast\phi\rangle\,,\qquad \forall\phi\,.
\end{equation}
Once the expression $M_\ast\phi$ vanishes on-shell (\ref{linear}),
the latter formula establishes the relation between conserved
currents and anti-self-adjoint operators.

We say that the variational theory (\ref{linear}) is of the
\textit{derived}-type if the wave operator is a finite-order
polynomial of another self-adjoint operator $W{}^a{}_b$, i.e.
\begin{eqnarray}\label{derived}M(\alpha;W)=\sum_{p=0}^{n}\alpha_{p}W^p\,,\qquad W^\dagger=W\,,\qquad p=2,\ldots,n\,,
\end{eqnarray}
where all the multi-indices are raised and lowered by the metric,
and $\alpha_p\,,\,p=0,\ldots,n\,,$ are some real constants. The
order of the polynomial is assumed to be irreducible, so the
coefficient at the highest order is nonzero, $\alpha_n\neq0$. In
accordance with the definition, each derived theory is defined by
two ingredients: the self-adjoint operator $W{}^a{}_b$ and
finite-order polynomial,
\begin{equation}\label{Charpol}
    M(\alpha;z)=\sum_{p=0}^n\alpha_pz^p\,,
\end{equation}
with $z$ being formal complex-valued variable. We call $M(\alpha;z)$
the \emph{characteristic polynomial}, while $W{}^a{}_b$ is referred
to as the \emph{primary wave operator}. Being considered in itself,
the primary operator defines the \emph{primary free field theory},
\begin{equation}\label{primary}
    W_{ab}\phi^b=0\,.
\end{equation}
As $W_{ab}$ is self-adjoint, the primary theory is variational. In
our article, we mostly deal with the class of theories, where the
primary wave operator does not involve higher derivatives. In this
setting, the higher derivative derived model (\ref{linear}),
(\ref{derived}) can be considered as originating from the lower
order primary theory (\ref{primary}).

In this paper, we consider the linear operator $X{}^a{}_b$ as a
symmetry\footnote{The notion of symmetry can be understood in
various ways. In this section, we provide the simple non-rigorous
understanding, which is sufficient for the purposes of this work.}
of linear equations (\ref{linear}) if it is interchangeable with the
wave operator of theory in the following sense:
\begin{eqnarray}\label{Sym}
    [X,M]{}^a{}_b=Y^a{}_cM{}^c{}_b\,,\qquad [X,M]{}^a{}_b=X{}^a{}_cM{}^c{}_b-M{}^a{}_c X{}^c{}_b\,,
\end{eqnarray}
with $Y{}^a{}_c$ being some other operator. Once this relation is
valid, the operator $X{}^a{}_b$ defines the linear transformation of
the fields such that leaves the mass shell (\ref{linear}) intact,
\begin{eqnarray}\label{}
    \delta_\xi
    (M\phi)^a=\xi(X+Y){}^a{}_b (M\phi)^b\approx0\,,\qquad \delta_\xi\phi^a=\xi X{}^a{}_b\phi^b\,,
\end{eqnarray}
where $\xi$ is the transformation parameter, being some constant.
The sign $\approx$ means equality modulo equations (\ref{linear}).
Each theory admits the trivial symmetries that read
\begin{equation}\label{}
    X{}^{a}{}_b=\widetilde{X}^a{}_cM{}^c{}_b\,,\qquad
    Y{}^a{}_b=[\widetilde{X},M]{}^a{}_b\,,
\end{equation}
where $\widetilde{X}^a{}_b$ can be any \textbf{}operator. The
corresponding transformations of fields vanish on-shell and do not
contain any valuable information about the dynamics of model. We
systematically ignore the trivial symmetries, considering any
symmetry modulo trivial one.

The non-trivial symmetries of linear system are known to form an
associative algebra. The multiplication operation is just a
composition of operators. It is easy to verify the latter fact: for
each pair of symmetries $(X_1){}^a{}_b,(X_2){}^a{}_b$, we get
\begin{eqnarray}\label{}
    &&[X_1X_2,M]{}^a{}_b=(X_1){}^a{}_c[X_2,M]{}^c{}_b+
    [X_1,M]{}^a{}_c(X_2){}^c{}_b\notag\\
    &&=((Y_1){}^a{}_d(X_2)^d{}_c+(X_1){}^a{}_d(Y_2){}^d{}_c+(Y_2)^a{}_c)
    M{}^{c}{}_b\,.
\end{eqnarray}
The elements of associative algebra that commute with the wave
operator form a subalgebra. In this paper, we mostly consider this
subalgebra as it is connected to conserved quantities.

Now, consider the derived equations (\ref{linear}), (\ref{derived}).
Assume that the primary model (\ref{primary}) has the symmetry
$X{}^a{}_b$ such that commutes with the primary wave operator,
\begin{equation}\label{W-Sym}
    [X,W]=0\,.
\end{equation}
Even if the primary theory admits a single symmetry, the associative
algebra of symmetries of derived theory has the two natural
generators -- $X$ and $W$ -- that commute with each other,
\begin{equation}\label{}
    [X,M]=[W,M]=0\,.
\end{equation}
 By composing the primary symmetry $X$ with the degree of the primary wave operator $W$
 we get the symmetry of the derived equations (\ref{derived}):
\begin{eqnarray}\label{Xp-Sym}
   [X_p,M]=0\,,\qquad (X_p){}^a{}_b=X{}^a{}_c(W^p){}^c{}_b\,,\quad p=0,\ldots,n-1\,.
\end{eqnarray}
Only the terms with $p=0,\ldots,n-1$ are relevant, because the
higher powers of primary operator can be absorbed by the wave
operator (\ref{derived}). In this way, for $p\geq n$, the symmetry
reduces on-shell to the symmetry with $p<n$. The generators $X_p$
can be assembled into the $n$-parameter series of \emph{derived
symmetries} of derived model,
\begin{equation}\label{XSer-Sym}
    X{}^a{}_b(\beta)=\sum_{p=0}^{n-1}\beta_p(X_p){}^a{}_b\,,
\end{equation}
with $\beta_p\,,\,p=0,\ldots,n-1\,,$ being real numbers. As is seen,
all the representatives of the series originate from one and the
same symmetry $X$ of primary model.

The symmetry (\ref{Sym}) preserves the action functional
(\ref{L-action}) if its operator is anti-self-adjoint and commutes
with the wave operator,
\begin{equation}\label{}
    [X,M]=0\,,\qquad X^\dagger=-X\,.
\end{equation}
The corresponding conserved current $j^\mu(\phi)$, being quadratic
in the fields, is defined by the condition
\begin{equation}\label{j-char}
    \int d^dx \partial_\mu j^\mu(\phi)=\langle\phi,XM\phi\rangle\,.
\end{equation}
This formula represents the Noetherian relationship between
symmetries and conservation laws. In the class of derived theories,
a single symmetry (\ref{W-Sym}) of primary model (\ref{primary})
determines the series of derived symmetries (\ref{XSer-Sym}). The
associated $n$-parameter series of conserved currents reads
\begin{eqnarray}\label{jSer}
    j^{\mu}(\beta)=\sum_{p=0}^{n-1}\beta_pj_p{}^\mu(\phi)\,,\qquad \int d^dx \partial_\mu
j_p{}^\mu(\phi)=\langle\phi,X_pM\phi\rangle\,,
\end{eqnarray}
where the quantities $j_p{}^\mu, p=0,\ldots,n-1$\,, come from
symmetries (\ref{Xp-Sym}). In particular, $j_0{}^\mu$ represents the
Noether conserved current for the symmetry (\ref{W-Sym}) of primary
model, while $j_p{}^\mu,p=1,\ldots,n-1\,,$ are other quantities.

The simplest possible symmetry of the free field theory is the
translation invariance. The translation generators $\partial_\mu$
are automatically anti-self-adjoint and they commute with the
primary wave operator (\ref{primary}),
\begin{equation}\label{}
    [\partial^\mu,W]=0\,.
\end{equation}
Once the primary wave operator is translation invariant, the derived
theory enjoy $n$-parameter series (\ref{XSer-Sym}) of derived
symmetries originating from this invariance:
\begin{eqnarray}\label{}
    X{}^{\mu a}{}_{b}(\beta)=\sum_{p=0}^{n-1}\beta_p(X_p){}^{\mu a}{}_{b}\,,\qquad (X_p){}^{\mu a}{}_{b}=(\partial^\mu W^p){}^a{}_b\,.
\end{eqnarray}
Each of these symmetries preserves the action (\ref{L-action}) of
the derived model (\ref{linear}), (\ref{derived}). The corresponding
conserved currents (\ref{jSer}) constitute the series of second-rank
energy-momentum tensors
\begin{eqnarray}\label{tSer}
    \mathit\Theta^{\mu\nu}(\beta)=\sum_{p=0}^{n-1}\beta_p
    T_p{}^{\mu\nu}(\phi)\,,\qquad \partial_\nu T_p{}^{\mu\nu}(\phi)=\langle \phi, \partial^\mu W^p
    M\phi\rangle\,.
\end{eqnarray}
This series includes the canonical energy-momentum as
$T_0{}^{\mu\nu}(\phi)$, while $T_p{}^{\mu\nu}, p=1,\ldots,n-1\,,$
are different conserved tensors associated with space-time
translation invariance of the model. In this way, the translation
invariance of the derived-type field theory results in the series of
the conserved tensors.

As the stability of higher derivative model is concerned, the
$00$-component of the conserved tensor is of interest. The
$00$-component of tensor (\ref{tSer}) reads
\begin{equation}\label{t00-Ser}
    \mathit\Theta^{00}(\beta)=\sum_{p=0}^{n-1}\beta_p
    T_p{}^{00}(\phi)\,.
\end{equation}
This expression is given by the sum of canonical energy $T_0{}^{00}$
and the other contributions $T_p{}^{00},p=1,\ldots,n-1$. Even if the
canonical energy is unbounded due to higher derivatives, the
quantity (\ref{t00-Ser}) can define bounded conserved charge. The
bounded conserved quantity, if it exists, stabilizes classical
dynamics of derived model (\ref{linear}), (\ref{derived}) at free
level.

The entries of conserved current series (\ref{jSer}) (or conserved
tensor series (\ref{tSer})) are independent in general, even though
it is not a theorem. This fact is supported by the following
observations: (i) the quantities $j_p{}^\mu $ are bilinear forms in
the fields $\phi^a$ and their space-time derivatives, and (ii) the
total number of derivatives involved in $j_p{}^\mu $ increases with
$p$. Once the symmetry $X$ and primary operator $W$ are the matrix
differential operators of orders $n_X$ and $n_W$, the conserved
current $j_p{}^\mu$ involves at most
\begin{equation}\label{}
    (p+n)n_W+n_X-1
\end{equation}
derivatives of fields. The term with the highest number of
derivatives contributes to the conserved quantity, and it is
on-shell nontrivial at least in the case without gauge symmetry and
constraints. Applying this argument to $j_{n-1}{}^\mu$, we conclude
that the highest order term can't come from the linear combination
of other currents with the lower number of derivatives. Thus,
$j_{n-1}{}^\mu$ is an independent conserved quantity. Proceeding
with the same argument to
 $j_{n-2}{}^\mu, j_{n-3}{}^\mu,\ldots, j_0{}^\mu$, we conclude that all these
quantities are not functions of each other. The maximal number of
independent entries in the conserved current series (\ref{jSer})
equals to the order $n$ of characteristic polynomial
(\ref{Charpol}). In the higher derivative scalar field model the
maximal number of independent conserved quantities is connected with
the space-time translation invariance\cite{KLS14}. In gauge
theories, some generators of series (\ref{jSer}) can be on-shell
trivial. For example, the extended Chern-Simons model of order $n$
admits $n-1$ independent conserved tensors (\ref{tSer})\cite{KKL}.
All the above applies to free theories. At the non-linear level, the
procedure of inclusion of interaction is sensitive to the choice of
representative in the conserved tensor series (\ref{jSer}) which is
conserved at the interacting level. This motivates us to work with
the full series of conserved quantities, even if some of its
generators are dependent.

The discussion above does not address the gauge symmetries. Now we
explain how the gauge symmetry is accounted for in the field
theories of derived type. The wave operator may have a null-vector,
\begin{equation}\label{WA-null}
    M{}^a{}_bR{}^b{}_\alpha=0\,.
\end{equation}
In this case, the derived theory  (\ref{linear}), (\ref{derived}) admits the gauge
transformation such that preserves the mass shell,
\begin{equation}\label{gs-primary}
    \delta_\varepsilon(M{}_{ab}\phi^b)=0\,,\qquad\delta_\varepsilon\phi^a=R{}^a{}_\alpha\varepsilon^\alpha,
\end{equation}
where $\varepsilon^{\alpha}$ are the gauge parameters, being
arbitrary functions of space-time coordinates. The operator
$R{}^a{}_\alpha$ is the gauge symmetry generator. In the present
work, we consider the class of derived models whose gauge symmetries
come from null-vectors of primary model, i.e.
\begin{equation}\label{}
     M{}^a{}_bR{}^b{}_\alpha=0\qquad\Leftrightarrow\qquad W{}^a{}_bR{}^b{}_\alpha=0\,.
\end{equation}
This property is not automatically satisfied for general derived
theory, so this is an additional assumption. The Podolsky's
electrodynamics and the extended Chern-Simons are consistent with
this assumption, for example.

\section{Consistent interactions of derived models}
In this section, we construct the class of consistent interactions
between two Poincar\'e invariant derived theories, with one of them
being gauge. At free level, these theories admit series of conserved
energy-momentum tensors. We construct the interaction such that
provides conservation of the deformation of certain representative
of energy-momentum tensor series. In so doing, we admit not
necessarily Lagrangian vertices.\footnote{For general problems of
inclusion of consistent interactions between not necessarily
Lagrangian theories see \cite{KLS13}.} To solve the problem we
proceed from the assumption that the two primary theories admit
consistent interaction vertex. Then, we seek for the way to lift
this vertex to the level of derived theory in such a way that does
not spoil the stability of the free derived theory.

Consider two subsets of fields $\phi^a\,,\mathit{\Phi}^i$ on
$d$-dimensional Minkowski space, with the multi-indices $a,i$
labeling the fields. The primary operators are denoted by $W_{ab}$
and $\mathcal{W}_{ij}$, so the primary field equations read
\begin{equation}\label{2-primary}
    W_{ab}\phi^b=0\,,\qquad \mathcal{W}_{ij}\mathit{\Phi}^j=0\,.
\end{equation}
These equations are variational, with the action functional being
the sum of actions (\ref{L-action}) of the fields $\phi^a$ and
$\mathit{\Phi}^i$. The primary operators are assumed to be invariant
under the space-time translations,
\begin{equation}\label{}
    [\partial^\mu,W]=[\partial^\mu,\mathcal{W}]=0\,.
\end{equation}
The Poincar\'e symmetry implies the existence of conserved
energy-momentum tensor
\begin{equation}\label{thetaphiPhi}
    \mathit{\Theta}^{\mu\nu}(\phi,\mathit{\Phi})=T^{\mu\nu}(\phi)+T^{\mu\nu}(\mathit{\Phi})\,,
\end{equation}
where $T^{\mu\nu}(\phi)$ and $T^{\mu\nu}(\mathit{\Phi})$ denote the
contributions of free fields $\phi$ and $\mathit{\Phi}$. The
quantity $\mathit{\Theta}^{\mu\nu}$ is determined by the formula
\begin{equation}\label{}
    \int d^dx\partial_\nu\mathit{\Theta}^{\mu\nu}(\phi,\mathit\Phi)=\langle\phi, \partial^\mu W\phi\rangle+
    \langle\mathit\Phi,\partial^\mu\mathcal{W}\mathit\Phi\rangle\,.
\end{equation}
The theory of the field $\mathit{\Phi}^i$ is assumed to be
gauge invariant. The gauge symmetry generator $R^i{}_\alpha$ of the
field is the null-vector of the primary operator $\mathcal{W}_{ij}$,
\begin{equation}\label{mathcalWR}
    \mathcal{W}_{ij}R^{j}{}_\alpha=0\,,
\end{equation}
see Eq. (\ref{WA-null}). The corresponding gauge transformation reads
\begin{equation}\label{gt-P}
    \delta_\varepsilon\mathit{\Phi}^i=R^i{}_\alpha\varepsilon^\alpha,\qquad \delta_\varepsilon\phi^a=0\,,
\end{equation}
where $\varepsilon^\alpha$ are gauge parameters. We do not allow gauge freedom for the
fields $\phi^a$.

We  assume that there exists a variational interaction vertex
between primary models such that (i) the gauge transformation
(\ref{gt-P}) is preserved at the interacting level and (ii) the
action functional is at most quadratic in $\phi$. The most general
action  such that meets these requirements has the form
\begin{equation}\label{S-primary-int}
    S[\phi,\mathit{\Phi}]=\frac{1}{2}(W_{ab}+\mathit{\Gamma}_{ab}(\beta \mathit{\Phi}))\phi^a\phi^b+\frac{1}{2}\mathcal{W}_{ij}\mathit{\Phi}^i\mathit{\Phi}^j\,,
\end{equation}
where the operator of vertex
$\mathit{\Gamma}_{ab}(\beta\mathit{\Phi})$ is a function of quantity
$\beta\mathit{\Phi}$, and $\beta$ is the coupling constant, which
can be arbitrary real number. We assume that
$\mathit{\Gamma}_{ab}(\beta\mathit{\Phi})$ is a finite-order
polynomial in $\mathit\Phi^i$, i.e.
\begin{equation}\label{}
    \mathit\Gamma_{ab}(\beta
    \mathit\Phi)=\sum_{k=1}^{k_{\text{max}}}\beta^k\mathit{\Gamma}^{(k)}{}_{ab}(\mathit{\Phi})\,,\qquad
    k_{\text{max}}<+\infty\,,
\end{equation}
with $\mathit{\Gamma}{}^{(1)}{}_{ab},
\mathit{\Gamma}^{(2)}{}_{ab},\ldots$ being linear, quadratic etc. in
the field $\mathit{\Phi}$. The expression $\beta^k$ means the $k$-th power of $\beta$. The presence of such decomposition
ensures that the equations of motion of the model are finite-order
polynomials in the coupling constant $\beta$. The Lagrange equations for the action functional
(\ref{S-primary-int}) read
\begin{eqnarray}\label{EL-primary-int}
   \partial_a S\equiv W_{ab}(\beta\mathit\Phi)\phi^b=0\,,\qquad \partial_i S\equiv \mathcal{W}_{ij}\mathit{\Phi}^j+\beta J_i(\phi,\beta\mathit\Phi)=0\,,
\end{eqnarray}
where
\begin{eqnarray}\label{W-inv}
W_{ab}(\beta\mathit\Phi)=W_{ab}+\mathit{\Gamma}_{ab}(\beta\mathit{\Phi})\,,\qquad
J_{i}(\phi,\beta\mathit\Phi)=
\partial_i\mathit\Gamma_{ab}(\beta\mathit\Phi)\phi^a\phi^b\,,
\end{eqnarray}
and the  derivatives $\partial_a, \partial_i$ are understood as
variational in $\phi^a,\Phi^i$. As is seen, in the model with
interaction, the operator of vertex
$\mathit\Gamma_{ab}(\beta\mathit\Phi)$ is added to free equations of
originally non-gauge field $\phi^a$, while the gauge field theory
get additional current-like contribution
$J_i(\phi,\beta\mathit\Phi)$. This structure of non-linear theory is
typical for lower spin fields. For example, it includes the minimal
coupling between the electromagnetic field and the charged matter
fields.

The infinitesimal gauge transformation of the action functional
(\ref{S-primary-int}) reads
\begin{equation}\label{gt-pA}
    \delta_\varepsilon\phi^a=\beta R{}^{a}{}_{b\alpha}\phi^b\varepsilon^\alpha\,,\qquad
    \delta_\varepsilon \mathit\Phi^i=R{}^{i}{}_{\alpha}\varepsilon^\alpha\,,
\end{equation}
where $R{}^i{}_\alpha$ is the null-vector of primary operator $\mathcal{W}_{ij}$
(\ref{mathcalWR}), and $R^a{}_{b\alpha}$ is some field-independent
structure function. We assume that it is antisymmetric in the
indices $ab$,
\begin{equation}\label{R-R}
 R{}_{ab\alpha}=-R{}_{ba\alpha}\,.
\end{equation}
The invariance of action (\ref{S-primary-int}) with respect to the
gauge transformations (\ref{gt-pA}) is equivalent to the condition
\begin{eqnarray}\label{G-W-R}
     (W+\mathit{\Gamma}(\beta\mathit{\Phi}))^{a}{}_cR^c{}_{b\alpha}-R^a{}_{c\alpha}(W+\mathit{\Gamma}(\beta\mathit{\Phi}))^c{}_b+\,R^i{}_\alpha\partial_i\mathit{\Gamma}{}^a{}_b(\beta\mathit{\Phi})=0\,.
\end{eqnarray}

As the gauge symmetries and gauge identities of Lagrangian theories
are connected by the Noether theorem, conditions (\ref{R-R}), (\ref{G-W-R}) determine the
selection rule for the structure functions
$\mathit{\Gamma}{}^{a}{}_{b},R^{a}{}_{b\alpha}$. The meaning of these conditions
is that the wave operator $W_{ab}(\beta\mathit\Phi)$ (\ref{W-inv}) is gauge invariant in the following sense:
\begin{eqnarray}\label{[R,W]}
    [W(\beta\mathit\Phi), R_\alpha]_{ab}&+&R^i{}_\alpha\partial_i\mathit{\Gamma}{}_{ab}(\beta\mathit{\Phi})=0\,,\notag\\\,
    [W(\beta\mathit\Phi), R_\alpha]_{ab}&=&W{}_{ac}(\beta\mathit\Phi)R^{c}{}_{b\alpha}-R_{ac\alpha}W{}^c{}_b(\beta\mathit\Phi)\,,\qquad\forall\beta\,.
\end{eqnarray}
The latter relation immediately ensures the gauge invariance of equations of motion of the field $\phi^a$.

The translation invariance of interacting theory
(\ref{S-primary-int}) is understood in the following sense:
\begin{eqnarray}\label{[p,W]}
   [W(\beta\mathit\Phi),\partial^\mu]{}_{ab}&+&\beta\partial^\mu\mathit\Phi^i\partial{}_i\mathit\Gamma{}_{ab}(\beta\mathit\Phi)=0\,,\notag\\\,
   [W(\beta\mathit\Phi), \partial^\mu]_{ab}&=&W{}_{ab}(\beta\mathit\Phi)\partial^\mu-\partial^\mu W{}_{ab}(\beta\mathit\Phi)\,,\qquad \forall\beta\,.
\end{eqnarray}
Once this condition is met, the action functional
(\ref{S-primary-int}) is preserved by the space-time translations.
The energy-momentum tensor of the non-linear theory has the
following structure:
\begin{equation}\label{ThetaSer}
    \boldsymbol{\mathit{\Theta}}^{\mu\nu}(\phi,\mathit{\Phi})=T^{\mu\nu}(\phi,\mathit\Phi)+T^{\mu\nu}(\mathit{\Phi})\,.
\end{equation}
The expression $T^{\mu\nu}(\phi,\mathit\Phi)$ is the gauge invariant
extension of the tensor $T^{\mu\nu}(\phi)$ (\ref{thetaphiPhi}) of
free theory, while $T^{\mu\nu}(\mathit{\Phi})$ is the
energy-momentum of free gauge field. The defining condition for the
conserved tensor reads
\begin{equation}\label{}
    \int d^d x\partial_\nu\boldsymbol{\mathit{\Theta}}^{\mu\nu}(\phi,\mathit\Phi)=
    \langle\phi, \partial^\mu W(\beta\mathit\Phi)\phi\rangle+
    \langle\mathit\Phi,\partial^\mu\mathcal{W}\mathit\Phi\rangle\,.
\end{equation}
The rhs of this expression is a total divergence because
$\partial^\mu W(\beta\mathit\Phi), \partial^\mu\mathcal{W}$ are
anti-self-adjoint operators.

Let us explain the meaning of conditions (\ref{R-R}), (\ref{G-W-R})
and its impact on the structure of interactions. Assume that the
gauge transformation of the gauge field $\mathit{\Phi}^i$ has a
null-mode $\varepsilon=\overline{\varepsilon}(x)$ such that
\begin{equation}\label{}
    R{}^i{}_\alpha\overline{\varepsilon}^\alpha=0\,.
\end{equation}
In this case, the action functional (\ref{S-primary-int}) is
preserved by the following transformation:
\begin{equation}\label{rigid-sym}
    \delta_{\overline{\varepsilon}}\phi^a=\xi
    R{}^a{}_{b\alpha}\overline{\varepsilon}^\alpha\,,\qquad
    \delta_{\overline\varepsilon}\mathit\Phi^i=0\,,
\end{equation}
with the parameter $\xi$ being a constant. In the free limit, this
transformation corresponds to certain internal symmetry of the free
model of the fields $\phi^a$, which localizes at the interacting
level. This symmetry does not follow from the relativistic
invariance of equations (\ref{2-primary}), so we have additional
prerequisite for constriction of interaction at free level. As we
deal with internal symmetry, this gives a restriction on the
multiplet of fields $\phi^a$ that are involved into the interaction.

The current-like term $J_i$ (\ref{W-inv}) originates from the
internal symmetry (\ref{rigid-sym}). As a consequence of relations
(\ref{R-R}), (\ref{G-W-R}), it is gauge invariant and satisfies
gauge identity,
\begin{equation}\label{R-J}
    \delta_\varepsilon J_i=0\,,\qquad R^i{}_\alpha J_i+R{}^a{}_{b\alpha}\partial_a S\, \phi^b=0\,.
\end{equation}
If the internal symmetry is the $U(1)$-transformation, and gauge
generator $R{}^i{}_\alpha$ is gradient, the current-like term
literally corresponds to the current of the $U(1)$-charge,
\begin{equation}\label{}
    J_i=J_\mu(x)\,,
\end{equation}
which meets condition (\ref{R-J}). Once current-like term $J_i$ is
added to the equations of motion of the gauge field, and the primary
operator $W$ is replaced by its gauge invariant extension
$W(\beta\mathit\Phi)$\,, the non-linear theory remains gauge invariant.
This means that interactions (\ref{EL-primary-int}) imitate the
electromagnetic-like couplings between gauge and matter fields.

Let us now turn to the details of inclusion of interactions between
derived theories. The most general ansatz for two derived theories
with the primary operators $W$ and $\mathcal{W}$ has the form
\begin{equation}\label{2-derived}
    M_{ab}(\alpha;W)\phi^b=0\,,\qquad \mathcal{M}{}_{ij}(A;\mathcal{W})\mathit{\Phi}^j=0\,.
\end{equation}
The characteristic polynomials of derived models (\ref{2-derived})
are the most general of orders $n$ and $N$, respectively,
\begin{equation}\label{2-Charpol}
    M(\alpha;z)=\sum_{p=0}^{n}\alpha_{p}z^p\,,
    \qquad
    \mathcal{M}(A;z)=\sum_{q=1}^{N}A_qz^q\,,
\end{equation}
where $z$ is a formal complex-valued variable. The real numbers
$\alpha_p\,,\, p=0,\ldots,n\,,$ and $A_q\,,\,q=1,\ldots,N\,,$ are model
parameters that distinguish different theories in
the considered class.

The derived theories (\ref{2-derived}) have the same space-time and gauge
symmetries as primary models. In particular, models (\ref{2-derived}), (\ref{2-Charpol}) are invariant under
gauge transformation (\ref{gt-P}). As for space-time translation,
a single symmetry of primary theory induces the $n+N$-parameter
series of derived symmetries such that
\begin{eqnarray}\label{}
    X^{\mu a}{}_b(\beta)=\sum_{p=0}^{n-1}\beta_p(\partial^\mu W^p){}^a{}_{b}\,,\qquad X^{\mu
i}{}_j(B)=\sum_{q=0}^{N-1}B_q(\partial^\mu\mathcal{W}^q){}^i{}_j\,,
\end{eqnarray}
where the real numbers $\beta_p,B_q$ are parameters. The
corresponding set of the conserved quantities is given by the series
of energy-momentum tensors
\begin{equation}\label{Theta1}
    \mathit{\Theta}^{\mu\nu}(\beta,B)=\sum_{p=0}^{n-1}\beta{}_p T{}_p{}^{\mu\nu}(\phi)+\sum_{q=0}^{N-1}B{}_qT{}_q{}^{\mu\nu}(\mathit{\Phi})\,,
\end{equation}
where the quantities $T{}_p{}^{\mu\nu}(\phi)$, $T{}_q{}^{\mu\nu}(\mathit\Phi)$ are
defined by the conditions
\begin{eqnarray}\label{d-Theta-free}
    \int d^d x\partial_\nu\mathit{\Theta}^{\mu\nu}(\beta,B)&=&\sum_{p=0}^{n-1}\beta{}_p \langle\phi,\partial^\mu W^pM\phi\rangle+\sum_{q=0}^{N-1}B{}_q
\mathit\langle\mathit\Phi,\partial^\mu\mathcal{W}^q\mathcal{M}\mathit\Phi\rangle\,.
\end{eqnarray}
By construction, the canonical energy-momentum
tensor of the model is included in the series as
\begin{equation}\label{}
\mathit{\Theta}^{\mu\nu}_{\text{can}}(\phi,\mathit{\Phi})=
T{}_0{}^{\mu\nu}(\phi)+T{}_0{}^{\mu\nu}(\mathit{\Phi})\,.
\end{equation}
Even though the canonical energy-momentum is almost always unbounded
due to higher derivatives, the other bounded quantities can present
in the series. Free higher derivative theory can be stabilized
by these quantities.

By interaction between two gauge models (\ref{2-derived}) we mean a
formal deformation of free theory with coupling parameter $\beta$
such that the non-linear equations of motion read
\begin{eqnarray}\label{E}
    &&E_{a}(\beta)\equiv M_{ab}\phi^b+\sum_{k=1}^\infty\beta^k
    \mathit{\Gamma}{}^{(k)}{}_{a}(\phi,\mathit{\Phi})=0\,,\notag\\
    &&\mathcal{E}_{i}(\beta)\equiv \mathcal{M}_{ij}
    \mathit{\Phi}^j+\sum_{k=1}^{\infty}\beta^k\mathit{\Gamma}^{(k)}{}_{i}(\phi,\mathit{\Phi})=0\,.
\end{eqnarray}
The structure functions $\mathit{\Gamma}{}^{(k)}$ are assumed to
have the homogeneity degree $k+1$ in the variables
$\phi^a,\mathit{\Phi}^i$, so $\mathit{\Gamma}{}^{(1)}$ are quadratic
in the fields, $\mathit{\Gamma}{}^{(2)}$ are cubic, and etc. No
action principle is required for interacting theory, so non-linear
equations (\ref{E}) can be non-Lagragian. We consider interaction
(\ref{E}) consistent if the number of gauge symmetries and gauge
identities is preserved by coupling and at least one representative
of the free series of energy-momentum tensors (\ref{d-Theta-free})
is still conserved at the non-linear level.

The $n+N$-parameter series of interaction vertices for free
equations (\ref{2-derived})  can be presented in the following form
\begin{eqnarray}\label{e-int-general}
    E_a\equiv\overline{M}{}_{ab}(\beta,B)\phi^b=0\,,\qquad\mathcal{E}_i\equiv
\mathcal{M}_{ij}\mathit{\Phi}^j+\sum_{p=0}^{n-1}\beta_p(\overline{J}_{p}){}_{i}(\beta,B)=0\,,
\end{eqnarray}
where the notation is used:
\begin{eqnarray}\label{Wbar}
    \overline{M}_{ab}(\beta,
    B)=\sum_{p=0}^{n}\alpha{}_{p}\overline{W}^p(\beta,B)\,,\qquad
    (\overline{J}_{p}){}_{i}(\beta, B)=\frac{1}{2}\partial_i\langle\phi,\overline{W{}}^p(\beta,B)\overline{M}(\beta,B)\phi\rangle\,,
\end{eqnarray}
and the line above means that the gauge field $\mathit\Phi^i$ is replaced by more general expression:
\begin{eqnarray}\label{Pbar}
    \overline{W}(\beta,B)=\left.W(\beta\mathit\Phi)\right|_{\beta\mathit\Phi=\overline{\mathit\Phi}(\beta,B)}\,,\qquad \overline{\mathit{\Phi}}{}^i(\beta,B)=\beta\mathit\Phi^i+\sum_{q=1}^{N-1}B_{q}(\mathcal{W}^q\mathit{\Phi}){}^i\,.
\end{eqnarray}
The coupling constants are real numbers
$\beta\,,\,\beta_{p}\,,\,B_q\,,\,p=0,\ldots,n-1\,,\,q=1,\ldots,N-1\,$.
Interactions (\ref{e-int-general}) generalize couplings of the
primary theory (\ref{EL-primary-int}) in the following sense: the
free wave operator $M$ of the fields $\phi^a$ is replaced by its
gauge invariant extension $\overline{M}(\beta,B)$, and the series of
current like terms $(\overline{J}_{p}){}_{i}(\beta, B)$ is added to
the free gauge equations.

The following facts ensure consistency of interaction (\ref{e-int-general}):
(i) the equations of motion are preserved by the gauge transformation
 (\ref{gt-pA}), while transformation law for the equations of motion reads
\begin{equation}\label{gt-eom}
    \delta_\varepsilon E^{a}=\beta R^{a}{}_{b\alpha}E^b\varepsilon^\alpha\,,\qquad \delta_\varepsilon \mathcal{E}^i=0\,;
\end{equation}
(ii) there are the gauge identities between equations of motion,
\begin{equation}\label{g-id}
    R{}^i{}_\alpha
    \mathcal{E}_i+\sum_{p=0}^{n-1}\beta_pR_{ab\alpha}(\overline{W}{}^p){}^b{}_cE^a\phi^c\equiv0\,;
\end{equation}
(iii) the second-rank conserved tensor $\boldsymbol{\mathit\Theta}^{\mu\nu}$ is defined by the condition
\begin{eqnarray}\label{d-Theta}
    \int d^dx \partial_\nu\boldsymbol{\mathit\Theta}^{\mu\nu}(\beta,B)=\sum_{p=0}^{n-1}\beta_p
    \langle\phi,\partial^\mu\overline{W}{}^pE\rangle
    +\sum_{q=0}^{N-1}B_q \langle\mathit\Phi, \partial^\mu \mathcal{W}^q \mathcal{E}\rangle\,,\qquad B_0\equiv\beta\,.
\end{eqnarray}
The proof of relations (\ref{gt-eom}), (\ref{g-id})  uses identities
\begin{equation}\label{App-ids}\begin{array}{c}
    (\text{a})\,\delta_\varepsilon (\overline{W}{}^p\phi)^a=\beta R^{a}{}_{b\alpha}(\overline{W}{}^p\phi)^b\varepsilon^\alpha\,;
    \quad(\text{b})\,
    R^i{}_\alpha(\overline{J}_{p}){}_{i}+R_{ab\alpha}(\overline{W}{}^p){}^b{}_cE^a\phi^c\equiv0\,;\\[3mm](\text{c})\,
    \delta_\varepsilon (\overline{J}_{p}){}_{i}=0\,,
\end{array}\end{equation} which hold for $p=0,\ldots,n-1$. We deduce these
relations in the Appendix A. The conserved tensor
$\boldsymbol{\mathit\Theta}^{\mu\nu}$ (\ref{d-Theta}) is a
deformation of a selected representative of free energy-momentum
series (\ref{d-Theta-free}), whose parameters are defined by
coupling constants. In Appendix B, we ensure the rhs of Eq.
(\ref{d-Theta}) is a total divergence and establish the structure of
conserved quantity.

Let us comment on the origin of interactions (\ref{e-int-general}).
The model (\ref{2-derived}) admits $n$-parameter series of gauge
invariant current-like terms
\begin{eqnarray}\label{J-1}
    J_i(\beta)=\sum_{p=0}^{n-1}\beta_p(J_{p}){}_{i}(\beta)\,,\qquad (J_{p}){}_{i}(\beta)=\frac{1}{2}\partial_i\langle\phi,W{}^p(\beta\mathit{\Phi})M(\beta\mathit{\Phi})\phi\rangle\,.
\end{eqnarray}
Once the primary wave operator $W$ is replaced by its gauge
invariant extension $W(\beta\mathit\Phi)$, the inclusion of such
current-like terms preserves gauge invariance. This means that
theory (\ref{e-int-general}) with the gauge-invariant primary
operator (\ref{W-inv}) and current-like term (\ref{J-1}) is
consistent. At this ways, we conclude that there is $n+1$-parameter
series of interactions with coupling constants
$\beta,\,\beta_0,\ldots,\beta_{n-1}\,.$ To involve coupling
constants $B_q$ into account, we perform the shift of the gauge
field $\beta\mathit\Phi\mapsto\overline{\mathit\Phi}(\beta,B)$
(\ref{Pbar}). This step preserves consistency of interactions and
inserts $N-1$ coupling parameters $B_1,\ldots,B_{N-1}$\, into the
field equations.

The non-linear theory (\ref{e-int-general}) is stable if the bounded
representative of free energy-momentum tensor series (\ref{Theta1})
is still conserved at the interacting level. The latter requirement
can be interpreted as the selection rule for admissible couplings.
The derived theories, whose characteristic polynomial
(\ref{Charpol}) has non-degenerate real roots, are usually stable at
free level\footnote{Several known examples confirm this observation
\cite{KKL},\cite{KKL18},\cite{AKL18}, though it is not a theorem at
the moment.}. This gives a good chance for existence of stable
interactions between such models. The theories with multiple real
or/and complex roots of the characteristic polynomial usually have
the degrees of freedom whose dynamics is unbounded by any conserved
quantity. These degrees of freedom are the true ghosts unless they
are suppressed by constraints or gauged out. The stability of
interactions in the models having no bounded conserved quantity at
free level is a subtle issue, which we do not consider in the
present paper.

The general representative in the class of interacting theories
(\ref{e-int-general}) is non-Lagrangian as the equations of motion
are not given by the variational derivatives of any functional. The
field equations are Lagrangian if the values of coupling parameters
meet the condition
\begin{equation}\label{Lag-vertex}
    \beta_0=\beta\,,
\end{equation}
and all other constants vanish. The action functional reads
\begin{equation}\label{}
    S=\frac{1}{2}\Big(\langle \phi,M(\beta\mathit\Phi)\phi\rangle+\langle\mathit{\Phi},\mathcal{M}\mathit{\Phi}\rangle\Big)\,.
\end{equation}
The corresponding conserved quantity is the canonical
energy-momentum tensor of the model. In the class of higher
derivative derived equations, the variational interaction vertex
usually suffers from the Ostrogradski instability. In contrast to
Lagrangian interactions, the non-Lagrangian couplings can preserve
bounded representative in the series of free conserved quantities
(\ref{Theta1}). These couplings define the stable non-linear theory.

Let us summarize the results of this section. Relations
(\ref{e-int-general})-(\ref{Pbar}) provide a general receipt of
construction of series of gauge invariant interaction vertices
between two derived theories,  if the gauge invariant  coupling is
known for primary models. The structure of interaction vertex
resembles couplings between electromagnetic field and charged
matter: the wave operator of originally non-gauge field is replaced
by its gauge invariant extension, while the current-like term is
added to equations of motion of gauge field. If the characteristic
polynomials of derived theories have the orders $n$ and $N$, there
are $n+N$ coupling constants. Any specific choice of coupling
constants selects certain representative in the series of free
conserved energy-momentum such that still conserves at the
interacting level. If the free theory has bounded conserved
quantities, it admits the stable interactions, though the stable
couplings are non-Lagrangian.

\section{Couplings between extended Chern-Simons and higher derivative scalar}
Following the general pattern of the previous section,  let us
consider the interactions between the gauge vector field and the
charged scalar. In this case, we have two subsets of fields on $3d$
Minkowski space,
\begin{equation}
\phi^a=(\text{Re}\,\phi)(x)+i(\text{Im}\,\phi)(x)\,, \qquad
\mathit{\Phi}^i=\mathit{\Phi}^\mu(x)\,.
\end{equation}
The primary operators read
\begin{equation}\label{primopSChS}
W=m^{-2}\partial_\mu\partial^\mu\,, \qquad
\mathcal{W}_{\mu\nu}=m^{-1}\varepsilon_{\mu\rho\nu}\partial^\rho\,.
\end{equation}
To make contact with the general scheme of previous section, we
identify the d'Alembertian as the primary operator for the scalar.
For the vector field, the Chern-Simons operator serves as the
primary operator. The constant $m>0$ has the dimension of mass. Here
we use it to make the primary operators dimensionless. The primary
equations of motion,
\begin{equation}\label{primSChS}
m^{-2}\partial_\mu\partial^\mu\phi=0\,, \qquad
m^{-1}\varepsilon_{\mu\nu\rho}\,\partial^\nu\mathit{\Phi}^\rho=0\,,
\end{equation}
have the gauge symmetry (\ref{gt-P}), with the gauge generator for
the vector field being gradient, $R^i{}_\alpha\equiv\partial^\mu
$\,.

The primary models (\ref{primSChS}) admit obvious variational
interaction vertex,
\begin{equation}\label{SSChS}
\displaystyle S[\phi,\mathit{\Phi}]=\frac{1}{2}\int\big[\phi^*D_\mu
D^\mu\phi+m\varepsilon_{\mu\nu\rho}\mathit{\Phi}^\mu\partial^\nu\mathit{\Phi}^\rho\big]\,d^3x\,,
\end{equation}
where $D_\mu$ is the covariant derivative,
\begin{equation}\label{Dmuphi}
D_\mu(\beta)\phi=\big(\partial_\mu-i\beta\mathit{\Phi}_\mu\big)\phi\,.
\end{equation}
For the complex conjugate field $\phi^*$, the covariant derivative
is given by the complex conjugation of (\ref{Dmuphi}). The action
functional (\ref{SSChS}) has form (\ref{S-primary-int}) with the
operator of vertex,
\begin{eqnarray}
\mathit{\Gamma}&\equiv&D_\mu
D^\mu-\partial_\mu\partial^\mu=(\partial_\mu-i\beta\mathit{\Phi}_\mu\big)(\partial^\mu-i\beta\mathit{\Phi}^\mu\big)-\partial_\mu\partial^\mu\,,
\end{eqnarray}
being the second-order polynomial in the coupling constant
$\beta$\,. The Lagrange equations for the action functional
(\ref{SSChS}) read
\begin{eqnarray}
\frac{\delta S}{\delta\mathit{\Phi}^\mu}\equiv
m\varepsilon_{\mu\nu\rho}\partial^\nu\mathit{\Phi}^\rho+i\beta\big(\phi^*D_\mu\phi-\phi
D_\mu\phi^*\big)=0\,,\qquad \frac{\delta S}{\delta
\phi}\equiv\frac{1}{2}D_\mu D^\mu\phi^*=0\,.
\end{eqnarray}
These equations have form (\ref{EL-primary-int}) with
\begin{eqnarray}
J_i\equiv J_\mu(\phi, \beta\mathit{\Phi}) =i(\phi^*D_\mu\phi-\phi
D_\mu\phi^*\big)\,,\qquad W(\beta\mathit{\Phi})=D_\mu D^\mu\,.
\end{eqnarray}
As the gauge field is a vector, the current-like term holds world
index. As we can see, the quantity $J_\mu(\phi, \beta\mathit{\Phi})$
is the charge current of complex scalar field. The infinitesimal
gauge transformation (\ref{gt-pA}) of the action functional
(\ref{SSChS}) reads
\begin{equation}\label{gtSChS}
\delta\phi=i\beta\phi^*\varepsilon\,, \qquad
\delta\mathit{\Phi}^\mu=\partial^\mu\varepsilon\,.
\end{equation}
The canonical energy-momentum tensor (\ref{ThetaSer}) for the action
(\ref{SSChS}) has the form
\begin{equation}\begin{array}{c}
\displaystyle
\boldsymbol{\mathit{\Theta}}_{\mu\nu}=F_{\mu}F_{\nu}-\frac{1}{2}\eta_{\mu\nu}F_{\lambda}F^{\lambda}+D_\mu\phi^*D_\nu\phi+D_\nu\phi^*D_\mu\phi
-\eta_{\mu\nu}D_\lambda\phi^*D^\lambda\phi\,, \\[3mm] F_\mu\equiv
m\varepsilon_{\mu\nu\rho}\partial^\nu\mathit{\Phi}^\rho\,.
\end{array}\end{equation} It is obvious that its $00$-component is bounded.
So, the theory of interacting Chern-Simons and charged scalar is
stable at the level of primary theories.

Now we consider interactions between two derived theories, whose
primary operators are given by (\ref{primopSChS}).
 The free equations (\ref{2-derived}) are chosen in the form
\begin{eqnarray}\label{derSChS}
&&M(\alpha; W)\,\phi\equiv m^2\prod\limits_{p=0}^{n-1}\big(m^{-2}\partial_\mu\partial^\mu+\alpha_p^2\big)\phi=0\,,\notag\\
&&\mathcal{M}_{\mu\nu}(A;\mathcal{W})\,\mathit{\Phi}^\nu\equiv\frac{m^2}{2}\sum\limits_{q=1}^NA_q(\mathcal{W}^q)_{\mu\nu}\mathit{\Phi}^\nu=0\,.
\end{eqnarray}
Here we assume that all the roots of characteristic polynomial
(\ref{Charpol}) for scalar field are different and positive. We
ignore all the other options because they do not result to the
stable theory at the free or interacting level. The extended
covariant derivative $\overline{D}_\mu$ is defined as follows:
\begin{equation}\label{overlineDmu}
\displaystyle
\overline{D}_\mu(\beta,B)\phi=\big(\partial_\mu-i\sum_{q=0}^{N-1}B_q(\mathcal{W}^q)_{\mu\nu}\mathit{\Phi}^\nu\big)\phi\,,
\qquad B_0\equiv \beta\,.
\end{equation}
The gauge invariant extension of the free wave operator $M(\alpha;
W)$ (\ref{derSChS}) reads
\begin{equation}
\overline{M}(\alpha; \overline{W}(\beta\mathit{\Phi}))\,\phi\equiv
m^2\prod\limits_{p=0}^{n-1}\big(m^{-2}\overline{D}_\mu\overline{D}^\mu+\alpha_p^2\big)\phi=0\,.
\end{equation}
The series of current-like terms (\ref{J-1}) can be represented in
the form
\begin{equation}
\displaystyle
(J_p)_\mu=i\big(\phi^{*(p)}\overline{D}_\mu\phi^{(p)}-\phi^{(p)}\overline{D}_\mu\phi^{*(p)}\big)\,,
\, p=0,\ldots,n-1\,,
\end{equation}
where the notation is used:\footnote{The fields $\phi^{(p)}$
describe irreducible components of reducible representation of the
Poincar\'e group that is described by the higher derivative theory
of complex scalar. The vectors $(J_p){}_{\mu}\,,\,
p=0,\ldots,n-1\,,$ are charge currents of irreducible components,
while the coupling constants $\beta_p,\,p=0,\ldots,n-1\,,$ can be
interpreted as the charges of components.}
\begin{equation}\label{phip}
\displaystyle \phi^{(p)}=\prod_{\substack{p'=0\\ p'\neq
p}}^{n-1}\frac{m^{-2}\overline{D}_\mu
\overline{D}^\mu+\alpha_{p'}^2}{\alpha_p^2-\alpha_{p'}^2}\phi\,,\,\,
p=0,\ldots,n-1\,.
\end{equation}
Let us notice that $(J_p)_\mu$ are the linear combinations of
$(J_{p}){}_{i}$ in (\ref{J-1}). We deal with objects $(J_p)_\mu$ for
reasons of convenience.

The $n+N$-parameter series of interaction vertices
(\ref{e-int-general}) for free equations (\ref{derSChS}) read
\begin{widetext}
\begin{eqnarray}\label{ESChS}
&&\displaystyle E\equiv m^2\prod\limits_{p=0}^{n-1}\big(m^{-2}\overline{D}_\mu(\beta,B)\overline{D}^\mu(\beta,B)+\alpha_p^2\big)\phi=0\,,\notag\\
&&\displaystyle
\mathcal{E}_\mu\equiv\frac{m^2}{2}\sum\limits_{q=1}^NA_q(\mathcal{W}^q)_{\mu\nu}\mathit{\Phi}^\nu+i\sum\limits_{p=0}^{n-1}\beta_p\big(\phi^{*(p)}\overline{D}_\mu(\beta,B)\phi^{(p)}-\phi^{(p)}\overline{D}_\mu(\beta,B)\phi^{*(p)}\big)=0\,,
\end{eqnarray}
\end{widetext}
where $\overline{D}_\mu(\beta,B)$ is the covariant derivative
(\ref{overlineDmu}), $\phi^{(p)}$ is given by (\ref{phip}), and
$\beta\,,\,\beta_p\,,\,B_q\,,\,p=0,\ldots,n-1\,,\,q=1,\ldots,N-1\,,$
are the coupling parameters.

Equations (\ref{ESChS}) are obviously gauge invariant with respect
to gauge symmetry (\ref{gtSChS}). The gauge identity between
equations (\ref{ESChS}) is seen in the form
\begin{equation}
\displaystyle
\partial^\mu\mathcal{E}_\mu-i\sum\limits_{p=0}^{n-1}\prod_{\substack{p'=0\\
p'\neq
p}}^{n-1}\frac{\beta_p}{\alpha_p^2-\alpha_{p'}^2}\big(\phi^{*(p)}E-\phi^{(p)}E^*\big)=0\,.
\end{equation}
The series of the energy-momentum tensors (\ref{Theta1}) for
non-linear theory (\ref{ESChS}) has the structure
\begin{equation}\label{ThetaSChS}
\boldsymbol{\mathit{\Theta}}_{\mu\nu}(\phi,\mathit{\Phi})=\sum_{p=0}^{n-1}\beta_p(T_p){}_{\mu\nu}(\phi,\mathit{\Phi})+\sum_{q=0}^{N-1}B_q(T_q)_{\mu\nu}(\mathit{\Phi})\,,
\end{equation}
where the quantities $(T_p){}_{\mu\nu},\,p=0,\ldots,n-1\,,$ and
$(T_q){}_{\mu\nu},\,q=0,\ldots,N-1\,,$ read
\begin{equation}\label{TqChS}\begin{array}{c}\displaystyle
(T_p){}_{\mu\nu}(\phi,\mathit{\Phi})=\overline{D}_\mu\phi^{*(p)}\overline{D}_\nu\phi^{(p)}+
\overline{D}_\nu\phi^{*(p)}\overline{D}_\mu\phi^{(p)}-\eta_{\mu\nu}\overline{D}_\lambda\phi^{*(p)}\overline{D}^\lambda\phi^{(p)}+
m^2\alpha_p^2\eta_{\mu\nu}\phi^{*(p)}\phi^{(p)}\,,\\[3mm]\displaystyle
(T_q){}_{\mu\nu}(\mathit{\Phi})=\sum_{r,s=1}^{N-1}\sum_{q'=1}^{N}C_{r,s}^{q,q'}A_{q'}\big(\mathit{\Phi}^{(r)}{}_{\mu}\mathit{\Phi}^{(s)}{}_{\nu}+\mathit{\Phi}^{(r)}{}_{\nu}\mathit{\Phi}^{(s)}{}_{\mu}-\eta_{\mu\nu}\mathit{\Phi}^{(r)}{}_\lambda\mathit{\Phi}^{(s)\lambda}\big)\,,
\\[3mm]\displaystyle
\mathit{\Phi}^{(q)}{}_{\mu}\equiv(\mathcal{W}^q){}_{\mu\nu}\mathit{\Phi}^{\nu}\,.
\end{array}\end{equation} The constants $C_{r,s}^{q,q'}$ are defined as
follows:
\begin{equation}\label{Crsqq}
C_{r,s}^{q,q'}= \left\{
\begin{array}{rll}
-1, &  q>q'\,, &q>r,s\,,\\
&  &q+q'=r+s+1\,;\\
1, & q<q'\,, &q\leq r,s\,,\\
& &q+q'=r+s+1\,;\\
0, & \text{otherwise}\,.&\\
\end{array}
\right.
\end{equation}
As is seen from (\ref{ThetaSChS}), the  quantity
$\boldsymbol{\mathit\Theta}^{\mu\nu}$ is the deformation of selected
representative of energy-momentum tensor series of free complex
scalar and Chern-Simons,
\begin{eqnarray}\label{ThetaSChS-free}
\boldsymbol{\mathit{\Theta}}^{\text{free}}_{\mu\nu}(\phi,\mathit{\Phi})&=&
\sum_{q=0}^{N-1}B_q(T_q)_{\mu\nu}(\mathit{\Phi})+
\sum_{p=0}^{n-1}\beta_p(\partial_\mu\phi^{*(p)}\partial_\nu\phi^{(p)}\notag\\
&+&\partial_\nu\phi^{*(p)}\partial_\mu\phi^{(p)}
-\eta_{\mu\nu}\partial_\lambda\phi^{*(p)}\partial^\lambda\phi^{(p)}+m^2\alpha_p^2\eta_{\mu\nu}\phi^{*(p)}\phi^{(p)})\,.
\end{eqnarray}
In this formula all the parameters of the conserved quantity are
fixed by the interaction, so a single representative of free series
of conserved quantities survives at the interacting level.

Let us discuss the stability of the theory (\ref{ESChS}) at the
interacting level. The $00$-component of (\ref{ThetaSChS}) can be
written in the form
\begin{eqnarray}\label{Theta00SChS}
\displaystyle \boldsymbol{\mathit{\Theta}}_{00}(\phi,\mathit{\Phi})&=&\frac{m^2}{2}\sum\limits_{r,s=1}^{N-1}C_{r,s}(A,B)\mathit{\Phi}^{(r)}{}_{\mu}\mathit{\Phi}^{(s)}{}_{\mu}\notag\\
&+&\sum\limits_{p=0}^{n-1}\beta_p\big(\overline{D}_\mu\phi^{*(p)}\overline{D}_\mu\phi^{(p)}+m^2\alpha_p^2\phi^{*(p)}\phi^{(p)}\big)\,.
\end{eqnarray}
where the matrix $C_{r,s}(A,B)$ is defined by
\begin{equation}\label{C-matrix}
\displaystyle
C_{r,s}(A,B)=\sum_{q=0}^{N-1}\sum_{q'=1}^{N}C_{r,s}^{q,q'}B_qA_{q'}\,.
\end{equation}
The expression (\ref{Theta00SChS}) is a quadratic form in the
variables $\mathit{\Phi}^{(q)}{}_{\mu}\,, \phi^{(p)}$\,. It is
bounded if $C_{r,s}(A,B)\,,\, r,\,s=1,\ldots,N-1\,,$ is a positive
definite matrix, and $\beta_p>0\,,\,p=0,\ldots,n-1\,,$ are positive
numbers.\footnote{The quantity $C_{r,s}(A,B))$ is the Bezout matrix
of the characteristic polynomial (\ref{Charpol}) of free extended
Chern-Simons theory and characteristic polynomial of derived
symmetry (\ref{Xp-Sym}), (\ref{XSer-Sym}). It is defined by the
relation
\begin{equation*}
\sum_{r,s=1}^{N-1}C_{r,s}(A,B)z^{r}u^{s}=\frac{M(A;z)N(B;u)-M(A;u)N(B;z)}{z-u}\,,
\end{equation*}
where
\begin{equation*}
 M(A;z)=\sum_{q=1}^NA_qz^q\,, \qquad
N(B;z)=\sum_{q=0}^{N-1}B_{q+1}z^{q+1}\,.
\end{equation*}
Accounting this definition, formula (\ref{Crsqq}) takes the form
\begin{equation*}
C_{r,s}^{q,q'}=\frac{\partial^2 C(A,B)}{\partial A_q
\partial B_{q'}}\,.
\end{equation*}
}

We consider the stability condition of interacting theory as the
natural selection rule for admissible couplings. In the sector of
scalar field, it restricts the coupling parameters $\beta_p$. In the
sector of gauge field, the matrix $C_{r,s}(A,B)$ can be positive
definite or indefinite depending on the values of the free model
parameters $A_q$ and coupling constants $B_q$. In the work
\cite{KKL} the following fact is noticed about the free extended
Chern-Simons theory: the series of conserved quantities
(\ref{TqChS}) includes the bounded representative if the
characteristic polynomial of the theory (\ref{derSChS}) has simple
real nonzero roots and a zero root of multiplicity one or two. In
all the other cases, the theory is unstable. For the viewpoint of
representation theory, the stability condition requires from the
free extended Chern-Simons theory to describe the unitary
representation of the Poincar\'e group. The roots of characteristic
equation determine the masses of irreducible representations in this
case. Once the roots of the characteristic equation of  extended
Chern-Simons theory are such that any conserved quantity is
unbounded in the series (\ref{TqChS}), and therefore the stable
interactions are impossible to include, the free theory corresponds
to the non-unitary representation. It seems quite natural that the
non-unitary theory cannot be stable at interacting level.

In general, the interaction (\ref{ESChS}) is non-Lagrangian. The
Lagrangian interaction vertex corresponds to the following values of
coupling parameters:
\begin{equation}
\beta=\beta_0\,, \,\,
B_1=\ldots=B_{N}=\beta_1=\ldots=\beta_{n-1}=0\,.
\end{equation}
The conserved quantity for variational couplings is the canonical
energy-momentum tensor.  As the canonical energy
$(\boldsymbol{\mathit{\Theta}}_\text{can}\,)_{00}$ is unbounded,
the Lagrangian interaction is unstable. The field equations
(\ref{ESChS}) include stable interactions, so the corresponding
coupling is non-Lagrangian. The coupling parameters that lead to
stable interactions are determined by conditions
$\beta_p>0\,,p=0,\ldots,n-1\,,$ and positive definiteness requirement for
the form $C_{r,s}(A,B)$ (\ref{C-matrix}).

\section{Conclusion and discussion of results}
At first, let us summarize the basic assumptions about their
consequences for the class of higher derivative field theories we
study in this paper. We consider inclusion of stable interactions in
certain class of higher derivative field theories which we term as
\textit{derived} models. The derived higher derivative theory
implies that the same field admits free field equation with the wave
operator $W$ without higher derivatives. We call this model a
\textit{primary} theory. The wave operator of the derived theory is
a $n$-th order polynomial in $W$. We call it a characteristic
polynomial of the derived theory. Every symmetry of primary theory
results in the $n$-parametric series of symmetries of the derived
theory. These symmetries are connected to the series of independent
conserved quantities. In particular, the translation invariance of
the primary model results in the $n$-parametric series of conserved
tensors. Under certain assumptions about the roots of characteristic
polynomial, the series can include bounded conserved quantities. The
canonical energy-momentum is included into the series, though it is
always unbounded. Once the free derived theory admits bounded
conserved quantities, it is considered stable.

We consider inclusion of interactions between two different derived
field theories. One of these is supposed to be gauge invariant, and
another one is non-gauge. We assume that the primary theories of
these two models admit consistent and stable interactions, with the
gauge symmetry remaining abelian at interacting level. We also
impose certain technical restrictions on the interaction between
primary theories as explained in the section 3. In particular, these
restrictions exclude the "gauging" scenarios of inclusion of
interactions\footnote{The "gauging" schemes mean that the gauge
symmetry of self-interacting gauge fields can change upon inclusion
of the non-gauge matter fields. Concerning the literature on
gauging, we can mention\cite{gaug1, gaug2} and references therein.
}. In principle, these restrictions could be relaxed without
breaking the general scheme of inclusion of stable interactions
in the derived-type theories, though this would make the
consideration much more cumbersome. Once the primary theories admit
consistent interactions, these can be lifted to the level of derived
theory. The lift is not unique, we get $n+N$-parametric series of
consistent interactions, where $n,N$ are the orders of
characteristic polynomials of the two models. At free level, the two
derived theories admitted $n+N$ parametric series of the conserved
tensors. Upon inclusion of interaction, only one tensor conserves
with any fixed set of the interaction parameters. The canonical
energy momentum is conserved if the interaction is Lagrangian. As
the canonical energy is unbounded, the Lagrangian interaction is
unstable.  Once the free derived theory admits bounded conserved
quantity, its conservation can be always preserved with appropriate
consistent and Poincar\'e invariant interaction, so the stability
persists at interacting level, though with the non-Lagrangian
vertex.

Also notice that even if the interaction is included in such a way
that none of the bounded conserved quantities (from the series
identified in section 2)  of the free theory remains conserving at
interacting level, these non-conserved currents can be still
informative. So they can be still relevant for Lagrangian
interactions. If the dynamics is known of bounded non-conserved
quantity, than we learn about the compact evolving phase-space
surface where the dynamics is confined. This might be useful for
identifying the isles of stability, and evaluating the velocity of
running away solutions. For discussion of these ideas and further
references see the recent work \cite{Nicolas}.

Even though the vertices of stable interactions are always
non-Lagrangian in this class of higher derivative theories, they can
still admit quantization. The matter is that these non-Lagrangian
field equations can admit Lagrange anchor which allows one to
quantize the non-Lagrangian field theory. The examples are provided
in the paper \cite{KLS14} where the some special interactions are
considered between the derived model with the quadratic
characteristic polynomial, and the theory without higher
derivatives. Another side of the existence of the Lagrange anchor is
that the theory should admit Hamiltonian formalism even though the
higher derivative equations are non-Lagrangian. The explicit
examples of Hamiltonian formalism for third-order non-Lagrangian
derived systems can be found in papers \cite{AKL18, KKL18}. The
constrained Hamiltonian formalism for $n$-th order derived theory
with the stable interactions constructed in Section 4 is developed
in the recent work \cite{AKL18-RFJ}.

\begin{acknowledgments} We thank
P.~Kazinski, A.~Kamenshchik, A.~Sharapov, and A.~Vikman for
discussions. The work of VAA is supported by Tomsk State University
Competitiveness Improvement Program. The work of DSC is supported by
the Russian Science Foundation grant 18-72-10123 in association with
the Lebedev Physical Institute of RAS. The work of SLL is supported
by the project 3.5204.2017/6.7 of Russian Ministry of Science and
Education.
\end{acknowledgments}


\section{Introduction}\label{sec:introduction} The quantum complexity of discrete problems is well studied. A speed-up of quantum algorithms over deterministic and randomized algorithms is shown for many
problems, starting from the factorization algorithm of Shor \cite{Shor}, followed by
database search algorithm of Grover \cite{Grover}. Other discrete problems, such
as discrete summation, computation of the mean, median and
$k$th-smallest element were also studied, see e.g.
\cite{Brassard},\cite{Brassard2}\cite{Durr},\cite{Grover2},\cite{Nayak}.
 \\
There is also a progress in studying the quantum complexity of continuous analogues of these problems. The first paper dealing with
the quantum complexity of a continuous problem was the work of Novak
\cite{Novak1}, where the integration of a function from the
H\"older class is considered.  Integration in other
function spaces (Lebesgue and Sobolev classes) space was also
investigated and a quadratic speed-up over the randomized setting
was shown (see \cite{Heinrich},\cite{Heinrich3}). The problem
of function approximation on a quantum computer was studied by
Heinrich \cite{Heinrich1},\cite{Heinrich2}. The problem of maximization of function form the H\"older class was investigated in~\cite{Gocwin} and the problem of finding the root of the function was studied in~\cite{GocwinKacewicz}. Also path
integration \cite{Wozniakowski} and differential equations
\cite{Kacewicz2005} on a quantum computer were
investigated, and a speed-up was established.
\\
In this paper we consider the problem of approximating the median of continuous distribution. This problem is a continuous equivalent of the problem of finding the median of the discrete sequence which is investigated in~\cite{Nayak}. We assume that the distribution is given by density function which belongs to a
classical H\"older class. We present almost matching upper and lower complexity bounds in the
quantum setting. We show that quantum computations yield a
speed-up compared to deterministic and randomized
algorithms over entire range of class parameters.
\\
We extend also the problem to the problem of computing of the vector of quantiles of the continuous distribution and show matching upper and lower complexity bounds for this problem.
\\
In Section \ref{sec:formulation} the problem of computing the mean is formulated and  necessary definitions are presented. The complexity bounds for this problem are presented in Section~\ref{sec:median}. Section~\ref{sec:quantiles} contains the results for the problem of computing of the vector of quantiles.
\section{Problem formulation and basic definitions}\label{sec:formulation}
Let $f:[0,1]\rightarrow {\bf R}$ be a real-valued, nonnegative function which integrates to 1.  We are interested in approximation the median of the probabilistic distribution with probabilistic density function $f$  with precision $\varepsilon>0$
in the sense of the absolute or the residual error criterion. Let $\xi:\ \int\limits_0^{\xi}f(t)dt=1/2$ be the median. We wish to compute a point $\hat{\xi}\in [0,1]$ for which $|\hat{\xi}-\xi|\leq \varepsilon$  (the absolute criterion), or
$\left|\int\limits_{0}^{\hat{\xi}}f(t)dt-\frac{1}{2}\right|\leq\varepsilon$ (the residual criterion). We consider a H\"older class of functions $f$ given by
\begin{eqnarray*}
F^{r,\rho}&=&\left\{f:[0,1] \rightarrow {\bf R} \,\, |\,\, f\in C^r([0,1]), |f^{(i)}(x)|\leq D,\,\, i=0,1,
\ldots r,\right.\\  & & \left.|f^{(r)}(x)-f^{(r)}(y)|\leq H |x-y|^{\rho} \;\mbox{ for } x,y\in[0,1] \right\},
\end{eqnarray*}
where $r\geq 0$, $0<\rho\leq 1$ and  $D,H$ are positive constants.
Let us define the class of density function from the H\"older class $F^{r,\rho}$ by $F^{r,\rho}(1)$, that is
\[
F^{r,\rho}(1)=\left\{f\inF^{r,\rho}\,\,|\,\, f(x)\geq 0\ \forall x\in[0,1],\ \int\limits_0^{1}f(t)dt=1\right\}.
\]
For the absolute criterion we will also need a technical assumption that the function $f$ is separated form zero. Let
\[
\tilde{F}^{r,\rho}(1)=\left\{f\inF^{r,\rho}(1)\,\,|\,\, f(x)\geq \gamma\ \forall x\in[0,1]\right\},
\]
where $\gamma$ is a positive constant.\\
We shall analyze the problem of approximation the median of the distribution with density function $f\inF^{r,\rho}(1)$ or $f\in\tilde{F}^{r,\rho}(1)$ in the quantum setting and for the comparison in the deterministic worst-case and in the randomized  settings.\\
To find the approximation $t^*$ we need some information about the function $f$.
In the deterministic setting, by
information we mean an operator $N$ that assigns to the function $f$ a vector $N(f)$
of length at most $n$, whose components are
the values of $f$ or its derivatives of order at most $r$ at some (adaptively chosen) deterministic points form $[0,1]$.
Let $(\Omega,\Sigma, {\rm P})$ be a probability space.
In the randomized setting, by information we mean a family of operators $N=\{N^\omega\}_{\omega\in \Omega}$.
Here, $N^\omega(f)$ is defined by the values of $f$ or its derivatives computed at randomly chosen points
from $[0,1]$. We assume that the successive evaluation points and the derivative numbers are chosen adaptively as functions of information values computed so far.
The decision whether to compute a successive piece of information is also taken based on the computed  values,
according to some termination criterion.
We assume that the length $n^\omega(f)$
of the information vector  is a measurable function and it satisfies
${\rm E}(n^{\omega}(f))\leq n$ for all $f\in F_d^{r,\rho}$.\\
In the quantum setting, information is gathered by applying $n$ times a (standard) quantum query operator
 about the function
$f$. For a detailed description of the quantum model of computation, quantum algorithms, and the quantum query operator, the reader is referred to \cite{Heinrich, Nielsen}. Roughly speaking, a quantum query is a unitary operator
defined by $f$ whose application plays a role of evaluating the value of a function or its derivative
in the standard deterministic setting. \\
The approximate solution $\hat{\xi}$ is obtained in the deterministic setting by an algorithm $\v$ that computes  $\hat{\xi}$ based on the information values $N(f)$, $\hat{\xi}=\v (N(f))$.
In the randomized setting, the algorithm is a family $\v=\{\v^{\omega}\}_{\omega\in \Omega}$,
and $\hat{\xi}=\hat{\xi}^{ (\omega)}=\v^{\omega}(N^{\omega}(f))$.
In the quantum setting,
the algorithm is defined by an application of a number of sequences of unitary operations, $n$ of them being quantum queries,
each of the sequence followed by a quantum measurement. Further possible operations performed on a classical computer
lead to $\hat{\xi}=\hat{\xi}^{(\omega)}$, see \cite{Heinrich,Nielsen}.
In the randomized and quantum settings, the point $\hat{\xi}=\hat{\xi}^{(\omega)}$
is a random number.\\
We shall consider the absolute and residual error criteria.
The superscript '$\omega$' in the definitions below in the deterministic setting is obviously irrelevant,
 and it will then be omitted.\\
We define the local absolute error of an algorithm $\v$ by
\[
{\rm e}_{\rm abs}^{\omega}(f,\v^\omega)=\left|\hat{\xi}-\xi\right|.
\]
The residual error is given by
\[
{\rm e}_{{\rm res}}^{\omega}(f,\v^\omega)=\left|\int\limits_0^{\hat{\xi}}f(t)dt-1/2\right|.
\]
 In the randomized and quantum settings, we
assume that ${\rm e}_{\rm abs}^{\omega}(f,\v^\omega)$ and ${\rm e}_{{\rm res}}^{\omega}(f,\v^\omega)$ are random variables.\\
Let the subscript '${\rm crit}$' stand for '${\rm abs}$' or '${\rm res}$'.
For $F$ being $F^{r,\rho}(1)$ or $\tilde{F}^{r,\rho}(1)$ the global error in the class $F$ in the deterministic, randomized and quantum settings
 is defined respectively by:\\
-- in the deterministic setting
\[
{\rm e}^{{\rm det}}_{\rm crit}(F,\v)=\sup_{f\in F}{\rm e}_{\rm crit}(f,\v),
\]
-- in the randomized setting
\begin{equation}\label{e_rand}
{\rm e}^{{\rm rand}}_{\rm crit}(F,\v)=\sup_{f\in F}\left(\int\limits_{\Omega}({\rm e}_{\rm crit}^{\omega}
(f,\v^\omega))^2\,d{\rm P}(\omega)\right)^{1/2},
\end{equation}
-- in the quantum setting
\[
{\rm e}^{{\rm quant}}_{\rm crit}(F,\v)=\sup_{f\in F}\; \inf\{\alpha:
{\rm P}(\,{\rm e}_{\rm crit}^{\omega}(f,\v^\omega)>\alpha \,)\leq 1/4\}.
\]
Hence, in all cases the worst case with respect to functions $f$ from $F$ is of interest.
In the randomized setting the error is defined as usual by the expected value, while
in the quantum setting by the demand to achieve the success probability of $\v$ at least 3/4. \\
The cost of an algorithm $\v$ in the deterministic and randomized settings is defined as an (expected)
 number of
computed information values, that is a number of evaluations of function $f$ or its derivatives
which an algorithm is based on. In the quantum setting, the (query) cost is meant as the number of
applications of the quantum query operator. We denote the cost of an algorithm $\v$
in the class $F$ by ${\rm cost}(F,\v)$ with suitable superscripts: '${\rm det}$', '${\rm rand}$' or '${\rm quant}$'
in the proper setting.
\newline
For $\varepsilon>0$, the $\varepsilon$-complexity in the class $F$ is defined
as a minimal cost of an algorithm that approximate the median of a distribution  with a density function from the class $F$ with the precision
at most $\varepsilon$,
\[
 {\rm comp}_{\varepsilon,{\rm crit}}^{\#}(F)=\inf_{\v}\{\, {\rm cost}^{\#}(F,\v):\;\; {\rm e}^{\#}_{\rm crit}
(F,\v)\leq \varepsilon \, \},
\]
where $ \#\in\{{\rm det},{\rm rand},{\rm quant}\}$. In the quantum setting, the complexity defined above is the
quantum query
complexity (the number of qubits can be arbitrary).
It is also possible to consider, for instance, the qubit complexity, defined as the minimal
number of qubits sufficient to achieve the error at most $\varepsilon$ (here the number of queries
can be arbitrary), or other mixed type complexities.\\
Our aim is to show possibly tight upper and lower bounds on ${\rm comp}_{\varepsilon}^{\#}(F)$.
Firstly, we recall in the next section a result on the problem of summation, integration and  solving initial value problems.

\section{Useful results}\label{sec:known}
We will recall in this section the useful known results, that will help us to prove the complexity bounds on the problem of approximating the mean and quantiles of continuous distribution.\\
We start with the problem of computing the mean of a sequence. We will use these results to prove the lower bounds. Let $g:\{0,1,\ldots,N-1\}\rightarrow [0,1]$ be a discrete function. We are interested in approximation of the mean of the sequence $(g(0),g(1),\ldots,g(N-1))$, that is the number $\displaystyle \frac{1}{N}\sum_{i=0}^{N-1}g(i)$ with precision $\varepsilon>0$. Consider the deterministic, randomized and quantum settings. The definition of the error,  the cost and the complexity are similar as in Section~\ref{sec:formulation}. Denote the $\varepsilon$-complexity of this problem by ${\rm comp}_{\varepsilon}(N)$ with a suitable superscript indicating the setting.  From the known complexity bounds for this problem (see \cite{Novak,Mathe} for the deterministic and the randomized setting, \cite{Grover2,Nayak} for the quantum setting) we have that
\begin{itemize}
\item in the deterministic setting   \begin{equation}\label{lower_detst}{\rm comp}_{\varepsilon}^{{\rm det}}(N)=\Theta (N(1-2\varepsilon));\end{equation}
\item in the randomized setting \begin{equation}\label{lower_rand}{\rm comp}_{\varepsilon}^{{\rm rand}}(N)=\Theta(\min\{N,\varepsilon^{-2}\});\end{equation}
\item in the quantum setting \begin{equation}\label{lower_quant}{\rm comp}_{\varepsilon}^{{\rm quant}}(N)=\Theta(\min\{N,\varepsilon^{-1}\}).\end{equation}
\end{itemize}

We will need also the upper bounds on the complexity of approximating the integral, that we will use to show the upper bounds. Let $g:[0,1]\rightarrow {\bf R}$ be a function from H\"older class $F^{r,\rho}$. Our aim is to approximate $\int_0^1g(x)dx$ with precision $\varepsilon$. Denote the complexity of this problem by ${\rm comp}_{\varepsilon}(Int,F^{r,\rho})$ (with suitable superscript). It is known (see e.g.~\cite{Novak} for the deterministic and randomized settings and \cite{Novak1} for the quantum setting), that the complexity of this problem is bounded by
\begin{itemize}
\item in the deterministic setting   \begin{equation}\label{int_det}{\rm comp}_{\varepsilon}^{{\rm det}}(Int,F^{r,\rho})=O(\varepsilon^{-1/{(r+\rho)}});\end{equation}
\item in the randomized setting \begin{equation}\label{int_rand}{\rm comp}_{\varepsilon}^{{\rm rand}}(Int,F^{r,\rho})=O(\varepsilon^{-1/{(r+\rho+1/2)}});\end{equation}
\item in the quantum setting \begin{equation}\label{int_quant}{\rm comp}_{\varepsilon}^{{\rm quant}}(Int,F^{r,\rho})=O(\varepsilon^{-1/{(r+\rho+1)}}).\end{equation}
\end{itemize}
In the randomized setting, besides of the average error, we will also need to know the distribution of the error. Let $\phi$ be the optimal randomized algorithm. Then, there exists positive constant $C$ such that for any $g\inF^{r,\rho}$, $\varepsilon>0$ and $\delta\in(0,1/2)$ we have
\[
{\rm P}\left(\left|\int_0^1g(x)dx-\phi(g)\right|>\varepsilon\right)\leq\delta
\]
with cost
\begin{equation}
{\rm cost}^{{\rm rand}}(\phi)\leq C \varepsilon^{-1/{(r+\rho+1/2)}}\log(1/\delta).
\end{equation}

We will also need some results on the initial-value problems in ordinary differential equations.  The problems have the form
\begin{equation}
\label{problem_initial}
\left\{ \begin{array}{l}
z'(x)=f(z(x)),\qquad x\in[a,b]\\
z(a)=\eta
\end{array}\right.,
\end{equation}
where $a<b$, $f:{\bf R}^d\rightarrow {\bf R}^d$, $z:[a,b]\rightarrow {\bf R}^d$, and
$f$ belongs to the H\"older class
{\setlength\arraycolsep{0pt}\begin{eqnarray*}
F_d^{r,\rho}=\{&&f:{\bf R}^d\rightarrow {\bf R}^d|\ f\in C^r({\bf R}^d), |D^{(i)}f(x)|\leq D,\ i=0,1,\ldots r,\\
&&\|D^{(r)}f(x)-D^{(r)}f(y)\|\leq H\|x-y\|^{\rho},\ \ x,y\in{\bf R}^d\},
\end{eqnarray*}}where $r\geq 0$, $\rho\in(0,1]$, $D$ and $H$ are positive constants. Denote the $\varepsilon$-complexity of this problem by ${\rm comp}_{\varepsilon}(IVP,F_d^{r,\rho})$. There are known complexity bounds of this problem (see \cite{Kacewicz87} for the deterministic setting, \cite{Heinrich_Milla,Daun} for the randomized setting and \cite{Kacewicz2005} for the quantum setting):
\begin{itemize}
\item in the deterministic setting   \begin{equation}\label{ivp_det}{\rm comp}_{\varepsilon}^{{\rm det}}(IVP,F^{r,\rho})=O(\varepsilon^{-1/{(r+\rho)}});\end{equation}
\item in the randomized setting \begin{equation}\label{ivp_rand}{\rm comp}_{\varepsilon}^{{\rm rand}}(IVP,F^{r,\rho})=O(\varepsilon^{-1/{(r+\rho+1/2)}});\end{equation}
\item in the quantum setting \begin{equation}\label{ivp_quant}{\rm comp}_{\varepsilon}^{{\rm quant}}(IVP,F^{r,\rho})=O(\varepsilon^{-1/{(r+\rho+1-\gamma)}}),\end{equation}
\end{itemize}
for arbitrary $\gamma\in(0,1)$ (constant in big-O notation in the quantum setting may depend on $\gamma$).

\section{Complexity bounds}\label{sec:median}
\subsection{Upper complexity bounds}\label{sec:upper}
First we will present a theorem that states the upper complexity bounds on the approximation of the median of absolutely continuous distribution in the deterministic, randomized and quantum setting for both: the absolute and the residual criterion.
\begin{thm}
There exist positive constants $C^{{\rm det}}_{{\rm abs}}$, $C^{{\rm rand}}_{{\rm abs}}$ and $C^{{\rm quant}}_{{\rm abs}}$ which depend on $r$, $\rho$, $D$ and $\gamma$, and $C^{{\rm det}}_{{\rm res}}$, $C^{{\rm rand}}_{{\rm res}}$  and $C^{{\rm quant}}_{{\rm res}}$ which depends on $r$, $\rho$ and $D$, such that for any $\varepsilon>0$ the $\varepsilon$-complexity of median approximation problem satisfies
\[
{\rm comp}_{\varepsilon,{\rm abs}}^{{\rm det}}(\tilde{F}^{r,\rho})\leq C^{{\rm det}}_{{\rm abs}} \left(\frac{1}{\varepsilon}\right)^{1/(r+\rho)}\, \log\frac{1}{\varepsilon},
\]
\[
{\rm comp}_{\varepsilon,{\rm abs}}^{{\rm rand}}(\tilde{F}^{r,\rho})\leq C^{{\rm rand}}_{{\rm abs}} \left(\frac{1}{\varepsilon}\right)^{1/(r+\rho+1/2)}\, \log^2\frac{1}{\varepsilon},
\]
\[
{\rm comp}_{\varepsilon,{\rm abs}}^{{\rm quant}}(\tilde{F}^{r,\rho})\leq C^{{\rm quant}}_{{\rm abs}} \left(\frac{1}{\varepsilon}\right)^{1/(r+\rho+1)}\, \log\frac{1}{\varepsilon}\log\log\frac{1}{\varepsilon},
\]
\[
{\rm comp}_{\varepsilon,{\rm res}}^{{\rm det}}(F^{r,\rho})\leq
 C^{{\rm det}}_{{\rm res}} \left(\frac{1}{\varepsilon}\right)^{1/(r+\rho)}\, \log\frac{1}{\varepsilon},
\]
\[
{\rm comp}_{\varepsilon,{\rm res}}^{{\rm rand}}(F^{r,\rho})\leq
 C^{{\rm rand}}_{{\rm res}} \left(\frac{1}{\varepsilon}\right)^{1/(r+\rho+1/2)}\, \log^2\frac{1}{\varepsilon}\log\log\frac{1}{\varepsilon},
\]
\[
{\rm comp}_{\varepsilon,{\rm res}}^{{\rm quant}}(F^{r,\rho})\leq
 C^{{\rm quant}}_{{\rm res}} \left(\frac{1}{\varepsilon}\right)^{1/(r+\rho+1)}\, \log\frac{1}{\varepsilon}\log\log\frac{1}{\varepsilon},
\]
\end{thm}
\begin{pf}
To proof this theorem we construct  algorithms $\v^{{\rm det}}$, $\v^{{\rm rand}}$ and $\v^{{\rm quant}}$. These algorithms use perturbed  bisection method for solving nonlinear equation $G(x)=0$ for $G(x)=\int\limits_0^x f(t)dt$. For the positive parameter $\varepsilon$ let $x_i$, $i=1,2,\ldots$ be the successive bisection points and $G_i$ be the approximation of $G(x_i)$ computed by optimal deterministic, randomized or quantum integration algorithm with precision $\varepsilon$ with probability at least $1-\delta^{{\rm rand}/{\rm quant}}$ in the randomized and quantum setting, respectively, where $\delta^{{\rm rand}}=\varepsilon^2/\lceil\log \varepsilon^{-1}\rceil$ and $\delta^{{\rm quant}}=1/(4\lceil\log\varepsilon^{-1}\rceil)$. That is
\begin{equation}\label{precision}
|G_i-G(x_i)|\leq \varepsilon \textrm{ with probability at least } 1-\delta^{{\rm rand}/{\rm quant}}.
\end{equation}
The successive interval is chosen based on these approximate values $G_i$. We finish the bisection procedure when for some $i_0$, $|G_{i_0}|\leq \varepsilon$ or $i_0=i_{max}:=\lceil\log \varepsilon^{-1}\rceil$. The algorithms return the value $\hat{\xi}:=x_{i_0}$.\\
Notice that in all previous steps $i=1,2,\ldots,i_0-1$ we have $|G_i|>\varepsilon$. So, due to (\ref{precision}) if $G_i\geq 0$ then $G(x_i)=G_i+G(x_i)-G(i)\geq G_i -\varepsilon\geq 0$ (with probability at least $1-\delta^{{\rm rand}/{\rm quant}}$ in the randomized and quantum settings). Similarly, if $G_i<0$, then with the same probability we have $G(x_i)<0$. So, with probability at least $(1-\delta^{{\rm rand}/{\rm quant}})^{i_{max}}$ last bisection interval contains the median. So, if the algorithm finishes when $i_0=i_{max}$, then  with probability at least $(1-\delta^{{\rm rand}/{\rm quant}})^{i_{max}}$ we have $\xi-\hat{\xi}\leq \varepsilon$. With the same probability, if the algorithm finishes when $|G_{i_0}|\leq\varepsilon$, then
\[
|G(x_{i_0})|\leq |G(x_{i_0})-G_{i_0}|+|G_{i_0}|\leq 2\varepsilon.
\]
So, we have
\begin{equation}\label{precision_solution}
|G(\hat{\xi})-G(\xi)|\leq 2\varepsilon \textrm{ or } |\xi-\hat{\xi}|\leq\varepsilon \textrm{ with probability at least } (1-\delta^{{\rm rand}/{\rm quant}})^{i_{max}}.
\end{equation}
We now state the error bounds. Consider first the residual criterion. Due to (\ref{precision_solution}), $|G(\hat{\xi})-G(\xi)|\leq 2\varepsilon$ or
$
|G(\hat{\xi})-G(\xi)|\leq \sup_{x\in[0,1]}|G'(x)|\,|\hat{\xi}-\xi|\leq D|\hat{\xi}-\xi|\leq D\varepsilon$. Thus,
\[
e_{{\rm res}}(f,\v^{{\rm det}/{\rm rand}/{\rm quant}})\leq \max\{2,D\}\varepsilon.
\]
This holds with probability at least $(1-\delta^{{\rm rand}})^{i_{max}}\geq 1-i_{max}\,\delta^{{\rm rand}}\geq 1-\varepsilon^2$ in the randomized setting and with probability at least $(1-\delta^{{\rm quant}})^{i_{max}}\geq 1-i_{max}\,\delta^{{\rm quant}}\geq 3/4$ in the quantum settings.

Hence we have
\begin{equation}\label{e_res_det}
{\rm e}^{{\rm det}}_{{\rm res}}(F^{r,\rho},\v^{{\rm det}})\leq\max\{2,D\}\varepsilon
\end{equation}
and
\begin{equation}\label{e_res_quant}
{\rm e}^{{\rm quant}}_{{\rm res}}(F^{r,\rho},\v^{{\rm quant}})\leq\max\{2,D\}\varepsilon.
\end{equation}
To get the error bound in the randomized setting, note that always $e_{{\rm res}}^{\omega}(f,\v^{{\rm rand}})\leq 1/2$. So we have in the randomized setting
\begin{multline*}
\left({\rm e}^{{\rm rand}}_{{\rm res}}(F^{r,\rho},\v)\right)^2=\sup_{f\inF^{r,\rho}}\int\limits_{\Omega}({\rm e}^{\omega}_{{\rm res}}(f,\v^{{\rm rand}}))^2d{\rm P}(\omega)\\
=\sup_{f\inF^{r,\rho}}\Big(\int\limits_{{\rm e}^{\omega}_{{\rm res}}(f,\v^{{\rm rand}})\leq \max\{2,D\}\varepsilon}({\rm e}^{\omega}_{{\rm res}}(f,\v^{{\rm rand}}))^2d{\rm P}(\omega)+\\
\int\limits_{{\rm e}^{\omega}_{{\rm res}}(f,\v^{{\rm rand}})>\max\{2,D\} \varepsilon}({\rm e}^{\omega}_{{\rm res}}(f,\v^{{\rm rand}}))^2d{\rm P}(\omega)\Big)
\leq \sup_{f\inF^{r,\rho}}\left(\max\{2,D\}\varepsilon^2+1/4\varepsilon^2\right)
\end{multline*}
Thus,
\begin{equation}\label{e_res_rand}
e^{{\rm rand}}_{{\rm res}}(F^{r,\rho},\v)/\leq(\max\{2,D\}+1/2)\,\varepsilon.
\end{equation}
Let us pass to the absolute criterion. Let us remind that we have additional assumption, that $f(x)\geq \gamma$ for all $x\in[0,1]$. Note that $G(\hat{\xi})-G(\xi)=f(c)(\hat{\xi}-\xi)$ for some $c\in[0,1]$. Hence, we have $|\hat{\xi}-\xi|\leq 1/\gamma |G(\hat{\xi})-G(\xi)|$. So, due to (\ref{precision_solution}), $|\hat{\xi}-\xi|\leq \max\{1,2/\gamma\}\varepsilon$ with probability at least $1-i_{max}\,\delta^{{\rm rand}/{\rm quant}}$ in the randomized and quantum settings. Hence, similarly as for the residual criterion, we have
\begin{equation}\label{e_abs_det}
{\rm e}^{{\rm det}}_{{\rm abs}}(\tilde{F}^{r,\rho},\v^{{\rm det}})\leq \max\{1,2/\gamma\}\varepsilon,
\end{equation}
\begin{equation}\label{e_abs_rand}
{\rm e}^{{\rm rand}}_{{\rm abs}}(\tilde{F}^{r,\rho},\v^{{\rm rand}})\leq (\max\{1,2/\gamma\}+1/2)\varepsilon
\end{equation}
and
\begin{equation}\label{e_abs_quant}
{\rm e}^{{\rm quant}}_{{\rm abs}}(\tilde{F}^{r,\rho},\v^{{\rm quant}})\leq \max\{1,2/\gamma\}\varepsilon.
\end{equation}
Let us state the cost bounds. In the deterministic setting, the cost of one bisection step for both criteria is of order
$O\left(\left(1/\varepsilon\right)^{1/(r+\rho)}\right)$.
The maximal number of steps is at most  $i_{max}=\lceil\log\varepsilon^{-1}\rceil$. So, the total cost is
\begin{equation}\label{cost_det}
{\rm cost}^{{\rm det}}(F^{r,\rho},\v^{{\rm det}})\asymp {\rm cost}^{{\rm det}}(\tilde{F}^{r,\rho},\v^{{\rm det}})=O\left(\left(\frac{1}{\varepsilon}\right)^{1/(r+\rho)}\log\varepsilon^{-1}\right).
\end{equation}
In the randomized setting the cost of one bisection step is
\begin{align*}
O\left(\left(\frac{1}{\varepsilon}\right)^{1/(r+\rho+1/2)}\log(1/\delta^{{\rm rand}})\right)&=O\left(\left(\frac{1}{\varepsilon}\right)^{1/(r+\rho+1/2)}\log\frac{\log\varepsilon^{-1}}{\varepsilon^2}\right)\\
&=O\left(\left(\frac{1}{\varepsilon}\right)^{1/(r+\rho+1/2)}\log\varepsilon^{-1}\right).
\end{align*}
So, the total cost in the randomized setting is
\begin{equation}\label{cost_rand}
{\rm cost}^{{\rm rand}}(F^{r,\rho},\v^{{\rm rand}})\asymp {\rm cost}^{{\rm rand}}(\tilde{F}^{r,\rho},\v^{{\rm rand}})=O\left(\left(\frac{1}{\varepsilon}\right)^{1/(r+\rho)}\log^2\varepsilon^{-1}\right).
\end{equation}
In the quantum setting, cost of one bisection step is
\[
O\left(\left(\frac{1}{\varepsilon}\right)^{1/(r+\rho+1)}\log(1/\delta^{{\rm quant}})\right)=O\left(\left(\frac{1}{\varepsilon}\right)^{1/(r+\rho+1)}\log\log\varepsilon^{-1}\right).
\]
So, the total cost is of order
\begin{equation}\label{cost_quant}
{\rm cost}^{{\rm quant}}(F^{r,\rho},\v)\asymp {\rm cost}^{{\rm quant}}(\tilde{F}^{r,\rho},\v)=O\left(\left(\frac{1}{\varepsilon}\right)^{1/(r+\rho+1)}\log\varepsilon^{-1}\log\log\varepsilon^{-1}\right).
\end{equation}
Comparing the error bounds (\ref{e_res_det}), (\ref{e_res_rand}), (\ref{e_res_quant}), (\ref{e_abs_det}), (\ref{e_abs_rand}), (\ref{e_abs_quant})  with the cost bounds (\ref{cost_det}), (\ref{cost_rand}), (\ref{cost_quant}) we get the desired complexity bounds.
\end{pf}

\subsection{Lower complexity bounds}\label{sec:lower}
The following theorem presents the lower bounds on the complexity of problem of approximating the median in the deterministic, randomized and quantum settings.
\begin{thm}
There exist positive constants $c^{{\rm det}}_{{\rm abs}}$, $c^{{\rm rand}}_{{\rm abs}}$, $c^{{\rm quant}}_{{\rm abs}}$, $c^{{\rm det}}_{{\rm res}}$, $c^{{\rm rand}}_{{\rm res}}$, $c^{{\rm quant}}_{{\rm res}}$ and $\varepsilon_0$, such that for any $\varepsilon\in(0,\varepsilon_0)$
\[
{\rm comp}_{\varepsilon,{\rm res}}^{{\rm det}}(F^{r,\rho}(1))\geq c^{{\rm det}}_{{\rm abs}}\left(\frac{1}{\varepsilon}\right)^{{1}/({r+\rho})},
\]
\[
{\rm comp}_{\varepsilon,{\rm res}}^{{\rm rand}}(F^{r,\rho}(1))\geq c^{{\rm rand}}_{{\rm abs}}\left(\frac{1}{\varepsilon}\right)^{{1}/({r+\rho+1/2})},
\]
\[
{\rm comp}_{\varepsilon,{\rm res}}^{{\rm quant}}(F^{r,\rho}(1))\geq c^{{\rm quant}}_{{\rm abs}}\left(\frac{1}{\varepsilon}\right)^{{1}/({r+\rho+1})},
\]
\[
{\rm comp}_{\varepsilon,{\rm abs}}^{{\rm det}}(\tilde{F}^{r,\rho}(1))\geq c^{{\rm det}}_{{\rm abs}}\left(\frac{1}{\varepsilon}\right)^{{1}/({r+\rho})},
\]
\[
{\rm comp}_{\varepsilon,{\rm abs}}^{{\rm rand}}(\tilde{F}^{r,\rho}(1))\geq c^{{\rm rand}}_{{\rm abs}}\left(\frac{1}{\varepsilon}\right)^{{1}/({r+\rho+1/2})},
\]
\[
{\rm comp}_{\varepsilon,{\rm abs}}^{{\rm quant}}(\tilde{F}^{r,\rho}(1))\geq c^{{\rm quant}}_{{\rm abs}}\left(\frac{1}{\varepsilon}\right)^{{1}/({r+\rho+1})},
\]
\end{thm}
\begin{pf}
We will proof the theorem above by reducing the problem of approximating the median of the distribution to the problem of approximating the mean of the discrete sequence.
Let $\varepsilon_1>0$ be a parameter to be specified later on.  The class $F^{r,\rho}$ contains $n=\Theta\left(\varepsilon_1^{-1/(r+\rho)}\right)$ functions $h_i$, $i=1,2,\ldots,n$, with disjoint supports in the interval $[0,1/4]$ and $g_i$, $i=1,2,\ldots,n$ with disjoint supports in the interval $[3/4,1]$ such that \[\int\limits_0^1h_i(x)dx=\int\limits_0^1g_i(x)dx=\varepsilon_1^{1+1/(r+\rho)},\quad
\max_{x\in[0,1]}h_i(x)=\max_{x\in[0,1]}g_i(x)=c\,\varepsilon_1,\quad i=1,2,\ldots n
\]
for some constant $c$
 (see~\cite{Novak},~p.~35).
Let $x_1,x_2,\ldots,x_n$ be a sequence of real numbers in $[0,1]$. Then the function $\displaystyle f_{\varepsilon_1}(x)=1+\sum_{i=1}^nx_ih_i(x)-\sum_{i=1}^nx_ig_i(x)$ belongs to the class $\tilde{F}^{r,\rho}(1)\subsetF^{r,\rho}(1)$. \\
Note that for sufficiently small $\varepsilon_1$ we have $2/3\leq f_{\varepsilon_1}(x)\leq 2$ for all $x\in[0,1]$. Let $\xi$ be the median of $f_{\varepsilon_1}$. Since $1/2=\int_0^{\xi}f_{\varepsilon_1}(x)dx=\xi\,f_{\varepsilon_1}(\eta)$ for some $\eta\in[0,1]$, thus $\xi=1/(2f_{\varepsilon_1}(\eta))$ and $1/4\leq \xi\leq 3/4$. This yields that
\begin{align*}
1/2&=\int\limits_0^{\xi}f_{\varepsilon_1}(x)dx=\int\limits_0^{\xi}\left(1+\sum_{i=1}^nx_ih_i(x)\right)dx=\int\limits_0^{1/4}\left(1+\sum_{i=1}^nx_ih_i(x)\right)dx+\int\limits_{1/4}^{\xi}1dx\\
&=1/4 + \varepsilon_1^{1+1/(r+\rho)}\sum_{i=1}^n x_i+\xi-1/4.
\end{align*}
Thus
\[
\sum_{i=1}^nx_i=\frac{1/2-\xi}{\varepsilon_1^{1+1/(r+\rho)}}.
\]
Suppose that algorithm $\v$ (deterministic, randomized or quantum) computes the median of the distribution with a density function $f$ with error at most $\varepsilon$ and cost $N$, for any function $f\in\tilde{F}^{r,\rho}(1)$, in particular for $f=f_{\varepsilon_1}$. Denote the result of the algorithm $\v$ for function $f_{\varepsilon_1}$ by $\hat\xi$.  Note that
\[
\left|\frac{1}{n}\sum_{i=1}^nx_i-\frac{1/2-\hat\xi}{n\varepsilon_1^{1+1/(r+\rho)}}\right|=\frac{|\hat\xi-\xi|}{n\varepsilon_1^{1+1/(r+\rho)}}=\frac{\left|\int_0^{\hat\xi}f_{\varepsilon_1}(x)dx-1/2\right|}{|f_{\varepsilon_1}(\eta)|n\varepsilon_1^{1+1/(r+\rho)}}
\]
for some $\eta\in[0,1]$. Hence, algorithm $\v$ also computes the mean $\displaystyle\frac{1}{n}\sum_{i=1}^n x_i$ with error at most $\displaystyle\frac{C\varepsilon}{n\varepsilon_1^{1+1/(r+\rho)}}$, where $C=1$ for residual criterion and $C=1/\gamma$ when the absolute criterion is used.

We now use lower complexity bounds for the problem of computing the mean of $n$ real numbers from $[0,1]$ presented in Section~\ref{sec:known}.\\
Consider first the deterministic setting. Recall that $n=\Theta\left(\varepsilon_1^{-1/(r+\rho)}\right)$. Let $G$ be a constant in the lower bound in the "$\Theta$" notation.
Take $\varepsilon_1=4C\varepsilon/G$. From (\ref{lower_detst}) we have that the cost of algorithm $\v$ is bounded by
\begin{eqnarray*}
\displaystyle
N&=&\Omega\left( n\left(1-\frac{2C\varepsilon}{n\varepsilon_1^{1+1/(r+\rho)}}\right)\right) =\Omega\left(n-\frac{2C\varepsilon}{\varepsilon_1^{1+1/(r+\rho)}}\right)\\
&=&\Omega\left(G\varepsilon_1^{-1/(r+\rho)}-\frac{1}{2}G\varepsilon_1\varepsilon_1^{-1-1/(r+\rho)}\right)=\Omega\left(\varepsilon^{-1/(r+\rho)}\right).
\end{eqnarray*}
\\
In the randomized setting, we take $\displaystyle\varepsilon_1=\varepsilon^{(r+\rho)/(r+\rho+1/2)}$. From (\ref{lower_rand}) we have
\[N=\Omega\left(\min\left\{n,\left(n\varepsilon_1^{(r+\rho+1/2)/(r+\rho)}\varepsilon^{-1}/C\right)^2\right\}\right) =\Omega\left(\varepsilon^{-\frac{2}{2(r+\rho)+1}}\right)=\Omega\left(\varepsilon^{-\frac{1}{r+\rho+1/2}}\right).\]
\\
In the quantum setting, we take $\displaystyle\varepsilon_1=\varepsilon^{(r+\rho)/(r+\rho+1)}$. From (\ref{lower_quant}) we have
\[N=\Omega\left(\min\left\{n,n\varepsilon_1^{(r+\rho+1)/(r+\rho)}\varepsilon^{-1}/C\right\}\right)= \Omega\left(\varepsilon^{-1/(r+\rho+1)}\right).\]
The cost of any algorithm for computing the median must be at least $N$. This yields the desired lower bounds on the complexity in the deterministic, randomized and quantum settings.
\end{pf}
Note that, since the lower and the upper complexity bounds match up to the logarithmic factor, the algorithm presented in Section~\ref{sec:upper} is almost optimal.
\section{Computing  of quantiles}\label{sec:quantiles}
\subsection{Problem formulation}
In this section we consider the problem of approximation of the vector of quantiles of absolute continuous distribution. Let $f$ be density function. Suppose that for $k\in{\bf N}$ we are given a vector $\alpha=(\alpha_1,\alpha_2,\ldots,\alpha_k)\in[0,1]^k$. Our aim is to approximate the vector of quantiles $\xi=(\xi_1,\xi_2,\ldots,\xi_k)$, such that $\int_0^{\xi_1}f(x)dx=\alpha_1,\ \int_0^{\xi_2}f(x)dx=\alpha_2,\ldots,\int_0^{\xi_k}f(x)dx=\alpha_k$. We assume here that function $f$ is separated from zero and belongs to the H\:older class $\tilde{F}^{r,\rho}(1)$ defined in Section~\ref{sec:formulation}. We use the absolute error criterion, thus the local error of algorithm $\v$ approximating the vector of quantiles $\xi$ for the density function $f$ is defined by:
\begin{equation}
{\rm e}(f,\v)=\max_{i=1,\ldots,k}|\hat\xi_i-\xi_i|,
\end{equation}
where $\hat\xi=(\xi_1,\xi_2,\ldots,xi_k)$ is the approximation returned by the algorithm $\v$. The global error in the class $\tilde{F}^{r,\rho}$, the cost and the complexity are defined similarly as in Section~\ref{sec:formulation}.
\subsection{Complexity bounds}
The following theorem presents the complexity bounds for the problem of approximation of the vector of the quantiles.

\begin{thm}
There exist positive constants $C^{{\rm det}}$, $C^{{\rm rand}}$, $C^{{\rm quant}}$, $c^{{\rm det}}$, $c^{{\rm rand}}$, $c^{{\rm quant}}$ and $\varepsilon_0$ which depend on $r$, $\rho$, $D$ and $\gamma$, such that for any $k\in {\bf N}$, $\alpha\in[0,1]^k$ and $\varepsilon>\varepsilon_0$ the $\varepsilon$-complexity of quantiles approximation problem satisfies
\[
c^{{\rm det}} \left(\frac{1}{\varepsilon}\right)^{1/(r+\rho)}\leq{\rm comp}_{\varepsilon,{\rm abs}}^{{\rm det}}(\tilde{F}^{r,\rho})\leq C^{{\rm det}} \left(\frac{1}{\varepsilon}\right)^{1/(r+\rho)},
\]
\[
c^{{\rm rand}} \left(\frac{1}{\varepsilon}\right)^{1/(r+\rho+1/2)}\leq{\rm comp}_{\varepsilon,{\rm abs}}^{{\rm rand}}(\tilde{F}^{r,\rho})\leq C^{{\rm rand}} \left(\frac{1}{\varepsilon}\right)^{1/(r+\rho+1/2)},
\]
\[
c^{{\rm quant}} \left(\frac{1}{\varepsilon}\right)^{1/(r+\rho+1)}\leq{\rm comp}_{\varepsilon,{\rm abs}}^{{\rm quant}}(\tilde{F}^{r,\rho})\leq C^{{\rm quant}} \left(\frac{1}{\varepsilon}\right)^{1/(r+\rho+1-\delta)},
\]
\end{thm}
\begin{pf}
Let $f\in\tilde{F}^{r,\rho}(1)$. Define the function $F:[0,1]\rightarrow[0,1]$ by $F(x)=\int_{0}^xf(x)dx$. Let $G=F^{-1}$ be the inverse of $F$. Then, $\xi=(\xi_1,\xi_2,\ldots,\xi_k)=(G(\alpha_1),G(\alpha_2),\ldots,G(\alpha_k))$. Note that $G$ is the solution of the initial-value problem
\[
G'(x)=\frac{1}{f(x)},\ \ x\in[0,1],\quad G(0)=0.
\]
Since $f\in\tilde{F}^{r,\rho}$, the function $1/f(x)$ is also H\"older continuous with the same parameters $r$ and $\rho$, but with different constant $\tilde D$ dependent on $D$, $r$ and $\gamma$. Let $l$ be the approximation of $G$ returned by the optimal algorithm solving initial-value problems in the H\"older class (deterministic, randomized or quantum). Then, the approximation of $\xi$ is defined by
\[
\tilde\xi=(\tilde\xi_1,\tilde\xi_2,\ldots,\tilde\xi_k)=(l(\alpha_1),l(\alpha_2),\ldots,l(\alpha_k)).
\]
The upper bounds on the complexity of the problem of approximating of the vector of quantiles follow directly from the upper bounds on the problem of solving the initial-value problems presented in Section~\ref{sec:known}.

The lower bounds follow from the lower bounds on the complexity of approximating of the median presented in Section~\ref{sec:lower}.
\end{pf}
\section{Remarks}
Note that the upper bounds for the problem of approximating the vector of quantiles improve slightly the bounds on the problem of approximating the median in the deterministic and randomized setting for the absolute error criterion. This bounds match the lower complexity bounds. All the other bounds are almost sharp.