text:
\section*{Introduction}

A well known problem is to determine which smooth  domains into an almost complex manifold $(M,J)$ are locally complete hyperbolic in the sense of Kobayashi. It seems natural to make some curvature assumptions on such domains as any $J$-pseudoconcavity boundary point is at finite Kobayashi distance to the interior of the domain \cite{iv-ro}.
I.Graham \cite{gra} gave asymptotic estimates of the Kobayashi pseudometric for strictly pseudoconvex bounded domains into $(\C^n,J_{st})$ and proved the complete hyperbolicity of those domains. In case $(M,J)$ is  any almost complex manifold, similar results were provided  by S.Ivashkovich-J.-P.Rosay \cite{iv-ro} and H.Gaussier-A.Sukhov \cite{ga-su}.

The situation is far from being so clear when the domains are weakly pseudoconvex as the geometry of their boundary is much more complicated; and the question of whether a smooth weakly pseudoconvex domain is locally complete hyperbolic or not is still open, even in $(\C^2,J_{st})$. However, D.Catlin \cite{ca} obtained local estimates similar to those obtained in \cite{gra} on smooth pseudoconvex domains of finite type in $(\C^2,J_{st})$, implying  their local complete hyperbolicity (see also \cite{be2}).

In the present paper we study this question for smooth $J$-pseudoconvex domains of finite type into a four dimensional almost complex manifold. In \cite{ber}, we described locally finite type domains $D=\{\rho<0\}$ where
$\rho$ is a smooth defining function for $D$ and $J$-plurisubharmonic (see Proposition $2.6$ in \cite{ber}).
More precisely, if $D=\{\rho<0\}\subset \C^2$ and if the origin $0 \in \partial D$ is of finite type $2m$, we proved that there is a change of coordinates in a neighborhood of the origin such that the structure $J$ and the
function $\rho$ can be locally written:
\begin{equation}\label{eqstrintro}
J=\left(\begin{array}{ccccc}
a_1& b_1 & 0 & 0 \\
c_1 & -a_1 & 0 & 0\\
0  & 0 & a_2 & b_2\\
0 & 0 & c_2 & -a_2\\
\end{array}\right)=J_{st}+O\left(|z_2|\right),
\end{equation}
and,
\begin{equation}\label{eqintro}
\rho=\Re e z_2+H_{2m}\left(z_1,\overline{z_1}\right)+\widetilde{H}(z_1,z_2)+
O\left(|z_1|^{2m+1}+|z_2||z_1|^m+|z_2|^2\right)
\end{equation}
where  $H_{2m}$ is a homogeneous polynomial of degree $2m$, subharmonic which  is not harmonic and
$\widetilde{H}(z_1,z_2)=\Re e \displaystyle \sum_{k=1}^{m-1} \rho_{k} z_1^k z_2.$

In this paper, we suppose that the harmonic term $\widetilde{H}(z_1,z_2)$ in the above expression (\ref{eqintro})
is identically zero and we prove:
\vskip 0,3cm

\begin{theo}\label{theo1}
Let $J$ be a smooth almost complex structure defined on $\R^4$.
Let $D=\{\rho<0\}$ be a  domain of finite type in $(\R^4,J)$,
where $\rho$ is a smooth defining function of $D$, $J$-plurisubharmonic
in a neighborhood of $\overline{D}$. We suppose that $J$ and $\rho$ satisfy respectively (\ref{eqstrintro}) and
(\ref{eqintro}). Moreover we assume that $\widetilde{H}(z_1,z_2)$ in (\ref{eqintro}) is identically zero. Then there exists a neighborhood $U$ of the origin for which $0$  is at infinite distance from points in $D\cap U$.
\end{theo}

\vskip 0,3cm

The proof of this theorem is inspired by \cite{iv-ro} and is based on the construction of good $J$-plurisubharmonic functions whose use is significant in almost complex manifolds.

\vskip 0,7cm

\begin{flushleft}
{\it Acknowledgments.} I would like to thank J. Bland and J.-P.Rosay for helpful discussions. I am particularly
indebted to J.-P.Rosay who pointed out an erroneous argument in the previous version of my paper.
\end{flushleft}

\section{Preliminaries}
We denote by $\Delta$ the unit disc of $\C$ and by $\Delta_{r}$
the disc of $\C$ centered at the origin of radius $r>0$.
\subsection{Almost complex manifolds and pseudoholomorphic discs}

An {\it almost complex structure} $J$ on a real smooth manifold $M$ is a $\left(1,1\right)$ tensor field
 which  satisfies $J^{2}=-Id$. We suppose that $J$ is smooth.
The pair $\left(M,J\right)$ is called an {\it almost complex manifold}. We denote by $J_{st}$
the standard integrable structure on $\C^{n}$ for every $n$.
A differentiable map $f:\left(M',J'\right) \longrightarrow \left(M,J\right)$ beetwen two almost complex manifolds is said to be
 {\it $\left(J',J\right)$-holomorphic}  if $J\left(f\left(p\right)\right)\circ d_{p}f=d_{p}f\circ J'\left(p\right),$
for every $p \in M'$. In case  $M'=\Delta \subset \C$, such a map is called a {\it pseudoholomorphic disc}.
If $f:\left(M,J\right)\longrightarrow M'$ is a diffeomorphism, we define an almost complex structure, $f_{*}J$,
on $M'$ as the
 direct image of $J$ by $f$ :
$$f_{*}J\left(q\right):=d_{f^{-1}\left(q\right)}f\circ J\left(f^{-1}\left(q\right)\right)\circ d_{q}f^{-1},$$
 for every  $q \in M'$.

\vspace{0,5cm}

The following lemma (see \cite{ga-su}) states that locally any almost
 complex manifold can be seen as the unit ball of
$\C^n$ endowed with a small smooth perturbation of the standard
 integrable structure $J_{st}$.
\begin{lem}\label{ilemloc}
Let $\left(M,J\right)$ be an almost complex manifold, with $J$ of class
 $\mathcal{C}^{k}$, $k\geq 0$.
Then for every point $p \in M$ and every $\lambda_0 > 0$ there exist a neighborhood $U$ of $p$ and
 a coordinate diffeomorphism $z: U \rightarrow \B$ centered a $p$ (ie $z(p) =0$) such that  the
direct image of $J$ satisfies $z_{*}J\left(0\right) = J_{st}$ and
$||z_*\left(J\right) - J_{st}||_{\mathcal{C}^k\left(\bar {\B}\right)}
 \leq \lambda_0$.
\end{lem}
This is simply done by considering a local chart $z: U \rightarrow \B$ centered a $p$ (ie $z(p) =0$), composing it
with a linear diffeomorphism to insure $z_{*}J\left(0\right) = J_{st}$ and  dilating coordinates.

\vspace{0,5cm}

So let  $J$ be an almost complex structure defined in  a neighborhood $U$ of the origin in $\R^{2n}$, and such
that $J$ is sufficiently closed to the standard structure in uniform norm on the closure $\overline{U}$ of $U$.
The $J$-holomorphy equation for a pseudoholomorphic disc
$u : \Delta \rightarrow U \subseteq \R^{2n}$ is given by
\begin{equation}\label{eqholo}
\frac{\partial u}{\partial y}-J\left(u\right)\frac{\partial u}{\partial x}=0.
\end{equation}
According to \cite{ni-wo}, for every $p \in M$, there is a neighborhood $V$ of zero in $T_{p}M$, such that for every $v \in V$,
there is a $J$-holomorphic disc $u$ satisfying $u\left(0\right)=p$ and $d_{0}u\left(\partial/\partial x\right)=v$.

\subsection{Levi geometry}

Let $\rho$ be a $\mathcal{C}^2$ real valued function on
a smooth almost complex manifold  $\left(M,J\right).$
We denote by $d^c_J\rho$ the differential
form defined by
$$d^c_J\rho\left(v\right):=-d\rho\left(Jv\right),$$
where  $v$ is a section of $TM$.
The {\it Levi form} of $\rho$ at a point $p\in M$ and a vector
$v \in T_pM$ is defined by
\begin{equation*}
\mathcal{L}_J\rho\left(p,v\right):=d\left(d^c_J\rho\right)(p)\left(v,J(p)v\right)=dd^c_J\rho(p)\left(v,J(p)v\right).
\end{equation*}

\vspace{0,5cm}

The next proposition is useful in order to compute the Levi form (see \cite{iv-ro}).
\begin{prop}\label{proplevi}\mbox{ }
Let $p\in M$ and $v\in T_pM$. Then
$$\mathcal{L}_J\rho\left(p,v\right)=\Delta \left(\rho \circ u\right)
 \left(0\right),$$
where $u : \Delta \rightarrow \left(M,J\right)$ is any $J$-holomorphic
 disc satisfying
$u\left(0\right)=p$ and $d_0u\left(\partial/\partial_x\right)=v$.
\end{prop}

If $\mathcal{L}_{J}\rho(p,v)\geq 0$ for every $p \in M$ and every $v \in T_pM$, we say that $\rho$ is
 {\it $J$-plurisubharmonic}. It is a well know fact that $\rho$ is $J$-plurisubharmonic if
and only if for every $J$-holomorphic disc $u: \Delta \rightarrow M$,
$\rho \circ u$ is subharmonic (see \cite{iv-ro}).

\subsection{Pseudoconvex domains of finite type}

In this section, we recall some facts about pseudoconvex domains of finite type in four dimensional almost complex manifolds  (See \cite{ber} for more detailed facts).

Let $D=\{\rho<0\}$ is a smooth domain in $\R^4$. Assume that $\rho$ is $J$-plurisubharmonic on a neighborhood of $\overline{D}$
where the structure $J$ is defined on a fixed neighborhood $U$ of $\overline{D}$.
Moreover we suppose that the origin  is a boundary point of $D$.
\begin{defi}\label{defco}
Let $u~: \left(\Delta,0\right)\rightarrow \left(\R^{4},0,J\right) $ be a $J$-holomorphic disc satisfying $u\left(0\right)=0$.
The order of contact $\delta_0\left(\partial D,u\right)$ with $\partial D$ at the origin is the degree of the first term
in the Taylor expansion of $\rho \circ u$. We denote by $\delta\left(u\right)$ the multiplicity of $u$ at the origin.
\end{defi}

We now define the (D'Angelo)  type and the regular type of the real hypersurface $\partial D$ at the origin.

\begin{defi}\label{deftyp}\mbox{ }

 \begin{enumerate}
\item The (D'Angelo) type  of $\partial D$ at the origin is defined by:
$$\Delta^1\left(\partial D,0\right):=\sup \left\{\frac{\delta_0\left(\partial D,u\right)}{\delta\left(u\right)}, \mbox{ }
u:\Delta\rightarrow \left(\R^{4},J\right) \mbox{  $J$-holomorphic nonconstant},  u\left(0\right)=0 \right\}.$$
The point $0$ is a point of finite (D'Angelo) type $2m$ if $\Delta^1\left(\partial D,0\right)=2m<+\infty$.
\item The   regular type of $\partial D$ at origin is defined by:
\begin{eqnarray*}
\Delta^1_{{\rm reg}}\left(\partial D,0\right)&:=&\sup \{\delta_0\left(\partial D,u\right), \mbox{ }
u:\Delta\rightarrow \left(\R^{4},J\right) \mbox{  $J$-holomorphic },\\
& & \hspace{6,5cm} u\left(0\right)=0, d_0u\neq 0 \}.
\end{eqnarray*}
\end{enumerate}
\end{defi}
The  type condition  as defined in part 1 of Definition \ref{deftyp} was introduced by
J.-P.D'Angelo \cite{d'a2} who proved that this coincides  with the regular type in complex manifolds
of dimension two.  It was proved in \cite{ber} that the (D'Angelo) type and  the regular type coincide in
four dimensional almost complex manifolds.

\vspace{0,5cm}

In the next proposition, we describe locally the almost complex structure $J$ and the  defining function $\rho$
(see \cite{ber} for a proof).
\begin{prop}\label{prop1}
Let $D=\{\rho<0\}$ is a smooth domain in $\R^4$. Assume that $\rho$ is $J$-plurisubharmonic on a neighborhood of $\overline{D}$ where the structure $J$ is defined on a fixed neighborhood $U$ of $\overline{D}$.
We suppose that the origin $0 \in partial D$ is a point of finite type $2m$.
Then there is a local change of coordinates in a neighborhood of the origin such that, in the new coordinates:
\begin{equation}\label{eqstr}
J=\left(\begin{array}{ccccc}
a_1& b_1 & 0 & 0 \\
c_1 & -a_1 & 0 & 0\\
0  & 0 & a_2 & b_2\\
0 & 0 & c_2 & -a_2\\
\end{array}\right)=J_{st}+O\left(|z_2|\right),
\end{equation}

\begin{equation}\label{eqform}
\rho=\Re ez_2+H_{2m}\left(z_1,\overline{z_1}\right)+\widetilde{H}(z_1,z_2)+O\left(|z_1|^{2m+1}+|z_2||z_1|^m+|z_2|^2\right)
\end{equation}
where  $H_{2m}$ is a homogeneous polynomial of degree $2m$, subharmonic which  is not harmonic and
$$\widetilde{H}(z_1,z_2)=\Re e \displaystyle \sum_{k=1}^{m-1} \rho_{k} z_1^k z_2.$$
\end{prop}

A crucial tool for the study of pseudoholomorphic curves into pseudoconvex domains of finite type is the local peak $J$-plurisubharmonic functions which existence was proved in \cite{ber}.
\begin{theo}\label{theoex}
Let $D=\{\rho<0\}$ be a domain of finite type in an almost complex
manifold $\left(M,J\right)$ of  dimension four. We suppose that
$\rho$ is a smooth defining function of $D$, $J$-plurisubharmonic on a neighborhood of $\overline{D}$.
Let $p \in \partial D$ be a boundary point. Then  there exist a function $\varphi_p$ and a neighborhood $U$ of $p$ such that $\varphi_p$ is
continuous up to $\overline{D}\cap U$ and satisfies:
\begin{enumerate}
\item $\varphi_p$ is $J$-plurisubharmonic on  $D\cap U$,
\item $\varphi_p\left(p\right)=0$,
\item $\varphi_p <0$ on  $\overline{D}\cap U\backslash\{p\}.$
\end{enumerate}
Such a function is called a local peak $J$-plurisubharmonic function at $p$.
\end{theo}

\subsection{The Kobayashi pseudometric}
The existence of local pseudoholomorphic discs proved in \cite{ni-wo}
allows to define the {\it Kobayashi pseudometric} $K_{\left(M,J\right)}$ for $p\in M$ and $v \in T_pM$ :
$$K_{\left(M,J\right)}\left(p,v\right):=\inf
\left\{\frac{1}{r}>0, u : \Delta \rightarrow \left(M,J\right)
\mbox{  $J$-holomorphic }, u\left(0\right)=p, d_{0}u\left(\partial/\partial x\right)=rv\right\}.$$

Let $d_{\left(M,J\right)}$ be the integrated pseudodistance of $K_{\left(M,J\right)}$ defined by:
$$d_{\left(M,J\right)}\left(p,q\right): =\inf\left\{\int_0^1
K_{\left(M,J\right)}\left(\gamma\left(t\right),\dot{\gamma}\left(t\right)\right)dt, \mbox{ }
\gamma : [0,1]\rightarrow M, \mbox{ }\gamma\left(0\right)=p, \gamma\left(1\right)=q\right\}.$$

We thus define :
\begin{defi}\label{defin}
\begin{enumerate}
\item The manifold $\left(M,J\right)$ is Kobayashi hyperbolic if the integrated pseudodistance $d_{\left(M,J\right)}$ is a distance.
\item Let $D\subset M$ be a domain in an almost complex manifold $\left(M,J\right)$.
A point $p \in \partial D$ is said to be at finite distance from $q \in D$ if there is a sequence of points $q_j \in D$ converging to $p$ and whose Kobayashi distances $d_{\left(D,J\right)}(q_j,q)$ to $q$ stay bounded.
Otherwise we say that the distance is infinite.
\end{enumerate}
\end{defi}

\vspace{0,7cm}

\section{Proof of Theorem \ref{theo1}}

In order to prove this theorem, we need the two following lemmas, where the dimension assumption is meaningful as it allows to find a coordinate system where the lines $\{z_1=c\}$ and $\{z_2=c'\}$ are almost complex submanifolds.

\begin{lem}\label{lemgrad}
Let $\Omega$ be an open subset
of $(\R^{4},J)$ and let $J$ be an almost complex structure satisfying (\ref{eqstr}).
Let $K$ be a compact subset of $\Omega$. There exists $\delta>0$
such that: for every $r\in [0,1)$
there exists a positive constant $C>0$
such that if $u=(u_1,u_2):\Delta \to \Omega$ is a $J$-holomorphic disc
with $u(\Delta)\subset K$, then
$$
|\nabla u_1(\zeta)| \leq C {\rm sup}_{t\in \Delta}|u_1(t)-u_1(0)|, \mbox{ and, } |\nabla u_2(\zeta)| \leq C {\rm sup}_{t\in \Delta}|u_2(t)-u_2(0)|,
$$
if ${\rm sup}_{t\in \Delta}\|u(t)-u(0)\|\leq \delta$ and
$\zeta \in \Delta_r$.
\end{lem}
This lemma is an anisotropic version of a result obtained by S.Ivashkovich, J.-P.Rosay in \cite{iv-ro}.
\begin{proof}
Depending on $u(0) \in K$, one can make a linear change of variables such that in the new coordinates
$J(u(0)) = J_{st}$. Set $\alpha_i:={\rm sup}_{t \in \Delta}~|u_i(t)-u_i(0)|$ for $i=1,2$
and consider the scaling map $\Lambda$ from $\R^4$ into itself defined
by  $\Lambda(z_1,z_2):=(\alpha_1^{-1} (z_1-u_1(0)),\alpha_2^{-1} (z_2-u_2(0)))$.
Notice that the $\Lambda_*J$-holomorphic disc $\Lambda\circ u$ satisfies
$\Lambda\circ u (\Delta) \subseteq \Delta\times \Delta$.
If $\alpha_1$ and $\alpha_2$ are small enough  then $\Lambda_*J$ is close to $J_{st}$ as the structure $J$ has a diagonal form.
It follows from Proposition 2.3.6 in \cite{jcsi}, that for $|z|\leq r$ one gets for $i=1,2$, $|\nabla (\Lambda \circ u)_i(z)|\leq C$
for some positive constant $C$. Hence $|\nabla u_i(z)|\leq C\alpha_i$, as desired.
\end{proof}

A straightforward computation leads to this very useful lemma:
\begin{lem}\label{lemlev}
Assume that $J$ is an almost complex structure on $\R^4$ satisfying (\ref{eqstr}).
Then the Levi form of $\Re e z_2$ at  $z$ and $v=\left(X_1,Y_1,X_2,Y_2\right)$ is equal to
\begin{eqnarray*}
\mathcal{L}_{J}\Re e z_2\left(z,v\right)&=& \left[\left(a_1-a_2\right)(z)\frac{\partial a_2}{\partial x_1}(z)
-c_2(z)\frac{\partial b_2}{\partial x_1}(z)+c_1(z)\frac{\partial a_2}{\partial y_1}(z)\right]X_1X_2+\\
&&\\
& &\left[\left(a_1+a_2\right)(z)\frac{\partial b_2}{\partial x_1}(z)
-b_2(z)\frac{\partial a_2}{\partial x_1}(z)+c_1(z)\frac{\partial b_2}{\partial y_1}(z)
\right]X_1Y_2+\\
&&\\
& &\left[-\left(a_1+a_2\right)(z)\frac{\partial a_2}{\partial y_1}(z)+
b_1(z)\frac{\partial a_2}{\partial x_1}(z)-c_2(z)\frac{\partial b_2}{\partial y_1}(z)
\right]Y_1X_2+\\
&&\\
& &\left[\left(a_2-a_1\right)(z)\frac{\partial b_2}{\partial y_1}(z)+
b_1(z)\frac{\partial b_2}{\partial x_1}(z)-b_2(z)\frac{\partial a_2}{\partial y_1}(z)
\right]Y_1Y_2+\\
&&\\
& & \left[\frac{\partial a_2}{\partial y_2}(z)-\frac{\partial b_2}{\partial x_2}(z)\right]\left(c_2(z)X_2^2-2a_2(z)X_2Y_2 -b_2(z)Y_2^2\right).
\end{eqnarray*}

\end{lem}

\vspace{1cm}

\begin{proof}[Proof of Theorem \ref{theo1}.]Let $U$ be a neighborhood of $0$ in $\R^4$.
We assume that, on $U$, the structure $J$ satisfies (\ref{eqstr}) and  that the defining function has the local expression
\begin{equation}\label{eqform2}
\rho=\Re e z_2+H_{2m}\left(z_1,\overline{z_1}\right)+O\left(|z_1|^{2m+1}+|z_2||z_1|^m+|z_2|^2\right)
\end{equation}
where  $H_{2m}$ is a homogeneous polynomial of degree $2m$, subharmonic which  is not harmonic.

\vspace{0.8cm}

Consider for a positive number $\delta>0$, the following anisotropic polydisc:
$$Q(0,\delta):=\{z \in \C^2, |z_1|<\delta^{\frac{1}{2m}}, |z_2|<\delta \}.$$
Notice that since the defining function $\rho$ satisfies (\ref{eqform2}), then for a sufficiently small
$\delta<1$, if $z \in Q(0,\delta)$ then we have ${\rm dist}\left(z,\partial D\right)\leq c\delta$
for some positive constant $c>0$.

\vspace{0.8cm}

Let $q'=(q_1',q_2') \in \partial D \cap U$ be a boundary point and let $\varphi_{q'}$ be a local peak $J$-plurisubharmonic function at
 the point $q'$. There is a positive constant $C_1$ such that
\begin{equation}\label{eqine}
-C_1\|z-q'\|\leq \varphi_{q'}\left(z\right) \leq -C_2\Psi_{q'}\left(z\right),
\end{equation}
where
$$\Psi_{q'}\left(z\right):=|z_1-q_1'|^{2m}+|z_2-q_2'|^2+|z_1-q_1'|^2|z_2-q_2'|^2$$
is a $J$-plurisubharmonic function on  $U$, shrinking $U$ if necessary.

\vspace{0.8cm}

Let  $u~: \Delta \rightarrow D\cap U$ be a $J$-holomorphic disc such that $u\left(0\right) \in Q(0,\delta)$
is sufficiently  close to the origin. In order to prove that the origin $0$ is at infinite distance
from points in $D\cap U$, we want to provide the following estimates
\begin{equation*}
|\nabla u_1\left(0\right)|\leq C\delta^{\frac{1}{2m}} \mbox{ and } \mbox{ } |\nabla u_2\left(0\right)|\leq C\delta
\end{equation*} for a positive constant
$C>0$.

\vspace{0.5cm}

Let $q'=(q_1',q_2') \in \partial D$ be the unique boundary point such that $q'=u(0)+(0,\delta_u)$ for some positive $\delta_u$. Notice that $\delta_u$ is asymptotically equivalent to ${\rm dist}\left(u\left(0\right),\partial D\right)$.
According to the $J$-plurisubharmonicity of $\Psi_{q'}$,  we have  for $|\zeta|< r$ where $0<r<1$:
$$\Psi_{q'}\left(u\left(\zeta\right)\right)\leq \frac{A}{2\pi}\displaystyle
\int_0^{2\pi}\Psi_{q'}\left(u\left(re^{i\theta}\right)\right)d\theta,$$
for an appropriate positive constant $A$.
Hence using (\ref{eqine}) and the $J$-plurisubharmonicity of the peak function $\varphi_{q'}$ we obtain:
$$\Psi_{q'}\left(u\left(\zeta\right)\right)\leq -\frac{A}{2\pi C_2}\displaystyle \int_0^{2\pi} \varphi_{q'}\left(u\left(re^{i\theta}\right)\right)d\theta
\leq -\frac{A}{C_2}\varphi_{q'}\left(u\left(0\right)\right).$$
Since
$$|u_2\left(\zeta\right)-q_2'|^{2}\leq \Psi_{q'}\left(u\left(\zeta\right)\right)$$
and according to (\ref{eqine}) and to the fact that $\|u(0)-q'\|=\delta_u$ is asymptotically equivalent to ${\rm dist}\left(u\left(0\right),\partial D\right)$,
we  obtain for a positive constant $C_3$:
$$|u_2\left(\zeta\right)-q_2'|^{2}\leq C_3{\rm dist}\left(u\left(0\right),\partial D\right).$$
Hence, for some other positive constant $C_3$ we have:
\begin{equation}\label{eqnorm1}
|u_2\left(\zeta\right)-u_2(0)|\leq C_3{\rm dist}\left(u\left(0\right),\partial D\right)^{\frac{1}{2}},
\end{equation}
for $|\zeta|< r$ where $0<r<1$.

For the same reason, we obtain the following estimate which hold for $|\zeta|< r<1$:
\begin{equation}\label{eqtan1}
|u_1\left(\zeta\right)-u_1(0)|\leq C_3{\rm dist}\left(u\left(0\right),\partial D\right)^{\frac{1}{2m}}.
\end{equation}

According to Lemma \ref{lemgrad}, inequalities (\ref{eqnorm1}) and (\ref{eqtan1}) imply:
\begin{equation}\label{eqnab}
|\nabla u_2(\zeta)|\leq B \sup_{|t|<r} |u_2\left(t\right)-u_2\left(0\right)|\leq
C_3B{\rm dist}\left(u\left(0\right),\partial D\right)^{\frac{1}{2}} \leq C_4\delta^{\frac{1}{2}}.
\end{equation}
and,
\begin{equation}\label{eqfin1}
|\nabla u_1(0)|\leq B \sup_{|t|<r} |u_1\left(t\right)-u_1\left(0\right)|\leq
C_3B{\rm dist}\left(u\left(0\right),\partial D\right)^{\frac{1}{2m}} \leq C_4\delta^{\frac{1}{2m}}.
\end{equation}
for  positive constants  $B$ and $C_4$. Notice that (\ref{eqfin1}) is the desired tangential estimate.

\vspace{0.5cm}

In order to obtain the normal estimate $|\nabla u_2(0)|$, we will construct a negative harmonic function.
To achieve this, we first need to control $\Re e u_2(\zeta)$ and $|\Delta {\Re e}~u_2(\zeta)|$ by $\delta$.

From  (\ref{eqnorm1}) and (\ref{eqtan1}) we obtain for $|\zeta|< r$:
\begin{equation}\label{eqdist3}
|u_1(\zeta)|\leq |u_1\left(\zeta\right)-u_1(0)|+|u_1(0)| \leq  C_3{\rm dist}\left(u\left(0\right),\partial D\right)^{\frac{1}{2m}}+\delta^{\frac{1}{2m}} \leq C_5\delta^{\frac{1}{2m}}.
\end{equation}
and,
\begin{equation}\label{eqdist4}
|u_2(\zeta)|\leq |u_2\left(\zeta\right)-u_2(0)|+|u_2(0)| \leq C_3{\rm dist}\left(u\left(0\right),\partial D\right)^{\frac{1}{2}} + \delta \leq C_{5}\delta^{\frac{1}{2}}.
\end{equation}
for a positive constant $C_{5}$. According to these estimates, and since the defining function satisfies
(\ref{eqform2}), we obtain
\begin{equation}\label{eqdist5}
\Re e  u_2(\zeta)\leq \left|H_{2m}\left(u_1(\zeta),\overline{u_1(\zeta)}\right)\right|+
O\left(|u_1(\zeta)|^{2m+1}+|u_2(\zeta)||u_1(\zeta)|^m+|u_2(\zeta)|^2\right)\leq C_{6}\delta.
\end{equation}
if $|\zeta|<r$, where $C_{6}>0$.
Due to Proposition \ref{proplevi} applied to the function $\Re e z_2$, we have:
\begin{equation}
\Delta {\Re e}~u_2(\zeta)=\mathcal{L}_{J}\Re e z_2\left(u(\zeta),\frac{\partial u}{\partial x}(\zeta)\right).
\end{equation}
According to Lemma \ref{lemlev}, (\ref{eqnab}), (\ref{eqdist4}) and to the fact that $J(z)=J_{st}+O(|z_2|)$, this implies
\begin{equation}\label{eqlap2}
|\Delta {\Re e}~u_2(\zeta)|\le C_{7}\delta
\end{equation}if $|\zeta|<r$, and a positive constant $C_{7}$.

We still denote by $\Delta {\Re e}~u_2$, the function equal to $\Delta {\Re e}~u_2$ on $\Delta_{r}$,
and $0$ elsewhere. And consider the following harmonic function
\begin{equation}\label{eqg2}
g(\zeta):=\Re e u_2(\zeta)-\left[\frac{1}{2\pi}\ln |\zeta|*\Delta {\Re e}~u_2\right](\zeta)
-(C_{6}+C_{7})\delta.
\end{equation}
From (\ref{eqlap2}), we have
$$\left|\frac{1}{2\pi}\ln |\zeta|*\Delta {\Re e}~u_2~\right|\leq C_{7}\delta.$$
Thus $g$ is a negative function and $|g(0)|\leq 2(C_{6}+C_{7})\delta$.
The classical Schwarz Lemma for negative harmonic functions gives
$|\nabla g(0)|\leq 2 |g(0)|$. Hence we obtain the following estimate:
$$|\nabla {\Re e}~u_2(0)|\leq |\nabla g(0)|+{\rm Sup}~|\Delta {\Re e}~u_2|
\leq 2|g(0)|+C_7\delta\leq (4C_{6}+5C_{7})\delta.$$

Moreover, the $J$-holomorphy equation for the disc $u$ implies:
\begin{eqnarray*}
a_2(u(0))\frac{\partial \Re e u_2}{\partial x}(0)+b_2(u(0))\frac{\partial \Im m u_2}{\partial x}(0)&=&\frac{\partial \Re e u_2}{\partial y}(0)\\
c_2(u(0))\frac{\partial \Re e u_2}{\partial x}(0)-a_2(u(0))\frac{\partial \Im m u_2}{\partial x}(0)&=&\frac{\partial \Im m u_2}{\partial y}(0),\\
\end{eqnarray*}
so that
$$|\nabla {\Im m}~u_2(0)|\leq C_{8} |\nabla {\Re e}~u_2(0)|,$$
for some positive constant $C_{8}$.

We finally get  the  normal estimate:
\begin{equation}\label{eqfin2}
|\nabla u_2(0)|\leq C\delta.
\end{equation}

\vspace{0.8cm}

Let us show how  (\ref{eqfin1}) and (\ref{eqfin2}) imply that the origin $0$ is at infinite distance from points in $D\cap U$. Consider the function:
$$\chi(z):=|z_1|^{4m}+|z_2|^2.$$
Since  $u(0) \in Q(0, \chi(u(0))^{\frac{1}{2}})$, it follows from the estimates
(\ref{eqfin1}) and  (\ref{eqfin2}) that
$|\nabla u_1(0)|\leq C\chi(u(0))^{\frac{1}{4m}}$ and $|\nabla u_2(0)|\leq C\chi(u(0))^{\frac{1}{2}}$ for some positive constant $C$ and thus:
\begin{eqnarray*}
|\nabla (\chi\circ u )(0)| &\leq & 4m|\nabla u_1(0)||u_1(0)|^{4m-1}+2|\nabla u_2(0)||u_2(0)|\\
&\leq & 4mC\chi(u(0))^{\frac{1}{4m}}(\chi(u(0))^{\frac{1}{4m}})^{4m-1}+2C\chi(u(0))^{\frac{1}{2}}\chi(u(0))^{\frac{1}{2}}\\
&\leq & (4mC+2C)(\chi\circ u )(0)\\
\end{eqnarray*}
which achieves the proof of Theorem \ref{theo1} by an integration argument (see Lemma 1.1 in \cite{iv-ro}).

\end{proof}


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

The nature and origin of the dark matter in the universe and the origin
of ultrahigh energy cosmic rays (UHECR) are certainly two of the most
interesting problems of astroparticle physics, and maybe even of physics,
at the beginning of the 21st century.

The three puzzles related to the origin of cosmic rays with
energies beyond the Greisen-Zatsepin-Kuzmin cutoff,
$E>E_{GZK}\simeq 4\times 10^{19}\,$eV, are:\\
{\bf 1.} Scattering off cosmic microwave
background photons limits the penetration depths of charged
particles at these energies to distances $<100\,$Mpc
\cite{GZK,stecker,new};\\
{\bf 2.} the distribution of arrival directions of UHECRs
does not seem to favor any known astrophysical sources within
the GZK cutoff length;\\
{\bf 3.} it seems extremely difficult to devise sufficiently efficient
astrophysical acceleration mechanisms which could accelerate particles
to energies $E>E_{GKZ}$, and at the same time overcome collisional
and radiation losses.

The so called top-down models of UHECRs combine both problems by proposing that
ultrahigh energy cosmic rays arise in the decay \cite{Hill,BKV,KR,BS,DFK} or
collisional annihilation \cite{BDK,DBK,DHW} of superheavy dark matter
particles. We will use both the acronym SHDM and
{\sc wimpzilla}\footnote{The word {\sc wimpzilla} was coined by Kolb,
Chung and Riotto in their investigations of possible origins of
superheavy dark matter particles \cite{rocky}.}
for superheavy dark matter particles.

SHDM decay can arise through direct decay of relic particles
(``{\sc wimpzilla}" decay) or through the annihilation of superheavy relic
bound states (``{\sc wimpzillium}" decay), but the underlying decay mechanism
plays no role for the predicted pattern of UHECR arrival directions.
This pattern differs strongly from the anisotropy pattern predicted by
collisional {\sc wimpzilla} annihilation. The reason for the
difference is that the anisotropy pattern predicted by decay models
is proportional to the {\sc wimpzilla} or {\sc wimpzillium} density
$n_X(\vek{r})$, and is dominated by the smooth background halo.
Collisional annihilation, on the other hand, predicts that the anisotropy
pattern should be proportional to $n_X^2(\vek{r})$, and is
constrained by unitarity limits on annihilation cross sections.
Therefore collisional annihilation can only work in dense cores of dark
matter substructure in the galactic halo. These dense cores were denoted as
{\sc wimpzilla} stars \cite{DHW}.
As a consequence, the collisional annihilation scenario predicts a pointlike
source distribution with increasing density towards the galactic center.

All modern acceleration (``bottom-up'')
models for UHECR origin assume powerful extragalactic sources, and therefore
predict that the pattern of observed arrival directions should not correlate
with the galactic halo. The different patterns of UHECR arrival patterns
predicted by all the different source models imply that a dedicated
large statistics experiment like the Pierre Auger observatory can easily
identify the correct model from its anisotropy signal.

Sec. \ref{sec:shdm} explains the origin of a mass window for the
generation of relic superheavy dark matter particles during inflation.
Sec. \ref{sec:flux} summarizes and updates the calculation of the
ultrahigh energy cosmic ray flux from collisional {\sc wimpzilla}
annihilation, and the origin of the unique pattern of arrival directions
predicted by that model. In Sec. \ref{sec:pattern} we compare the
different anisotropy signatures predicted by all contemporary proposals
for the origin of ultrahigh energy cosmic rays, and Sec. \ref{sec:conc}
contains our conclusions.

\section{Superheavy dark matter from inflation}\label{sec:shdm}

Particle creation as a consequence of non-adiabatic expansion was
discovered already in the late 1950s and 1960s, see \cite{parker} and
references there. However, the effect was found to be negligible in
radiation or dust dominated epochs, and therefore this mechanism for
non-thermal particle creation was rediscovered and garnered much more
interest only after the necessity for inflation was realized.
The effect is usually considered in terms of the Bogolubov transformation
between in and out vacua in an expanding universe
\cite{rocky,CKR,KT,GPRT,CCKR}. Particle creation during
preheating after inflation can also arise as a consequence of a direct
coupling between the inflaton and other matter fields \cite{lev1,lev2}.

Here we will expand on a simple discussion given in \cite{DHW} to see
how particle production in an inflationary universe can
be understood by studying the evolution equations of weakly coupled
scalar fields in the expanding universe. The reasoning outlined here is
not intended to compete in any way with the traditional operator methods
to study particle production from spatial expansion, but it may provide a
complementary and helpful view.

If interactions with other matter fields can be neglected, a scalar field
in a Friedmann--Robertson--Walker background with metric
\begin{equation}\label{eq:frw}
ds^2=-dt^2+a^2(t)\!\left(
\frac{dr^2}{1-kr^2}+r^2d\vartheta^2+r^2\sin^2\vartheta\,d\varphi^2\right)
=-dt^2+a^2(t)\!\left(
d\vek{x}^2+k\frac{(\vek{x}\cdot d\vek{x})^2}{1-k\vek{x}^2}\right)
\end{equation}
satisfies
\begin{equation}\label{eq:phimot}
\ddot{\phi}(\vek{x},t)+3\frac{\dot{a}(t)}{a(t)}\dot{\phi}(\vek{x},t)
-\frac{1}{a^2(t)}\!\left[\left(\delta^{ij}-k x^i x^j\right)
\partial_i\partial_j-3kx^i\partial_i\right]\phi(\vek{x},t)
+m^2\phi(\vek{x},t)=0,
\end{equation}
and the corresponding energy density per unit of comoving volume is
\begin{eqnarray}\label{eq:rho}
\varrho(\vek{x},t)&=&\sqrt{-g}\,T^{00}
\\ \nonumber
&=&\frac{a^3(t)}{2\sqrt{1-kr^2}}\!\left[
\dot{\phi}^2(\vek{x},t)+m^2\phi^2(\vek{x},t)
+\frac{1}{a^2(t)}\!\left(\delta^{ij}-k x^i x^j\right)
\partial_i\phi(\vek{x},t)\cdot\partial_j\phi(\vek{x},t)\right].
\end{eqnarray}
The violation of time translation invariance in (\ref{eq:frw})
implies violation of energy conservation, of course:
\begin{eqnarray}\label{eq:rhodot}
\dot{\varrho}(\vek{x},t)&-&a(t)\partial_i\!\left[\dot{\phi}(\vek{x},t)
\frac{\delta^{ij}-k x^i x^j}{\sqrt{1-kr^2}}\partial_j\phi(\vek{x},t)\right]
\\ \nonumber
&=&\frac{a^2(t)\dot{a}(t)}{2\sqrt{1-kr^2}}\!\left[
3\!\left(m^2\phi^2(\vek{x},t)-\dot{\phi}^2(\vek{x},t)\right)
+\frac{1}{a^2(t)}\!\left(\delta^{ij}-k x^i x^j\right)
\partial_i\phi(\vek{x},t)\cdot\partial_j\phi(\vek{x},t)\right].
\end{eqnarray}
{\it A priori} the sign of the expression on the right hand side
is indefinite, e.g. a rapidly evolving field of low mass and with
small spatial fluctuations loses energy during spatial expansion
as long as it evolves rapidly enough. However, Eq. (\ref{eq:phimot})
tells us that spatial fluctuations (which should never have been large
anyway for the Friedmann--Robertson--Walker {\it ansatz} to work) are
soon negligible in the inflationary expanding universe
$a(t)\propto\exp(Ht)$.
This leaves us with the simple equation
\begin{equation}\label{eq:phi2}
\ddot{\phi}(t)+3H\dot{\phi}(t)
+m^2\phi(t)\simeq 0.
\end{equation}
Approximately constant $H$ during inflation then yields for
the time evolution of the comoving energy density
\begin{eqnarray*}
\varrho(t)&\simeq&
\frac{1}{2}a^3(t)\left(\dot{\phi}^2(t)
+m^2\phi^2(t)\right)
\\
&\simeq&
A_+\exp\!\left(t\sqrt{9H^2
-4m^2}\right)+A_-\exp\!\left(-t\sqrt{9H^2
-4m^2}\right)+B.
\end{eqnarray*}
This implies a growing mode in the comoving energy density of weakly coupled
states with $m<1.5H\simeq 10^{14}\,$GeV.
What is special about the superheavy particles is that their
comoving energy density is conserved after inflation, because the
behavior of massive $(m>t^{-1})$ weakly coupled states
in the subsequent radiation and dust dominated backgrounds preserves
their energy.
The asymptotic solution for weakly coupled massive states
with $m>t^{-1}$ in such a background yields (with $\ell=3$ for dust,
$\ell=4$ for radiation)
\[
\phi(t)\propto
t^{-3/\ell}\cos\!\left(mt+\varphi\right),
\]
\[
\varrho(t)\propto a^3(t)t^{-6/\ell}\propto t^0,
\]
and this implies in particular that the comoving density of massive
particles freezes out
at the end of inflation $(t\simeq 10^{-36}\,\mbox{s})$ if
\[
m>t^{-1}\simeq 10^{12}\,\mbox{GeV}.
\]
These considerations indicate a mass window for direct gravitational production
of superheavy relic particles during inflation
\[
10^{12}\,\mbox{GeV}<m<10^{14}\,\mbox{GeV}.
\]

After inflation intrinsically unstable matter of mass $M$ will
decay on time scales $\tau\lesssim M_{Planck}^2/M^3$, where the upper
bound assumes that the particles couple to their decay products at least
with gravitational strength. This means that {\it intrinsically unstable}
superheavy particles should not be relic particles. This difficulty
has motivated the collisional annihilation scenario for ultrahigh energy
cosmic rays from superheavy dark matter \cite{BDK}, because the unitarity
limit on reaction cross sections
\begin{equation}\label{eq:sigma1}
\langle\sigma_A v\rangle\lesssim\frac{4\pi}{M^2v}
\end{equation}
implies that superheavy dark matter without direct decay channels
will still be around.

\section{Calculation of the flux}\label{sec:flux}

We consider decay or annihilation of superheavy
dark matter particles of mass $M_X\ge 10^{12}\,$GeV.
The spectral fluxes at our location $\vek{r}_\odot$
from decay or collisional annihilation of the
dark matter particles of density $n_X(\vek{r})$ are then
\begin{equation}\label{eq:jdec}
j_d(E)=\frac{d\mathcal{N}(E,M_X)}{dE}
\int d^3\vek{r}\,\frac{1}{4\pi|\vek{r}_{\odot}-\vek{r}|^2}
\frac{n_X(\vek{r})}{\tau_d}
\end{equation}
and
\begin{equation}\label{eq:jann}
j_a(E)=\frac{d\mathcal{N}(E,2M_X)}{dE}
\int d^3\vek{r}\,\frac{\nu}{16\pi|\vek{r}_{\odot}-\vek{r}|^2}
n_X(\vek{r})^2 \langle\sigma_A v\rangle,
\end{equation}
respectively. Here
$d\mathcal{N}(E,E_{in})$ is the number of particles
in the energy interval $[E,E+dE]$ emerging from a decay or annihilation
event of initial energy $E_{in}$.
$d\mathcal{N}(E,E_{in})/dE$ is related to fragmentation
functions via
\[
\frac{d\mathcal{N}(E,E_{in})}{dE}=\sum_i\frac{1}{\sigma_A}
\frac{d\sigma^{(i)}}{dE}
=\frac{1}{E_{in}}\sum_i F^{(i)}(x,E_{in}^2),
\]
where $x=E/E_{in}$ and $F^{(i)}(x,E_{in}^2)$ is the differential
number of particles of species $i$ generated in the prescribed
$x$-range in a decay or annihilation event with $s=E_{in}^2$.
The factor $\nu$ in Eq. (\ref{eq:jann}) equals 4 if the $X$-particles
are Majorana particles, and 1 otherwise.

$N$-body simulations predict that usually about $f_{cl}\sim 5-10$\% of
dark matter halos should exist in substructure. Therefore in the decay
models the ultrahigh energy cosmic ray flux will be dominated by the
smooth background halo. Evaluating (\ref{eq:jdec}) e.g. for a
Navarro-Frenk-White halo \cite{NFW}
\[
n_X(r)=4n_X(r_s)r_s^3/r(r+r_s)^2
\]
yields
\begin{equation}\label{eq:jdecnfw}
j_d(E)=4\frac{d\mathcal{N}(E,M_X/2)}{dE}
\frac{n_X(r_s)}{\tau_d}\frac{r_s^3}{r_s^2-r_{\odot}^2}
\ln\!\left(\frac{r_s}{r_{\odot}}\right).
\end{equation}
The scale radius $r_s$ for the Milky Way is not well known,
but for the whole range of say $5\,\mathrm{kpc}\le r_s\le
50\,\mathrm{kpc}$ comparison between (\ref{eq:jdecnfw})
and the ultrahigh energy cosmic ray flux from AGASA \cite{agasa}
yields that decaying superheavy dark matter should only
make a small contribution to the galactic dark matter
density. Of course, the problem is to explain a lifetime
$\tau_d\ge 10^{17}\,$s for unstable superheavy dark matter.

Eq. (\ref{eq:jann}) can also be evaluated analytically for
a Navarro-Frenk-White halo, but the unitarity limit for
s-wave annihilation
\begin{equation}\label{eq:sigma}
\langle\sigma_A v\rangle=\xi\times\frac{4\pi}{M_X^2v}
=\xi\times 4.40\times 10^{-43}\,\mbox{m}^3/\mbox{s}
\times \left(\frac{10^{12}\,\mathrm{GeV}}{M_X}\right)^2
\times \frac{100\,\mathrm{km/s}}{v},\quad
\xi\le 1,
\end{equation}
implies that any ultrahigh energy cosmic ray flux from
collisional annihilation in the background halo is
negligible. Therefore collisional {\sc wimpzilla} annihilation
can only make a noticeable contribution to the ultrahigh energy
cosmic ray flux if it originates in dense cores of
dark matter subclumps of the galactic halo \cite{BDK}.
These dense cores were denoted as {\sc wimpzilla} stars
in \cite{DHW}, and a simplified estimate of the cosmic ray flux
from {\sc wimpzilla} stars yields
\begin{equation}\label{eq:flux1}
j_a(E)\simeq\nu\frac{N_{cl}\overline{V}_{core}}{16\pi d^2}
\frac{d\mathcal{N}(E,2M_X)}{dE}
n_X^2 \langle\sigma_A v\rangle
\simeq\nu\frac{0.1f_{cl}M_{halo}}{16\pi d^2 M_X}
\frac{d\mathcal{N}(E,2M_X)}{dE}
n_X \langle\sigma_A v\rangle.
\end{equation}
Here $d^{-2}=\langle r^{-2}\rangle$ is a mean inverse distance squared
for our separation from galactic {\sc wimpzilla} stars. In our estimate
in \cite{DHW} we had used $d\simeq 10\,$kpc, which seems reasonable:
e.g. our inverse mean distance squared to visible galactic substructures,
the globular clusters, corresponds to $d=7.3\,$kpc.

We have parametrized the unknown annihilation cross section already in
terms of the unitarity bound (\ref{eq:sigma}). Now we would like to
compare e.g. with the UHE cosmic ray spectrum in Ref. \cite{BDK},
Fig. 1 (note that the cross section stated there would only be required
for collisional annihilation in the smooth background halo
and does not apply to annihilation in dark matter subclump cores).
For the comparison we parametrize the core density $n_X$ of
{\sc wimpzilla} stars in terms of the solar density
\[
n_X=\eta\times\frac{\varrho_\odot}{10^{12}\,\mbox{GeV}}
=\eta\times 7.89\times 10^{17}\,\mbox{m}^{-3}.
\]
This yields a flux at $10^{11}\,$GeV of order
\begin{equation}\label{eq:e3j}
E^3 j(E)\Big|_{E=10^{11}\,\mathrm{GeV}}
=5.74\times 10^{25}\,\mathrm{eV}^2\,\mathrm{m}^{-2}\,
\mathrm{s}^{-1}\times\frac{f_{cl}}{0.05}\times
\frac{M_{halo}}{2\times 10^{12}M_\odot}\times
\left(\frac{7.3\,\mathrm{kpc}}{d}\right)^2\times\nu
\times\frac{\eta\xi}{10^{-3}}.
\end{equation}
This complies with the spectral fit for the flux per solid angle
in Ref. \cite{BDK}, Fig. 1, if $\eta\xi<1$, as would be expected
for consistency of the model:
If the {\sc wimpzillas} are not Majorana particles ($\nu=1$)
the required value for the fit is
$\eta\xi\simeq 9.8\times 10^{-4}$, and otherwise the required
value would be $\eta\xi\simeq 2.5\times 10^{-4}$.

\section{Expected anisotropy patterns beyond $20\,$EeV}\label{sec:pattern}

Fits of fragmentation functions to the ultrahigh energy cosmic ray spectrum
indicate that in the collisional annihilation scenario cosmic rays above
$20\,$EeV should be dominated by the fragmentation products of {\sc wimpzilla}
annihilation \cite{BDK}. Therefore the expected anisotropy patterns discussed
in this section only apply in this energy range.

The purpose of this section is to emphasize that the anisotropy pattern
observed by the Pierre Auger observatory above $20\,$EeV will provide a
crucial direct test of the different proposals for the origin of ultrahigh
energy cosmic rays: Bottom-up acceleration of charged paricles
in AGNs (see e.g. \cite{RB}) or gamma ray bursts \cite{GRB},
ultrahigh energy neutrinos travelling over cosmological
distances and creating Z-bursts \cite{zburst},
{\sc wimpzilla} or {\sc wimpzillium} decay
\cite{Hill,BKV,KR,BS,DFK}, or collisional annihilation in {\sc wimpzilla}
stars \cite{BDK,DHW}.\\

\begin{tabular}{|l|l|}
Origin of UHECRs & Expected anisotropy pattern\\
 & \\ \hline & \\
Collisional annihilation in {\sc wimpzilla} stars & About 1000 pointlike
 sources with increasing\\
 & density towards the galactic center.\\
 & No correlation with galactic SNRs.\\
 & \\ \hline & \\
{\sc wimpzilla} or {\sc wimpzillium} decay &
Dominated by uniform increase towards galactic center. \\
 & \\ \hline & \\
Z-bursts & Approximately isotropic distribution with only weak\\
 & correlation to structure within $150\,$Mpc. \\
 & \\ \hline & \\
Bottom-up acceleration & Correlation with local superstructure
within $150\,$Mpc. \\ & \\ \hline
\end{tabular}
{}\\

Once fully operational, the Pierre Auger observatory should see more than
300 events per year with energies above $40\,$EeV. If the collisional
annihilation scenario is correct, the observatory should see a large number
of multiplets within its angular resolution, with increasing
density towards the galactic center.

\section{Conclusions}\label{sec:conc}

Collisional annihilation of superheavy dark matter particles as a source
model for ultrahigh energy cosmic rays has the advantage of naturally
explaining the absence of a GZK cutoff in the spectrum without correlation
to local AGNs, without the need of postulating an extremely powerful and
efficient acceleration mechanism, and without the need to explain
extremely long lifetimes of intrinsically unstable particles. It has
the disadvantage of having to postulate formation of a few relatively
dense cores of dark matter subhalos or subhalo remnants.

The two different top-down scenarios of decay or collisional annihilation
imply qualitatively different anisotropy signals in cosmic ray arrival
directions, and these also differ from the anisotropy signals predicted
by the various bottom-up scenarios. What constitutes a statistically
significant dataset to confirm or reject the different anisotropy
patterns depends on the different patterns, of course. However, for the
collisional annihilation scenario the Pierre Auger observatory should see
a large number of multiplets within its angular resolution already after
about two years of full operation, because the UHECRs should arise from a
limited number of pointlike sources in our galactic halo. Based on
current particle physics extrapolations for chemical composition of
cosmic rays to very high energies, bottom-up acceleration scenarios seem
much more popular in the cosmic ray community. However, collisional
{\sc wimpzilla} annihilation is not an overly speculative theory, and
based on its expected anisotropy pattern it will be confirmed or rejected
in a very short time frame. It would therefore seem prudent to await what
the Pierre Auger observatory will tell us about patterns of UHE cosmic ray
arrival directions before jumping to any conclusions about the validity
of different models.

redpajama_set_name: RedPajamaArXiv


text:
\section{INTRODUCTION} \label{sec1}

Radiotherapy is one of the most widely used modalities for palliative and therapeutic cancer care \cite{steelbook}.
Its performance depends on energy deposition, and the resulting effects on tumor and surrounding healthy tissue.
While the dose profile is controlled through the choice of radiation and delivery method \cite{radonc}, its effects on tissue can be manipulated using chemical radiosensitizers \cite{lehnert} and/or tumor heating \cite{hyperrt}.
The use of radiosensitizers enhances cell kill and allows for optimization of radiation dosage and scheduling such that the collateral damage to the normal tissue can be minimized.
Cytotoxic effects of heat shocks and their role in radiosensitization for radiotherapy have been of interest since the 1970s.
\vspace{-5pt}

There are mechanistic differences in how heat and radiation damage cells \cite{radbio}.
This allows for additive use of hyperthermia and radiotherapy.
The role of hyperthermia in radiochemosensitization, without the introduction of systemic toxicity, is also well known \cite{peeken}.
This allows for the synergistic use of hyperthermia and radiotherapy.
Therefore, hyperthermia can be used with radiotherapy both additively and synergistically.
Several randomized trials have demonstrated clinical benefits of using hyperthermia in combination with radiotherapy \cite{radbio}.
Hyperthermia is a fully approved part of cancer care in European countries like Germany and the Netherlands.
\vspace{-5pt}

Radiation dosimetry has traditionally been treated as single physics \cite{attix}.
However, the modeling of tumor response to radiotherapy could benefit from the coupling of dose deposition with other relevant physics.
When patients are exposed to heat and radiation, tumors undergo constitutive changes.
Taking these changes into account simultaneously with the radiation dose deposition, should provide a more realistic and effective approach to model treatment.
For example, hyperthermia increases radiosensitivity in certain types of tumor cells, making tumors more susceptible to radiation treatment \cite{hyperrt}.
Therapies such as Combined-Hyperthermia-Radiotherapy (CHR) make use of this effect.\vspace{-5pt}

To appropriately model the overall effect, the synergy between features such as tumor heating, radiation dosimetry, and tumor dynamics should be considered simultaneously.
Furthermore, there are cases in which the target organs are in constant motion (e.g.~lungs), thereby constantly shifting the location of the tumor (in center-of-mass frame).
This motion should be accounted for in real time while treating tumors using external radiation beams and mathematically modeling their response.
Coupled/multiphysics modeling could prove to be a powerful tool to determine tumor dynamics in such systems.
Other treatment modalities include a combination of radiotherapy with immunotherapy \cite{radimm}, and radiotherapy with chemotherapy \cite{radche}.
For such combined treatments, the accurate determination of tumor response should also be performed using multiphysics models.
\vspace{-5pt}

Several studies have presented work on mathematical tumor dynamics modeling.
Ordinary and partial differential equations are used to model these systems in \cite{altrock} and \cite{radmichor}.
Radiotherapy-specific tumor dynamics models have been introduced in \cite{rockne} and \cite{sachs}.
Moreover, internal radiation dosimetry has been extensively treated using the Monte Carlo method \cite{rogers}.
In this work, we introduce the idea of a novel multiphysics-based approach to mathematically model CHR systems.
This model couples radiation transport (and dosimetry), heat transfer, and cell population dynamics.
To our knowledge, this is the first attempt to develop a multiphysics model for tumor response to CHR.
\vspace{-5pt}

This is by no means an exhaustive study.
The test problems in this preliminary work \textit{do not} model real-world treatment scenarios; the governing equations and treatment specifics are heavily simplified.
These are not the subject of investigation in this report; our goal is to provide a first step towards building a fully-coupled multiphysics model for this type of treatment, with the focus of this paper being the introduction to the coupled approach.
\vspace{-5pt}

The remainder of this paper is organized as follows.
In \cref{sec:gov_eq}, we introduce relevant physics, governing equations, coupling, and numerical approach.
In \cref{sec3}, we furnish our example problems and provide numerical results.
We conclude this paper with a brief discussion in \cref{sec4}.
\vspace{-5pt}

\section{GOVERNING EQUATIONS}
\label{sec:gov_eq}

The model to describe tumor dynamics during CHR treatment incorporates three physics: 1) radiation transport to model internal dosimetry, 2) heat transfer to incorporate the effects of tumor heating, and 3) cell population dynamics to determine how the tumor responds to the combination of hyperthermia with radiotherapy.
Our preliminary model is one-dimensional in space, monoenergetic, and assumes isotropic scattering.
\vspace{-5pt}

\subsection{Radiation Transport and Dosimetry}\vspace{-5pt}

The following monoenergetic Boltzmann transport equation with isotropic scattering is used to model the flux distribution \cite{lm}:
\begin{equation}
\label{transport1}
\frac{1}{v}\frac{\partial \psi}{\partial t} + \mu\frac{\partial \psi}{\partial x} + \sigma_t(x) \psi =  \frac{\sigma_{s}(x)}{2} \phi(x, t) + \frac{Q(x, t)}{2},
\end{equation}
where $\psi = \psi(x,\mu, t)$ represents the angular flux at position $x$ and along direction $\mu$ at time $t$; $v$ is the photon speed; $\phi$ is the scalar flux; $\sigma_t$ is the total macroscopic cross-section; $ \sigma_{s}$ is the $0^{th}$ moment of the scattering cross-section; and $Q$ is an internal source.
The dose is calculated using the following equation \cite{martin}:
\begin{subequations}
\begin{equation}
\label{dose1}
D(x, t) = \left(\frac{\mu_{en}}{\rho}\right)\Phi(x, t) E,
\end{equation}
\begin{equation}
\label{fluence}
\Phi(x, t) = \int_{0}^{t} dt' \phi(x, t').
\end{equation}
\end{subequations}
Here, $D$ represents the dose; $\frac{\mu_{en}}{\rho}$ is the mass energy absorption coefficient; $E$ is the energy per photon; and $\Phi$ is the fluence at time $t$.
We calculate the effective dose $D_E$ using the standard linear-quadratic model \cite{radbio}:
\begin{equation}
\label{edose1}
D_E(x, t) = \alpha D(x, t) + \beta D^2(x,t),
\end{equation}
where $\alpha$ and $\beta$ are radiobiology parameters that determine relative contribution of each term in the sum toward the total radiation effect \cite{radbio}.
The cell survival probability $S$ is chosen in such a way that a larger dose results in smaller survival probability, as given by \cite{rockne}:
\begin{equation}
\label{survival}
S(x,t) = e^{-D_E(x, t)}.
\end{equation}

\subsection{Hyperthermia and Heat Transfer}\vspace{-5pt}

Thermal medicine is slowly increasing its presence in cancer care. Several techniques, such as radiofrequency and ultrasound ablation, have been developed to provide focused heating inside tumors \cite{radbio}.
While the use of a more sophisticated heat transport model is warranted, we use simple heat conduction to model hyperthermia \cite{cannon}:
\begin{equation}
\label{heat}
\rho c_p \frac{\partial T(x, t)}{\partial t} = \frac{\partial }{\partial x} \kappa \frac{\partial T(x, t)}{\partial x} + q(x),
\end{equation}
where $c_p$, $\rho$, and $\kappa$ are the specific heat capacity, the material density, and the thermal conductivity of tissue respectively (assumed constant).
Moreover, $q$ is the volumetric heat source.
\vspace{-5pt}

\subsection{Cell Population Dynamics}\vspace{-5pt}

We assume that tumor growth depends on three factors: proliferation, invasion, and cell kill induced by treatment.
Proliferation determines the balance between cell division and cell loss due to natural death; cell proliferation is higher in tumors.
Invasion essentially determines the transport of tumor cells in the area of interest; continuum-based approaches to model tumor dynamics often use Fick's law to describe invasion.
Two proliferation models are widely considered in the literature: exponential and logistic.
Though less suitable for problems considering realistic tumor growth time-frames (months to years), we use the exponential model because our test problems consider tumor dynamics over a short period of thirty minutes.
The following balance equation determines tumor cell concentration, $c(x,t)$, over time \cite{rockne}:
\begin{equation}
\label{td1}
\frac{\partial c}{\partial t} = \frac{\partial}{\partial x}\left(I\frac{\partial c}{\partial x}\right) + pc(x,t) - R(x, t)c(x,t).
\end{equation}
The equation above equates the rate of change of tumor cell concentration with the sum of net dispersal rate, net proliferation rate, and loss of tumor cells due to therapy.
Here, $I$ is the diffusion coefficient representing motility of tumor cells; $p$ is the proliferation rate; and $R$ represents the effect of therapy on cell kill \cite{rockne,harpold}.
Under stand-alone radiotherapy,
\begin{equation}
\label{rrt}
R=1-S.
\end{equation}
Patient-specific invasion and proliferation rates can be obtained using MRI images \cite{rockne}.
We define how supplemental heat affects radiotherapy next.
\vspace{-5pt}

\subsubsection{Radiosensitivity parameter} \vspace{-5pt}

As mentioned in \cref{sec1}, heat affects the radiosensitivity of tissue.
Traditionally, this is expressed using the thermal enhancement ratio (TER) \cite{radbio}, which defines the ratio of dose required to produce a given level of biological damage with and without heat.
Experiments have measured TER for various types of tumors \cite{gillette}.
Some experiments have seen TER values as high as $4.3$ under favorable environments \cite{radbio}.
In this paper, we express this in a different way: we introduce a radiosensitivity parameter, $\xi$, to collectively account for the additive and synergistic effects of hyperthermia.
We define this parameter as the ratio of the biological damage caused by a given amount of radiation dose to tissue with and without heating.
Setting the relevant biological effect to be $R$,
\begin{equation}
\label{ter}
\xi(T_{avg}, t) = \frac{R_{CHR}}{R_{RT}},
\end{equation}
where $R_{CHR}$ quantifies the combined effect of heat and radiation on irradiated tissue, and $R_{RT}$ represents the effect of stand-alone radiotherapy.
$T_{avg}$ is the time-averaged value of temperature at time $t$.
Combining \cref{rrt,ter}, we obtain
\begin{equation}
\label{rchr}
R_{CHR} = \xi R_{RT} = \xi(1-S).
\end{equation}
At normal body temperature, $\xi =1$ irrespective of time; we also set $\xi$ at time zero to be unity.
We arbitrarily assume $\xi=2.5$ for our system when the tumor is heated to an average temperature of $45 \degree$C over $30$ minutes.
In order to define $\xi$ elsewhere, we employ bilinear interpolation \cite{nm} between the data points $\xi(37, 0) = 1$, $\xi(37, 30)=1$, $\xi(45, 0)=1$  and $\xi(45, 30)=2.5$. This returns:
\begin{equation}
\label{rsp}
\xi(T_{avg},t) = 1 - 0.23125T_{avg} + 0.00625(T_{avg}\times t).
\end{equation}
Several simplifying assumptions have been made in this model.
We do not consider temperature dependent cross-sections for transport.
Material density change with heating is ignored.
Moreover, we assume that the decay heat is negligible compared to the constant volumetric heat source in the heat equation.
These simplifications allow us to get a one-way coupled system.
We plan to eliminate these assumptions and model this system with a fully-coupled framework in a later paper.
\vspace{-5pt}

\section{ONE-WAY COUPLED CHR MODELING}
\label{sec3}

The interaction of \cref{transport1,dose1,edose1,survival,td1,heat,rrt,rchr,rsp} determines the tumor response over time.
While our problem is set in slab geometry and is one-way coupled, real-world problems are in general more complex.
Transport alone can depend on seven independent variables \cite{lm}, making its full coupling with other physics difficult.
Moment-based nonlinear acceleration schemes, such as nonlinear diffusion acceleration (NDA), aid in accelerating the convergence of slowly converging physics and allow us to isolate the angular flux from the coupled set of physics by introducing a discretely consistent low-order (LO) equation  \cite{patelFBR}.
The system of equations for nonlinear diffusion acceleration is the following:
\begin{subequations}
\begin{equation}
\label{ho1}
\frac{1}{v}\frac{\partial \psi^{HO}}{\partial t} + \mu\frac{\partial \psi^{HO}}{\partial x} + \sigma_t \psi^{HO} =  \frac{\sigma_{s}(x)}{2} \phi^{LO} + \frac{Q}{2},
\end{equation}
\begin{equation}
\label{lo}
\frac{1}{v}\frac{\partial \phi^{LO}}{\partial t} + \frac{\partial}{\partial x}\left(-D\frac{\partial \phi_{LO}}{\partial x} + \hat{D}\phi^{LO}\right) + \sigma_a \phi^{LO} = Q,
\end{equation}
\begin{equation}
\label{c1}
\hat{D} = \frac{J^{HO} + D\frac{\partial \phi^{HO}}{\partial x}}{\phi^{HO}},
\end{equation}
\begin{equation}
\label{c2}
J^{HO} = \int_{-1}^1 \mu \psi d\mu\,.
\end{equation}
\end{subequations}
Here, $J$ is the current, $D = \frac{1}{3 \sigma_t}$ is the diffusion coefficient, and $\hat{D}$ is the consistency parameter.
The high-order (HO) transport equation (\ref{ho1}) couples with the diffusion-like LO equation (\ref{lo}) via the closure relations described in \cref{c1,c2}, forming the high-order low-order (HOLO) system of equations.
Nonlinear diffusion acceleration is achieved by simultaneous solution of \cref{ho1,lo,c1,c2}.
We discretize the HO equation using the standard backward Euler method in time and diamond-difference/discrete ordinates scheme in space/angle.
The heat, cell population dynamics, and the LO equations employ backward Euler and central finite difference for time/space discretization.
\vspace{-5pt}

Although the schemes used in this paper are well-known, we briefly present them for completeness.
The discretized transport equation is written as:
\begin{subequations}
\begin{equation}
\label{rtdisc}
\frac{1}{v}\frac{\psi_{i,n,k}-\psi_{i,n,k-1}}{\Delta t} + \mu_n \frac{\psi_{i_{1+\frac{1}{2}}, n, k} -\psi_{i_{1-\frac{1}{2}}, n, k}}{\Delta x} + \sigma_{t, i}\psi_{i,n,k} = \frac{\sigma_{s, i}}{2}\phi_{i, k} + \frac{Q_{i,k}}{2}
\end{equation}
\begin{equation}
\label{closure1}
\psi_{i, n, k} = \frac{\psi_{i_{1+\frac{1}{2}}, n, k}+\psi_{i_{1-\frac{1}{2}}, n, k}}{2}.
\end{equation}
\end{subequations}
Here, $i$, $n$, and $k$ represent spatial, angular, and temporal indices respectively.
Integer index values represent cell-centers and fraction values represent cell-edges.
\Cref{dose1,edose1,survival} can be easily discretized in the same fashion.
The LO equation uses the following discretization \cite{knollJFNK}:
\begin{align}
\label{lodisc}
\frac{1}{v} \frac{\phi_{i,k}-\phi_{i,k-1}}{\Delta t} -D \frac{\phi_{i+1, k} -2\phi_{i, k}+ \phi_{i+1, k}}{\Delta x^2} +&\frac{\hat{D}_{i+\frac{1}{2}}}{2 \Delta x}(\phi_{i+1,k} + \phi_{i, k}) \\
& - \frac{\hat{D}_{i-\frac{1}{2}}}{2 \Delta x}(\phi_{i,k} + \phi_{i-1,k}) + \sigma_a \phi_{i,k}  = Q_{i,k}.\nonumber
\end{align}
The standard centered-space-backward-time discretization for the heat equation is given by
\begin{equation}
\label{heatdisc}
\rho c_p \frac{T_{i, k} - T_{i, k-1}}{\Delta t} = \kappa \frac{T_{i+1, k} - 2T_{i, k} + T_{i-1, k}}{\Delta x^2} + q_{i},
\end{equation}
while the cell population dynamics equation is written as
\begin{equation}
\label{tddisc}
\frac{c_{i, k} - c_{i, k-1}}{\Delta t} = I \frac{c_{i+1, k} - 2C_{i, k} + c_{i-1, k}}{\Delta x^2} + p c_{i, k} -R_{i,k}c_{i,k}.
\end{equation}
The coupling and solution can be represented via the following Newton step:
\begin{equation}
\left[\begin{array}{cccccc}
J_{\phi^{LO} \phi^{LO}} & 0 & 0 & 0 & 0 & 0 \\
J_{\phi^{LO}\Phi} & J_{\Phi\Phi} & 0 & 0 & 0 & 0 \\
0 & J_{\Phi D_E} & J_{D_E D_E} & 0 & 0 & 0  \\
0 & 0 & 0 & J_{T T} & 0 & 0 \\
0 & 0 & J_{D_E R} & J_{T R} & J_{R R} & 0\\
0 & 0 & 0 & 0 & J_{R c} & J_{c c}
\end{array}\right] \left[\begin{array}{c}
\delta \phi^{LO} \\
\delta \Phi \\
\delta D_E \\
\delta T \\
\delta R \\
\delta c
\end{array}\right] = -\left[\begin{array}{c}
F_{\phi^{LO}} \\
F_{\Phi} \\
F_{D_E} \\
F_T \\
F_R \\
F_c
\end{array} \right],
\end{equation}
where $F_{\phi^{LO}}$,$F_{\Phi}$, $F_{D_E}$, $F_T$, $F_R$, and $F_c$ are the residual forms of the relevant equations for flux, fluence, effective dose, temperature, effect of therapy on cell kill, and tumor cell concentration.
Since the Jacobian matrix above is not  block-diagonal, the physics are interrelated and therefore warrant a multiphysics treatment.
The one-way coupling (represented by the lower-triangular Jacobian matrix) allows us to solve the problem in a serial fashion, in which we solve transport first to determine flux, fluence, dose, effective dose, and the survival probability (in that order).
Next, we solve the heat equation to determine the temperature, followed by heat-adjusted effect of radiation on tissue.
Finally, we calculate the tumor cell concentration using parameters obtained by solving the previous equations and move on to the next time step.
This kind of modeling allows us to continuously track tumor cell density over time.
For each of the test problems, we choose a uniform time step of $\Delta t = 0.025$ $s$ and $270$ spatial cells.
The solutions obtained are converged to at least four digits precision for this discretization.
We set the angular discretization for the transport equation to eight angles.
The convergence tolerance for transport solves is set to $10^{-6}$.
The heat and tumor dynamics equations are solved using MATLAB's backslash function \cite{matlab}.
\vspace{-5pt}

\subsection{Test Problem}\vspace{-5pt}

We consider the evolution of tumor-cell concentration distribution in a three-region slab over a period of thirty minutes.
The first and the third regions, each $14$ mm thick, represent healthy soft-tissue.
The second soft-tissue region (in the middle) is $2$ mm thick, and has uniformly distributed tumor cells with a concentration of $2^{20}$ $\frac{cells}{cm^3}$.
This region maintains an exponentially decaying internal radiation source with a half life of $75$ days.
The average photon energy of the source is $0.4$ MeV.
In order to induce hyperthermia, this region also has a constant volumetric heat source.
Specific values of the radiation and the heat source will be provided in the next subsection.
We use vacuum boundaries for transport, and open boundaries for the heat and tumor dynamics equations.
We assume the parameters presented in \cref{tab1} \cite{rockne,warrell,kaeri}.
We follow existing literature and assume similar continuum properties for both healthy and tumor tissue.
We represent the transport in healthy and tumor tissues using the same parameters.
The total cross-section (total attenuation coefficient) is $0.1053578$ $cm^{-1}$; and the scattering cross-section is $0.105324$ $cm^{-1}$ \cite{kaeri}.
\begin{table}[!htb]
  \centering
  \caption{\bf{Parameters}}
\begin{tabular}{ |c|c|c|c|c|c|c|c|c|}\hline\label{tab1}
 \textbf{I}  & \textbf{p} & \pmb{$\alpha$}  & \pmb{$\frac{\alpha}{\beta}$}  & $\bm\kappa$ & \pmb{$c_p$} & \pmb{$\rho$}  & \pmb{$\frac{\mu_{en}}{\rho}$} \\
  $\mathrm{[mm/yr]}$ & [/yr]  & [kg/J]&[$\mathrm{J/kg}$] &[W/m.K] &[kJ/kg.K] &[g/cc] & [$\mathrm{cm^2/g}$]  \\
  \hline
  4.29 & 35.13 & 0.203 & 10 & 0.51 & 3.68 & 1.0 & 0.0325\\ \hline
\end{tabular}
\end{table}

\subsection{Numerical Results}\vspace{-5pt}

First, we demonstrate that the use of a scale-bridging HOLO algorithm returns the same tumor response profile to CHR treatment as stand-alone source iteration \cite{lm}. We also analyze how increasing radiation source strength affects the tumor response.
Specifically, we consider three initial source strengths for the tumor region: $Q_0 = 1.4 \times 10^5$, $2.8 \times 10^5$, and $5.6 \times 10^5$ $\frac{photons}{cm^3 s}$.
Everything else is the same as previously described.\vspace{-10pt}
\begin{figure}[bht]
\centering
\begin{subfigure}{.48\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{solver_comparison.jpg}
  \caption{HOLO vs. Source Iteration}
  \label{fig:sub1}
\end{subfigure}
\begin{subfigure}{.48\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{rad_comparison.jpg}
  \caption{Tumor response to different sources}
  \label{fig:sub2}
\end{subfigure}
\caption{Tumor response to stand-alone radiotherapy}
\label{fig:test}
\end{figure}
\vspace{-20pt}

We observe in \cref{fig:test} that both methods return tumor cell concentration within the prescribed tolerance.
This means it is valid to introduce HOLO acceleration into the multiphysics algorithm.
Moreover, according to \cref{fig:sub2}, increasing source strength improves the overall cell kill as expected.
We note that this study is not prescriptive.
\vspace{-5pt}

Next, we consider tumor dynamics in three scenarios: 1) untreated tissue (no radiation or heat source), 2) slab with an initial isotropic internal radiation source of $5.6\times10^{5}$ $\frac{photons}{cm^3 s}$ in the tumor region, and 3) slab with an initial isotropic internal radiation source of $5.6\times10^{5}$ $\frac{photons}{cm^3 s}$ and a heat source of $\frac{1}{3}$ $\frac{W}{cm^3}$ in the tumor region.
\vspace{-5pt}

\Cref{fig:fig1} shows that the introduction of radiation in the tumor begins to kill cancer cells and reduces the total tumor cell concentration in the slab after thirty minutes of treatment.
Introduction of heat enhances the cell kill; this is evident from the further reduction in tumor cell concentration.\vspace{-5pt}

\begin{figure}[h]
\centering
  \includegraphics[width=.58\linewidth]{conc_comparison.jpg}
   \caption{Tumor cell concentration comparison}
  \label{fig:fig1}
\end{figure}

Finally, we compare tumor response to CHR treatment with different heat sources.
We vary the strength of the volumetric heat source in the tumor region of the model slab: $q = \frac{1}{27}$, $\frac{1}{9}$, and $\frac{1}{3}$ $\frac{W}{cm^3}$, and we fix $Q_0$ at $5.6\times10^{5}$ $\frac{photons}{cm^3 s}$.
We find from \cref{fig:fig2} that increasing the strength of the heat source indeed results in a higher tumor cell kill.
This is attributed to increased radiosensitivity in tumor tissue due to the introduction of hyperthermia.
\begin{figure}[!h]
\centering
  \includegraphics[width=.58\linewidth]{heat_comparison.jpg}
   \caption{Tumor response to different heat sources}
  \label{fig:fig2}
\end{figure}

\section{CONCLUSIONS}\label{sec4}

We presented a one-way coupled multiphysics model to determine tumor response to CHR treatment.
We observed that the introduction of radiation results in tumor cell kill.
Moreover, introduction of heat enhances radiosensitivity of the tumor cells, which increases the cell kill.
We made several simplifying assumptions and chose a simple geometry to demonstrate the concept behind this multiphysics model.
We plan to eliminate these simplifying assumptions in subsequent papers such that real-world treatments can be modeled.
Specifically, we will introduce anisotropic scattering, higher dimensions, and energy dependence into our transport model.
This will call for an efficient nonlinear HOLO scheme for forward-peaked transport problems; we will endeavour to develop such a method.
We also plan to extend this model to address random tumor media.
Moreover, we understand the importance of developing a benchmark to verify future models and code implementations.
This will constitute a major portion of our future work in CHR modeling.
\vspace{-5pt}

\section*{ACKNOWLEDGMENTS}

J.~K.~Patel and R.~Vasques acknowledge support under award number NRC-HQ-84-15-G-0024 from the Nuclear Regulatory Commission.
The statements, findings, conclusions, and recommendations are those of the authors and do not necessarily reflect the view of the U.S. Nuclear Regulatory Commission.
\vspace{-5pt}

 \bibliographystyle{mandc}


text:
\section{Introduction}
A Boolean function $f$ is of the form $f:\{0,1\}^{n}\rightarrow\{0,1\}$ (or equivalently $f:\mathbb{V}_{2}^{n}\rightarrow \mathbb{V}_{2}$). The output of the Boolean function $f$ can be represented as a string $s$ of ones and zeros. It can also be represented as a multivariate polynomial over $GF(2)$. This polynomial can be expressed as an exclusive disjunction ($\oplus$) of a constant $a_{0}$ and one or more conjunctions of the function argument. This is called the Exclusive Sum-Of-Product (ESOP) representation. An alternative representation of the ESOP form is known as the Algebraic Normal Form (ANF). The general ANF for a function $f(x_{1},..,x_{n})$ over $n$-variables can be written as,
\begin{equation*}
\begin{split}
f(x_{1},..,x_{n}) = & a_{0} \oplus a_{1}x_{1} \oplus \cdots \oplus a_{i}x_{i} \oplus \cdots \oplus a_{n}x_{n} \\
& \oplus \cdots \oplus a_{1,2,...,n}x_{1}x_{2}\cdots x_{n}\\
\end{split}
\end{equation*}

\subsubsection*{Reversible and Irreversible Boolean functions}
An $n$-variable vectorial Boolean function is termed \textit{reversible} if all its output patterns map uniquely to an input pattern and vice-versa. It can be expressed as an $n$-input, $n$-output bijection or alternatively, as a Boolean permutation function over the truth value set $\{0, 1, \ldots, 2^{n-1}\}$. An \textit{irreversible} Boolean function
$f_{irr}:\{0,1\}^{n}\rightarrow\{0,1\}^{m}$ with $n \ne m$ can also be made reversible with the help of additional output lines, which adds distinguishing patterns to the irreversible output. Correspondingly, additional inputs are added. If an input line is constant-initialized and the constant is recovered after the circuit execution then, it is termed as \emph{ancilla}. Otherwise, it is termed as \emph{garbage}.

\subsection{Reversible Logic Synthesis} \label{sec:rev_synth}
Reversible Boolean logic synthesis is achieved with the help of reversible logic gates. The gates are characterized by their implementation cost in quantum technologies, which is denoted as the Quantum Cost (QC). A reversible gate library is a complete set of reversible gates which can be used to implement any reversible circuit. For example the set of NOT, CNOT, controlled-$V$ and Controlled-$V^+$, known as NCV, is a reversible gate library widely used in the literature. Recently, there has been a significant research activity towards the realization of quantum circuits using Clifford+T gates, considering the importance of fault tolerance in quantum computing. Efficient synthesis for NCV circuits~\cite{mmd_swap}, Clifford+T circuits~\cite{amy_middle} and mapping of NCV circuits to Clifford+T gates~\cite{ncv_to_clifford} have been proposed in the literature. Few gates from these libraries are outlined below. For detailed discussion on primitive quantum gates and their universality, readers may refer to~\cite{elem_gates, barenco_2bit}.

\begin{itemize}[leftmargin=*]
 \item \textbf{NOT gate}: This is represented using the matrix $\bigl(\begin{smallmatrix} 0&1 \\ 1&0 \end{smallmatrix}\bigr)$.
 \item \textbf{CNOT gate}: CNOT$(a,b)$=$(a,a\oplus b)$. This gate can be generalized with $Tof_n$ gate, where first $n-1$ variables are used as control lines. Generalized Toffoli gates form an universal reversible logic gate library, termed as MCT gate library.
 \item \textbf{Hadamard gate}: The unitary transformation for Hadamard gate is, $\frac{1}{\sqrt{2}}\bigl(\begin{smallmatrix} 1&1 \\ 1&-1 \end{smallmatrix}\bigr)$.
 \item \textbf{T gate}: The unitary transformation for this gate is, $\bigl(\begin{smallmatrix} 1&0 \\ 0&e^{\frac{i\pi}{4}} \end{smallmatrix}\bigr)$.
 \item \textbf{Phase gate}: Also denoted by $S$ and performs the unitary transformation of $\bigl(\begin{smallmatrix} 1&0 \\ 0&i \end{smallmatrix}\bigr)$. NOT, CNOT, Hadamard, $S$ and $S^+$ gates together form the Clifford gate family.
 \item \textbf{Toffoli gate}/\textbf{CCNOT gate}: The gate flips the target if both the control lines are set: $Tof(a,b,c) = (a,b,c\oplus ab)$.
\end{itemize}

\subsection{Cost Models} \label{sec:cost}
For evaluating the performance of the synthesis tools and benchmark circuit implementations, different cost models have been proposed in the literature. The most basic cost model is the number of reversible logic gates needed for the implementation. These logic gates could be, with rather large number of control lines, such as, for Multiple-Control Toffoli gates (MCT) or Mixed-Polarity Multiple Control Toffoli (MPMCT) gates. To closely reflect the cost in terms of quantum gates, the QC value is computed for each of these gates, which is nothing but the number of 2-qubit gates~\cite{maslov_web} needed to implement these circuits. In recent fault-tolerant Quantum circuit implementations, the cost is estimated in terms of T-count, corresponding to the Cliffort+T realization. Logical depth or T-depth.

Finally, considering the cost and difficulty of implementing large number of qubits, one performance indicator is the total number of lines or ancilla/garbage count. It is well-known that there exists a trade-off between gate count/QC and ancilla count~\cite{integration_tradeoff,anupam_tradeoff}. It is trivial to show that with increasing ancilla count, MCT gates with smaller control lines can be realized, thereby establishing a clear trade-off between QC and ancilla. However, it not clear if such a trade-off exists for gate count vs ancilla.

\section{Background} \label{sec:existing_bounds}
For theoretical as well as practical interests, several researchers worked towards establishing the bounds on the gate count for reversible circuits. For selective classes of Boolean functions, much lower upper bound compared to generic Boolean functions has been reported in~\cite{constructive}. For reversible circuits of practical interest, e.g., modulo-exponentiation, addition, reversible circuit implementation with parameterized upper bounds have been established~\cite{takahashi_adder,markov_survey}. We outline the major results in the gate count upper bound for different synthesis methods.

\subsection{Upper Bound for Transformation-based Synthesis Methods} \label{sec:mmd_bound}
The gate cost upper bound for transformation-based synthesis methods, for $n$-variable reversible logic circuits, have been reported as following for two different gate libraries.
\begin{lemma}
The upper bound of gate count for MMD~\cite{mmd} using MPMCT library is $(n-1)2^n + 1$~\cite{mmd_tcad,mmd_swap}.
\end{lemma}

\begin{lemma}
The upper bound of gate count for MMD~\cite{mmd} using MPMCT(+Fredkin gate) library is $(n-2)2^n + 2 + n$~\cite{mmd_swap}.
\end{lemma}

\subsection{Upper Bound for BDD-based Synthesis Methods} \label{sec:bdd_bound}
The Binary Decision Diagram (BDD)-based synthesis method, originally proposed in~\cite{wille_bdd}, observed an gate count upper bound of $3.2^n$. Note that this method generates high number of ancilla, which is in the order of $O(2^n)$. By utilizing a bi-conditional decomposition, for selective Boolean functions, a tighter upper bound is proposed in~\cite{bbdd}, however, the worst-case upper bound remains the same for generic Boolean functions. Ancilla-free reversible logic synthesis for BDDs has been proposed recently~\cite{mathias_ancilla_free_bdd}, which does not report any theoretical bounds so far.

\subsection{Upper Bound for ESOP-based Synthesis Methods} \label{sec:esop_bound}
The tightest upper bound for reversible circuits have been reported in~\cite{nabilla_bound} following a mix of Young subgroups and Exclusive Sum-Of-Products (ESOP)-based synthesis approach. In the first phase of the synthesis with young subgroups, it is established that an $n$-variable reversible function can be implemented with a circuit of $(2n-1)$ single target gates. Further, it was shown that, an $n$-variable single target gate can be realized with at most $29\cdot2^{n-8}$ MPMCT gates, if $n \geq 8$. Therefore, the gate count upper bound for the reversible circuit is $29\cdot2^{n-8}\cdot(2n-1)$ MPMCT gates, if $n \geq 8$. The upper bound of the single target gate implementation is established by using a result from~\cite{gaidukov}. Note that this method does not require any ancilla line.

For the rest of the paper, we will closely follow the above synthesis method and focus on the complexity of a single-output Boolean function. We will assume the single target distribution obtained from Young subgroup method~\cite{nabilla_bound} and concentrate on the functional decomposition approach of Boyar et al~\cite{boyar_bool} to check for potential improvements.

\subsection{Multiplicative Complexity} \label{sec:mult_complexity}
In parallel with the efforts for enumerating the reversible circuit complexity, progress in the analysis of classical Boolean functions in terms of multiplicative complexity is reported~\cite{boyar_bool}. The primary application for this analysis has been towards the design of cryptographic primitives, though.

\begin{definition}
The multiplicative complexity $C_\wedge(f)$ of a Boolean function is the minimum number of multiplications (AND-gates) that are sufficient to evaluate the function over the basis $(\wedge, \oplus, \neg)$.
\end{definition}

It has been established that the multiplicative complexity of functions having degree $d$ is at least $d-1$~\cite{schnorr_mult_c} and the multiplicative complexity of a randomly selected $n$-variable Boolean function is at most $2^{\frac{n}{2} + 1} - n/2 -2$ ($n$ even)~\cite{boyar_bool}.

Intuitively, multiplicative complexity of a Boolean function has some correspondence with the Toffoli gate. The number of AND gates in a single term of ESOP expression represents the number of control lines in a Toffoli gate. In the following sections, we follow the same approach of the derivation of multiplicative complexity and explore if it reduces the theoretical upper bound of gate count compared to the currently known bounds.

\section{Upper Bounds via Multiplicative Complexity Analysis} \label{sec:tighter_bound}
In this section, we first revisit the function decomposition approach used in~\cite[Section 4]{nabilla_bound} from a mathematical perspective. We show that their work can be described as solving a linear recurrence of the form: $$f(n) = 2f(n-1) + 2$$ with the initial conditions $f(5) = 10$, $f(4) = 4$.

The reason for adopting the recurrence relation, $f(n) = 2f(n-1) + 2$ is described in second paragraph of the proof of \cite[Theorem 1]{nabilla_bound}, so we skip that part. The case for $f(5)$ is described in~\cite{wille_sat}, which is cited in~\cite{nabilla_bound}.

Now, assume that, we want to solve: $f(n) = 2f(n-1) + 2$.
We substitute $n$ by $n+1$ to yield, $f(n+1) = 2f(n) + 2$.
Now, subtracting, we get $f(n+1) - f(n) = 2f(n) -2f(n-1)$ $\implies f(n+1) = 3f(n) - 2f(n-1)$.

Hence, the characteristic equation is: $x^2=3x-2 \implies (x-2)(x-1) =0 \implies x=1,2$, and hence, $f(n) = \alpha\cdot1^n+\beta\cdot2^n$, where $\alpha, \beta$ are constants to be determined from the initial conditions.

Now, given that, $f(5)=10, f(4)=4$,
we generate the linear equations: $f(5)= \alpha + 32\beta = 10, f(4)=\alpha+16\beta=4$, which gives, $\alpha=-2, \beta=\frac38=3\cdot2^{-3}$.

Therefore, $f(n) = -2+3\cdot2^{-3}\cdot2^n = 3\cdot2^{n-3}-2$.

\subsection{Based on Function Decomposition}
A recursive decomposition procedure can be used for systematic implementation of reversible Boolean functions. From a given ESOP form, three following functional decompositions can be applied.
\begin{equation}
f = \overline{x_i}\cdot f_{x_i=0} \oplus x_i\cdot f_{x_i=1} \label{eqn-shannon}
\end{equation}
\begin{equation}
f = f_{x_i=0} \oplus x_i\cdot f_{x_i=2} \label{eqn-pdavio}
\end{equation}
\begin{equation}
f = f_{x_i=1} \oplus \overline{x_i}\cdot f_{x_i=2} \label{eqn-ndavio}
\end{equation}
where $f_{x_i=2} = f_{x_i=0} \oplus f_{x_i=1}$. These are known as Shannon, positive Davio and negative Davio decompositions respectively. Here, we adopt the positive Davio decomposition. The Shannon decomposition leads to a $2$-MCT gate implementation, whereas both positive Davio and negative Davio can be realized with a single MCT or MPMCT gate respectively. Given a Boolean function $f \mt{n}$ in its ESOP representation, we implement ($\wedge$ and $\oplus$ gates) it with MCT gate library only. We denote the MCT gate complexity of $f$ by \tc{f}.

In~\cite[Lemma 5]{boyar_bool}, it was proved that it requires at most $2^n - n - 1$ AND ($\wedge$) gates for realizing all the positive minterms of $n$ Boolean variables. In that, authors proposed a circuit of alternating $\wedge$ and $\oplus$ gates. It is evident that this set of all possible minterms can also be implemented with only MCT gates, when $0$-initialized ancilla lines are used. A representative implementation is shown in the following Fig.~\ref{minterm}. This requires $2^n - n - 1$ MCT gates and total $2^n-1$ lines.

\begin{figure}[htb]
\centering
$
\Qcircuit @C=0.7em @R=.5em {
\lstick{x_1} & \ctrl{1} & \ctrl{2}   & \qw      & \ctrl{1} & \qw & \rstick{x_1}\\
\lstick{x_2} & \ctrl{2} & \qw        & \ctrl{1} & \ctrl{1} & \qw & \rstick{x_2} \\
\lstick{x_3} & \qw      & \ctrl{2}   & \ctrl{3} & \ctrl{4} & \qw & \rstick{x_3} \\
\lstick{0}   & \targ    & \qw        & \qw      & \qw      & \qw & \rstick{f = x_1x_2}\\
\lstick{0}   & \qw      & \targ      & \qw      & \qw      & \qw & \rstick{f = x_1x_3}\\
\lstick{0}   & \qw      & \qw        & \targ    & \qw      & \qw & \rstick{f = x_2x_3}\\
\lstick{0}   & \qw      & \qw        & \qw      & \targ    & \qw & \rstick{f = x_1x_2x_3}
}
$
\caption{Implementation of Minterms}\label{minterm}
\end{figure}

\begin{lemma} \label{lemma_xor}
If all the positive minterms of $n$ variables are computed, then, for realizing $r$ different $n$-variable Boolean functions over the basis $(\wedge, \oplus, 1)$, at most $r(\frac{2^n-1}{2} - 1)$ MCT gates are needed.
\end{lemma}

\begin{proof}
Every function $\{f_i | 1 \leq i \leq r\}$ from the above set of $r$ functions can use a set of positive minterms to prevent formation of sharing co-factors. To begin with, we assume that the MCT of all minterms (except constant function) is computed, which requires $2^n-2$ MCT gates. Further, the minterm $0$ is made available as an ancilla line. Therefore, for computing each function, we need to look into the complete set and the zero set. In either case, the target function will never require more than half of the minterms to undergo exclusive-OR operation. However, if we assume that $f_1$ to need half of the minterms and $f_2$ to need the other half, there is a possibility that $f_3$ will require a few more minterms from the $f_1$ and a few minterms from $f_2$. In that way, storing the intermediate functions will not be useful and we get an upper bound of $(\frac{2^n-1}{2} - 1)$ MCT gates for each of the functions.
\end{proof}

We start with an idea similar to~\cite[Theorem 6]{boyar_bool}. Following positive Davio decomposition, $\forall f \mt{n}$, it is possible to find $f_1,f_{0} \mt{n-1}$ such that, $$ f(x_1,x_2,\ldots,x_n) = x_1f_1(x_2,\ldots,x_n) \oplus f_{0}(x_2,\ldots,x_n).$$ Alternatively, it can be stated that, $f$ is divided by $x_1$, with $f_1$ being the quotient and $f_{0}$ being the remainder. For example,  if $f(x_1,x_2,x_3) = x_1x_2x_3 \oplus x_1x_2 \oplus x_2x_3 \oplus x_1 \oplus x_2\oplus 1$, then $f_1(x_2,x_3)=x_2x_3\oplus x_2\oplus1$ and $f_{0}(x_2,x_3) = x_2x_3\oplus x_2\oplus1$. Also, notice that both $f_0, f_1$ are of (at most) $n-1$ variables.

Note that, if we can implement $f_1$ and $f_0$, then we can readily implement $f$ just by using one MCT gate (using $x_1, f_1$ as the control lines and $f_0$ as the target line, as depicted in Fig.~\ref{stp-1}).

\begin{figure}[htb]
\centering
$
\Qcircuit @C=0.7em @R=.5em {
\lstick{x_1} & \ctrl{1} & \qw & \rstick{x_1}\\
\lstick{f_1} & \ctrl{1} & \qw & \rstick{f_1} \\
\lstick{f_0} & \targ & \qw & \rstick{f = x_1f_1 \oplus f_0}
}
$
\caption{Implementation of $f$ from $f_1$ and $f_0$ as inputs}\label{stp-1}
\end{figure}

Thus, we make the following observation: $$\tc f \leq 1 + \tc{\{f_1, f_0\}},$$ where \tc{\{f_1,f_0\}} denotes the MCT gate complexity of $f_1,f_0$, combined. Therefore, we can express the MCT gate complexity of an $n$-input function $f$ in terms functions of lower number of input variables. We can recursively use this decomposition:
\begin{align*}
{f} = &{x_1f_1\oplus f_{0}} \\
    = &{x_1(x_2f_{11}\oplus f_{10})\oplus (x_2f_{01}\oplus f_{00})} \\
    = &{x_1(x_2(x_3f_{111}\oplus f_{110})\oplus x_3f_{101}\oplus f_{100})} \\
      & \oplus {(x_2(x_3f_{011}\oplus f_{010})\oplus x_3f_{001}\oplus f_{000})}
\end{align*}
It can be noted that if we conduct this decomposition in steps of $4$ factors, instead of steps of $2$, exactly the same number of gates would be needed at the end. This is due to the fact that the $4$-factor decomposition, if implemented in reversible circuit, would need $3$ MCT gates.

In the above decomposition, at $k$-th step, we require $2^k-1$ additional MCT gates other than the ones needed by the new functions of lower input variables. Thus, after $k$-th step, we can write:
\begin{align}
&\tc{f} \notag\\& \leq 1 + \tc{\{f_1,f_0\}} \notag\\ & \leq 3 + \tc{\{f_{11}, f_{10}, f_{01},f_{00} \}} \notag
\\ & \hspace{.29\textwidth}\vdots  \notag\\ & \leq (2^k-1) + \tc{\{f_{i_1,i_2,\cdots,i_k}|{i_1,i_2,\ldots,i_k}\in \{0,1\}^k\}} \label{eqn-1}
\end{align}

Note that, all functions $\{f_{i_1,i_2,\cdots,i_k}|$ ${i_i,i_2,\ldots,i_k}\in$ $\{0,1\}^k\}$ are (at most) of $n-k$ variables, hence, $\{f_{i_1,i_2,\cdots,i_k}|$ ${i_1,i_2,\ldots,i_k}\in$ $\{0,1\}^k\}$ denotes some circuit with $n-k$ variables.

Therefore, using the lemma~\ref{lemma_xor} (here, $r = 2^k$) and the MCT gate count for the minterm constructions, in equation~\ref{eqn-1}, we can write
\begin{align}
& \tc{f} \notag\\
& \leq (2^k-1) + \tc{\{f_{i_1,i_2,\cdots,i_k}|{i_1,i_2,\ldots,i_k}\in \{0,1\}^k\}} \notag\\
& \leq (2^k-1) + \{2^{n-k} - (n-k) - 1\} + 2^k\{\frac{2^{n-k}-1}{2} - 1\} \notag\\
& \leq (2^k-1) + \{2^{n-k} - (n-k) - 1\} + 2^{n-1} - 3\cdot2^{k-1}
\label{eqn-AC}
\end{align}

Assuming $k=n/2$ to be the point, where we stop the functional decomposition and construct the functions from the available positive minterms, we obtain
\begin{align}
\tc{f} \leq 2^{n-1} + 2^{n/2-1} - n/2 - 2
\label{eqn-finalonegate}
\end{align}

This leads to an overall MCT gate complexity of
\begin{align}
(2n-1)(2^{n-1} + 2^{n/2-1} - n/2 - 2)
\label{eqn-finalMCT}
\end{align}
by following the Young subgroup decomposition. Corresponding bound for the garbage count is $O(2^{n/2})$, which is introduced for the minterm construction and for the realization of the functions once the decomposition procedure is stopped.

\begin{figure}[htbp]
\resizebox{0.9\columnwidth}{!}{
\centering
\includegraphics[width=\columnwidth]{graph.pdf}
}
\caption{Upper Bound of MCT Gate Complexity}\label{fig:complexity}
\vspace{-0.3cm}
\end{figure}

\subsection{Discussion}
It can be noted that the MCT gate complexity, derived through the approach of~\cite{boyar_bool} is worse compared to that of the transformation-based method (for $n<8$) and the ESOP-based method (for $n\geq 8$) (see Fig. \ref{fig:complexity}). This leads to the conclusion that introducing ancilla/garbage lines does not lead to any reduction in the gate count complexity, if functional decomposition approach is followed. On the other hand, different complexity bounds in different synthesis approaches indicate that it might be beneficial to adopt a hybrid synthesis approach for reversible circuit realization.

\bibliographystyle{IEEEtran}


text:
\section{Introduction}

The LHCb collaboration \cite{LHCbB} claims two candidates for pentaquark states
which are seen in the decay $\Lambda^0_b\to K^-J/\psi\, p$.
These are narrow peak $5/2^?(4450\pm 4)$ with a width $\Gamma=(39\pm 24)$ MeV
and a broad state  $3/2^?(4380\pm 38)$ with
$\Gamma=(205\pm 94)$ MeV.
In the $(K^-\,p)$-channel a set of the observed $\Lambda$-states \cite{PDG} is taken into
account in the fit procedure of ref. \cite{LHCbB} while an existence of resonances in the
$(K^-J/\psi)$-spectrum is turned down. But
in the quark-diquark schemes the adoption of the pentaquark states means an
adoption of the tetraquark states as well. Then the $(K^-J/\psi)$-spectrum should be
filled up by exotic mesons, supposing these exotic states exist.

In this note we emphasize that the strange exotic mesons with hidden charm,
$(c\bar c\, s\bar u)$, (masses and spins have been estimated in \cite{ala-1}) can decay
by process $(c\bar c\, s\bar u)\to K^-J/\psi$ thus imitating contributions from the
$pJ/\psi$-channel.

\begin{figure}
\centerline{\epsfig{file=mjpsik_LHCb.eps,height=0.50\textwidth}}
\caption{\label{f1}
Fig.11 of the ref. [1]: Projections onto $m_{J/\psi K}$ in various intervals of
$m_{Kp}$ for the reduced model fit (cFit) with two $P^+_c$ states of $J^P$ equal to
$3/2^-$ and $5/2^+$: (a) $m_{Kp}< 1.55$ GeV, (b) $1.55<m_{Kp}< 1.70$ GeV,
(c) $1.70<m_{Kp}< 2.00$ GeV, (d) $2.00<m_{Kp}< 2.00$ GeV. The data are shown as (black)
squares with error bars, while the (red)
circles show the result of the
fit. The individual resonances are given in the legend
. }
\end{figure}

\section{Exotic meson states}

Strange exotic mesons with hidden charm were discussed in \cite{ala-1} as composite diquark-antidiquark systems.
In terms of the quark model we have found the masses and spins of such states.
Results for the $(cs\cdot\bar c\bar s)$ and $(cq\cdot\bar c\bar s)$ systems read:
\be \label{33}
\begin{tabular}{l|l|l||l|l|l}
&$(cs)\cdot(\bar c\bar s)$ & mass  & & $(cq)\cdot(\bar c\bar s)$& mass  \\
$J^{PC}$&$I=0$&MeV  & $J^{P}$& $ I=1/2$& MeV   \\
\hline
$2^{++}$& $ A_{(cs)}\cdot A_{(\bar c\bar s)}$&$4630\pm 50$ &
$2^{+}$& $ A_{(cq)}\cdot A_{(\bar c\bar s)}$& $4530\pm 50$\\
$1^{+-}$& $ A_{(cs)}\cdot A_{(\bar c\bar s)}$&$4400\pm 50$ &
$1^{+}$& $ A_{(cq)}\cdot A_{(\bar c\bar s)}$& $4300\pm 50$\\
$0^{++}$& $ A_{(cs)}\cdot A_{(\bar c\bar s)}$&$4280\pm 50$ &
$0^{+}$& $ A_{(cq)}\cdot A_{(\bar c\bar s)}$& $4180\pm 50$\\
$1^{++}$& $ [A_{(cs)}\cdot S_{(\bar c\bar s)}]_{sym}$&$4380\pm 50$ &
$1^{+}$& $ A_{(cq)}\cdot S_{(\bar c\bar s)}$& $4280\pm 50$\\
$1^{+-}$& $ \{S_{(cs)}\cdot A_{(\bar c\bar s)}\}_{asym}$&$4330\pm 50$ &
$1^{+}$& $ S_{(cq)}\cdot A_{(\bar c\bar s)}$& $4230\pm 50$\\
$0^{++}$& $ S_{(cs)}\cdot S_{(\bar c\bar s)}$&$4140\pm 50$ &
$0^{+}$& $ S_{(cq)}\cdot S_{(\bar c\bar s)}$& $4040\pm 50$\\
 \end{tabular}
\ee
Here $A_{(cs)}$ and $ S_{(\bar c\bar s)}$ are axial ($1^+$) diquark and scalar ($0^+$)
antidiquark correspondingly. The masses of the states $(A_{(cs)}\cdot A_{(\bar c\bar s)})$,
$(A_{(cs)}\cdot S_{(\bar c\bar s)})$, $(S_{(cs)}\cdot S_{(\bar c\bar s)})$ are
determined basing on the Belle \cite{belle-ss}, CDF \cite{cdf-ss}, CMS \cite{CMS-ss},
D0 \cite{D0-ss} data by the diquark - antidiquark model \cite{ala-1}. This model
treats the exotic meson states
 as two-component composite systems with\\
(i) diquark-antidiquark component $(cs)\cdot(\bar c\bar s)$, \\
(ii) meson-meson component $(c\bar s)\cdot(s\bar c)$. \\
The masses of exotic mesons with open strangeness,
$(cq)\cdot(\bar c\bar s)$ and  $(cs)\cdot(\bar c\bar q)$, are fixed by
strange quark weighting, it is accepted to be
$\Delta m_s =100$ MeV (this point is discussed in \cite{ala-1}).

\section{Interplay of of exotic meson and baryon states}

The $(K^-J/\psi)$-spectra  measured in \cite{LHCbB}
 are shown on Fig. \ref{f1} (figure 11 of ref. \cite{LHCbB}),
the mass distributions have structures which nicely coincide with masses of the
$(cq)\cdot(\bar c\bar s)$-states shown in the last column of eq. (\ref{33}).
Let us emphasize that Fig. \ref{f1}d demonstrates data events from the mass region
with a minimal contribution of the $\Lambda$-resonances thus providing argument
for search
of resonances in the $(K^-J/\psi)$ channel.

We observe a crossing of contributions of the baryon state  $3/2^?(4380\pm 38)$ and
meson states with $J^P=0^+,1^+$ in the mass interval
$M_{(cq)\cdot(\bar c\bar s)}\sim (4180-4300)$ MeV,
and that require a study of their interplay.
The meson states can imitate a signal of the broad baryon resonance at mass near 4380 MeV.

\section{Conclusion}

A wide discussion accompanies the observation of structures in the
$p\,J/\psi $ channel
\cite{ala-2,maia-polo,mikh,liu-wang,meis-olle,miro,he,lebe,srr,x-meis,kuba-wolo,wang,he2,burns}.
In line with
\cite{ala-2} we discuss a version which suggests that
the narrow peak
 $5/2^?(4450\pm 4)$ is the genuine pentaquark state with $J^P=5/2^-$
 and the broad bump  $3/2^?(4380\pm 38)$
is mainly resulted due to exotic mesons $1^+$ in the $(K^-J/\psi)$ channel
with masses $\sim 4189-4300$ MeV.
We conclude to study the broad pentaquark states one needs to determine
the character of the irregularities in the $(K^-J/\psi)$ channel.

We are grateful to Y.I. Azimov and J. Nyiri for helpful discussions.
 The work was supported by grants RSGSS-4801.2012.2,
 RFBR-13-02-00425 and RSCF-14-22-00281 .