\section{Introduction and discussion of the results}
The heat semigroup  $e^{t\Delta}$  associated to standard Laplacian has been studied by many authors in PDEs and physics. In fact, the study of the heat semigroup pervades throughout mathematical analysis and physics, being indispensable in many situations. Moreover, there has been increasing interest in applications of the fractional Laplacian to the mathematical modelling of various physical phenomena, see e.g.\ \cite{ng, NL} and the references therein. The Hermite operator (also known as quantum harmonic oscillator) $H=-\Delta+|x|^2$ plays a vital role in quantum mechanics and analysis \cite{folland, st}. Nevertheless, there are only few mathematical papers which deal with fractional powers of Hermite operator $H^{\beta} \ (\beta>0),$ see e.g.\ \cite{bhimani1,EC,st1}.

 The  spectral decomposition  of $H$ on $\mathbb R^d$ is given by
\begin{equation}\label{eq spec dec}
H=\sum_{k=0}^{+\infty}(2k+d)P_k,\qquad P_k f=\sum_{|\alpha|=k}\langle f,\Phi_\alpha\rangle\Phi_\alpha,
\end{equation}
where $\langle\cdot,\cdot\rangle$ is the inner product in $L^2(\bR^d)$ and $\Phi_\alpha$, $\alpha\in \mathbb{N}^{d}$, are the normalised Hermite functions, forming an orthonormal basis for  $L^2(\bR^d)$. Observe the regularising effect of $P_k$, which takes temperate distributions into Schwartz functions: $P_k:\mathcal{S}'(\bR^d)\to\mathcal{S}(\bR^d)$. Since 0 is not in the spectrum of $ H $ we can define the fractional powers $ H^\beta$ for any $ \beta \in \mathbb{R} $ by means of the spectral theorem, namely
$$  H^\beta f = \sum_{k=0}^\infty (2k+d)^\beta P_kf.$$
We remark  that $ H^\beta $ is a densely defined unbounded operator for $ \beta > 0 $. We thus define the heat semigroup associated to $H^{\beta}$  $(\beta>0)$ by
\begin{equation}\label{eq heat semi}
e^{-tH^\beta}f=\sum_{k=0}^{+\infty} e^{-t(2k+d)^\beta} P_kf.
\end{equation}
In this note we study the behaviour of this semigroup on modulation spaces. In order to define these spaces, we recall the   definition  of short-time Fourier transform (STFT - also known as the Bargmann transform \cite{tataru}) of $f\in \mathcal{S}'(\bR^d)$ with respect to a fixed window function $ g\in \mathcal{S}(\mathbb R^d)\setminus\{0\}$:
\[
V_g f(x,\xi) \coloneqq \int_{\bR^d} e^{-i \xi \cdot y}f(y) \, \overline{g(y-x)} \, dy,\quad x,\xi\in\bR^d.
\]
We then measure the phase-space content of $f$ by means of mixed Lebesgue (quasi-)norms $L^{p,q}(\bR^d \times \bR^d)$ of $V_g f$, leading to the so-called modulation spaces $ M^{p,q}$ \cite{F1,gro book,wang}:
\[
\|f\|_{M^{p,q}}:=\| \|V_g f(x,\xi)\|_{L^p_x}\|_{L^q_\xi},
\]
with $0<p,q\leq\infty$; see Section \ref{PR} for a more general definition involving weights. Heuristically, one can think of a function in $M^{p,q}$ as having the local regularity of a function whose Fourier transform is in $L^q$ and decaying at infinity as a function in $L^p$. We have, in particular, $M^{2,2}=L^2$. Modulation spaces can be equivalently designed as a family of Besov-type spaces with the dyadic geometry in frequency replaced by a decomposition in isometric boxes.  \par
We now state our main result, concerning the action of $e^{-tH^\beta}$ on such spaces.
\begin{theorem}\label{mainthm} Let $\beta>0$, $0< p_1,p_2,q_1,q_2\leq\infty$ and
\[
 \frac{1}{\tilde{p}}\coloneqq \max\Big\{\frac{1}{p_2}-\frac{1}{p_1},0\Big\},\ \ \frac{1}{\tilde{q}}\coloneqq \max\Big\{\frac{1}{q_2}-\frac{1}{q_1},0\Big\},\quad \sigma\coloneqq  \frac{d}{2\beta} \Big(\frac{1}{\tilde{p}}+\frac{1}{\tilde{q}}\Big).
 \] Then
\begin{equation}\label{eq mainthm}
\|e^{-tH^\beta} f\|_{M^{p_2,q_2}}\leq C(t)\|f\|_{M^{p_1,q_1}}
\end{equation}
for every $t> 0$, where
\begin{equation}\label{optimalconstant}
C(t)=\begin{cases}
C_0e^{-td^\beta} & t\geq 1\\
C_0t^{-\sigma}& 0<t\leq 1
\end{cases}
\end{equation}
for some $C_0>0$.
\end{theorem}
To the best of our knowledge the result is new even in the case $\beta=1$, $p_1=q_1=p_2=q_2$. \par
Let us give a flavour of the heuristics behind the behaviour of $ C(t).$ First, by testing the above estimate on the ground state of $H$ we see that the decay at infinity is absolutely sharp for every choice of the exponents and also that the exponent of $t$ for $t$ small can never be positive. Indeed, if  $ f = \Phi_0$ then $P_0f=f$ and  $P_k f=0$ for $k\not=0$, so that  $e^{-tH^\beta} f=e^{-td^\beta} f$.

As the modulation spaces $ M^{p,q} $ increase when one of the exponents increases while the other one is kept fixed, it is sufficient to prove the result with $p_2$ replaced by $\min\{p_1,p_2\}$ and similarly for $q_2$, namely we can assume $p_2\leq p_1$, $q_2\leq q_1$. Now, in somewhat sloppy terms, the effect of the map $e^{-t H^\beta}$ in phase space is to damp the content of a function in a way that roughly amounts to the multiplication by the function
\[
F_t(x,\xi)=e^{-t(|x|^2+|\xi|^2)^\beta}.
\] The operator norm of this pointwise multiplication operator, as a map $L^{p_1,q_1}\to L^{p_2,q_2}$ (mixed-norm Lebesgue spaces in $\bR^d\times \bR^d$), can be computed by H\"older's inequality, namely $\|F_t\|_{L^{\tilde{p},\tilde{q}}}=C(t)$, as given in the statement ($0<t\leq1$). Of course the negative exponent of $t$, for $t$ small, is the price to pay for the {\it phase-space smoothing}, when passing from $p_1$ to $p_2<p_1$ or from $q_1$ to $q_2<q_1$. \par
We emphasize that the exponents are just assumed to be positive, without further conditions, similarly to the refined estimates for the heat semigroup in Besov spaces \cite{wang240} (see also \cite[Sec.\ 2.2]{wang}) and in real Hardy spaces \cite{chen}.  \par
We also observe that the same estimates hold  when $M^{p,q}$ is  replaced by the so-called Wiener amalgam space $W^{p,q}\coloneqq \mathcal{F} M^{p,q}$, i.e.\ the image of $M^{p,q}$ under the Fourier transform, endowed with the obvious norm. This follows at once by applying the above estimates to $\mathcal{F}^{-1} f$ and using the fact that $\mathcal{F}^{-1}$ commutes with the spectral projections $P_k$ (because $\mathcal{F}^{-1}\Phi_\alpha=i^{|\alpha|}\Phi_\alpha$) and therefore with $e^{-t H^\beta}$ by \eqref{eq heat semi}.

Similar estimates for the fractional heat semigroup $e^{-t(-\Delta)^\beta}$ were obtained in \cite[Thm.\  3.1]{chen}; see also \cite{wang} for related results. However, our analysis follows a completely different pattern and is necessarily less elementary, because the operator $e^{-tH^\beta}$ is not a Fourier multiplier. Actually, except for the case $\beta=1$, there is not even an explicit integral formula for $e^{-tH^\beta}$ and we rely on the theory of fractional powers and heat kernels of globally elliptic pseudodifferential operators in $\bR^d$ (see e.g.\ \cite{NR,shubin} for a general account).\par
From the above fixed-time estimates one could deduce Strichartz (space-time) estimates by a standard machinery, via the $TT^\ast$ method or real interpolation (see e.g.\ \cite{tao} and \cite[Thm.\  3.1]{chen} respectively; see also \cite{pierf}). However, in the case of the Hermite semigroup one should be able to obtain a broader range of space-time estimates beyond those derived by the fixed-time estimates, see for instance \cite[Cor.\ 2.1]{wang} for the heat semigroup. Hence we prefer to postpone a systematic study of Strichartz estimates, including some interesting related topics - in fact, it seems that the techniques of this paper could be successfully applied to obtain some new $L^p$ estimates as well. Instead, here we focus on some direct applications of Theorem \ref{mainthm} to the heat equation for $H^\beta$ with a nonlinearity of power type, providing some results which are definitely a consequence of the trapping effect of the quadratic potential in $H$ and do not hold for the corresponding heat equation without potential.

Specifically, we  consider the Cauchy problem for the nonlinear heat equation\footnote{The subsequent arguments and results can be trivially modified for other algebraic nonlinearities such as $\lambda u^k$,
	$\lambda\in\mathbb{C}$, $k\in\mathbb{N}$, $k\geq 2$.} associated to $H^{\beta}:$
\begin{equation} \label{HE}
\begin{cases} \partial_t u+H^\beta u=\lambda |u|^{2k} u\\
u(0,x)=u_0(x)
\end{cases}
\end{equation}
with $(t,x)\in(0,+\infty)\times\bR^d$, where $k\in \mathbb{N}$, $k\not=0$, $\lambda\in\mathbb{C}$ and $\beta>0$. \par
As an application of Theorem \ref{mainthm}, we   establish global well-posedness for \eqref{HE}. Specifically, we  have following result.
\begin{theorem}\label{mt2}
Let $ 1\leq p \leq \infty$, $1\leq q \leq \frac{2k+1}{2k}$ and
\[
\frac{1}{q}+\frac{\beta}{kd}>1.
\]
  Define the subspace $X\subset L^\infty([0,+\infty),M^{p,q})$ of elements $u$ satisfying
\[\|u\|_{X}\coloneqq  \left\| e^{td^{\beta}}  \|u(t, \cdot)\|_{M^{p,q}} \right\|_{L_t^{\infty}([0, \infty))}< \infty .\]
\begin{itemize}
\item[(a)] There exists $\varepsilon>0$ such that for any $u_0\in M^{p,q}$  satisfying  $\|u_0\|_{M^{p,q}} \leq  \varepsilon,$  the problem \eqref{HE} has
	 a unique global solution
	 \[u \in L^{\infty}([0, \infty), M^{p,q}). \]
\item[(b)] We have $u \in C([0,\infty), M^{p,q})$  if $p<\infty$. \par\noindent
\item[(c)] If $\varepsilon$ is small enough the above solution enjoys exponential decay in time; specifically,   $u \in X.$	 \par\noindent
\item[(d)] The same results hold with $M^{p,q}$ replaced by $W^{q,p}$, for the same range of exponents.
\end{itemize}
\end{theorem}
The following remarks are in order (see Example \ref{ex} for details).

\begin{remark}\label{fe}
Let $f_\alpha(x)=|x|^{-\alpha}$, $0\leq\alpha<d$, $x\in\bR^d$.
\begin{itemize}
\item[(i)]
The hypothesis  $1\leq q\leq \frac{2k+1}{2k}$ in Theorem \ref{mt2} is natural.  Indeed, $f_\alpha\in M^{\infty,q}\subset W^{q,\infty}$, $q=\frac{2k+1}{2k}$, if $0\leq \alpha<\frac{d}{2k+1}$. But in the limiting case $\alpha=\frac{d}{2k+1}$ we
have $|f_\alpha(x)|^{2k} f_\alpha(x)=|x|^{-d}$, which is not even locally integrable. Notice that when $\beta>\frac{kd}{2k+1}$, under the assumption $1\leq q\leq\frac{2k+1}{2k}$ the additional condition $\frac{1}{q}+\frac{\beta}{kd}>1$ is automatically satisfied.
\item[(ii)] To give a flavour of the the type of singularities and oscillations at infinity admitted for the initial data, we observe that, for example, one can take as initial datum $u_0(x)=\varepsilon f_\alpha(x)(1+c\cos |\xi|^2)$, $c\in\mathbb{C}$, $0\leq\alpha<\min\{d/(2k+1),\beta/k\}$, with $\varepsilon$ small enough. Also, if $\chi$ is any smooth function with compact support in $\bR^d$ and $\Lambda$ is any lattice in $\bR^d$, one can consider $u_0(x)=\varepsilon\sum_{\mu \in\Lambda} f_\alpha(x-\mu)\chi(x-\mu)$, with the same restrictions on $\alpha$ and $\varepsilon$. Observe that of course $f_\alpha\not\in L^p$ for every $p$ (if $\alpha\not=0$), thus Theorem \ref{mt2} reveals that we can control initial data beyond $L^p$ (there is an enormous literature on nonlinear heat equations  with Cauchy data in $L^p$; see e.g.\ \cite{ms, vas} and the references therein).
\item[(iii)] The results in Theorem \ref{mt2} look interesting because they are in sharp contrast to what happens for the standard heat equation with the above nonlinearity and $\lambda=1$ (focusing case), where for real constant initial data (hence in $M^{\infty,1}$), even small, one has always blow-up in finite time for every $k\not=0$, as one sees at once by solving the ordinary differential equation $u_t=u^{k+1}$, $u$ real (again the literature in this connection is large; see  \cite{vas} for a comprehensive survey).
\item[(iv)] In Theorem \ref{mt2} we suppose $p,q\geq 1$. In fact, we are interested in well-posedness in the lowest regularity/biggest spaces; moreover, the proof relies on the Minkowski integral inequality. However, the problem of the persistence of regularity in quasi-Banach spaces, as well as weighted variants of the above results, although not of primary interest, could be worth investigating.
\end{itemize}

\end{remark}

The proof of Theorem \ref{mt2} can be adapted (in fact simplified) to prove local well-posedness for \eqref{HE} in the same spaces without any smallness assumption on the initial data. We leave the details to the interested reader and we limit ourselves to briefly state the result as follows
(with the necessary clarifications in the case $p=\infty$, as in Theorem \ref{mt2}).
\begin{theorem}\label{LW}
Let $ 1\leq p \leq \infty$,  $ 1\leq q \leq \frac{2k+1}{2k}$ and $\frac{1}{q}+\frac{\beta}{kd}>1.$  Then \eqref{HE}  is locally well-posed in $M^{p,q}$ and $W^{q, p}$.
\end{theorem}

We observe that this result implies, in particular, local well-posedness in $M^{p,1}$, $1\leq p\leq\infty$. Actually this special case, when $0<\beta\leq 1$, follows directly from \cite[Thm.\  1.1]{nicola} - it was also re-obtained in \cite{EC} (in fact, if $\beta\leq 1$, the operator $H^\beta$ is a pseudodifferentual operator with a real-valued Weyl symbol, bounded from below and with bounded derivatives of order $\geq 2$, therefore satisfying the assumptions in \cite[Thm.\ 1.1]{nicola}; see Proposition \ref{prop shubin} below). However  Theorem \ref{LW} applies to every $\beta>0$  and for a range of $q>1$, allowing more singular initial data (in contrast, $M^{\infty,1}$ contains only continuous functions).\par
Modulation spaces have been widely applied in the study of nonlinear PDEs. The local and global well-posedness  for the  heat equation associated to   Laplacian  in weighted modulation spaces  goes back to the  work of Iwabuchi  \cite{ib}. In \cite{FD}, authors have proved ill-posedness for the fractional heat equation  and in \cite[Thm.\ 1.1]{bhimani1}  finite time blow-up  has been established in some modulation spaces.  On the other hand, Bhimani et al.\ proved in \cite{bhimani3}   global well-posedness for the Hartree-Fock equations associated to harmonic oscillator in some modulation spaces; see also \cite{RM}.  There is a large literature dealing with the analysis of PDEs on modulation spaces; we refer to the surveys \cite{BOP,RSW} and the monograph \cite{wang}, and the references therein (see also \cite{BO,CNR}); we also mention the article \cite{kt} for results on the Hermite operator obtained using phase-space methods. However, the study of nonlinear {\it global} well-posedness in modulation and Wiener amalgam spaces is very limited when the corresponding linear propagator is not a Fourier multiplier and, in fact, new interesting phenomena can occur, as observed in Remark \ref{fe} {\it (iii)}.
\par\bigskip
In short, the paper is organised as follows. We collect some background material on modulation spaces and a number of preliminary results in Section 2. Section 3 is devoted to the proof of Theorem \ref{mainthm}, while in Section 4 we provide the proof of Theorem \ref{mt2}.

\section{Preliminary results}\label{PR}
We write $|x|^2 = x\cdot x$ for $x \in \bR^d$, where $x\cdot y$ is the inner product on $\bR^d$. We denote by $\mathcal{S}(\bR^d)$ the Schwartz class of rapidly decaying smooth functions on $\bR^d$ and by $\mathcal{S}'(\bR^d)$ the space of temperate distributions. The bracket $\langle f,g\rangle$ stands for the inner product of $f,g \in L^2(\bR^d)$ as well as for the action of $f \in \mathcal{S}'(\bR^d)$ on $g \in \mathcal{S}(\bR^d)$; in both cases we assume it to be conjugate-linear in the second entry.

Given $x,\xi \in \bR^d$, the translation operator $T_x$ and the modulation operator $M_{\xi}$ are defined as
\[ T_x f(y) \coloneqq f(y-x),\quad M_\xi f(y) \coloneqq e^{i \xi \cdot y} f(y), \quad f \in \mathcal{S}(\bR^d), \]
and can be extended to temperate distributions by duality. The composition $\pi(z) = M_{\xi}T_x$, $z = (x,\xi) \in \bR^{2d}$, is referred to as a time-frequency shift.

The short-time Fourier transform (STFT) of $f \in \mathcal{S}'(\bR^d)$ with respect to a window function $g \in \mathcal{S}(\bR^d)\setminus\{0\}$ is defined by
\[ V_g f(x,\xi) \coloneqq \langle f,\pi(x,\xi)g \rangle = \int_{\bR^d} e^{-i \xi \cdot y}f(y)\overline{g(y-x)}dy. \]

Modulation spaces were introduced by Feichtinger \cite{F1} in the '80s. They consist of functions enjoying suitable summability/decay conditions on the phase-space side. Consider a weight $m(x,\xi)$ in phase space, i.e.\ a continuous and strictly positive function in $\bR^d\times\bR^d$ with at most polynomial growth - in fact, we will often use the polynomial weight $v_s(x,\xi) \coloneqq (1+|x|+|\xi|)^s$, $(x,\xi) \in \bR^d\times\bR^d$, $s\in\mathbb{R}$ \footnote{Actually we need a further technical assumption which will be always verified in the following - hence the reader could ignore this issue - namely $m$ is required to be $v_s$-moderate for some $s\geq0$; see \cite{gro book}.}. Let $0< p, q \le \infty$ and $g \in \mathcal{S}(\bR^d)\setminus\{0\}$; the modulation space $M^{p,q}_m(\bR^d)$ is the set of all $f\in \mathcal{S}'(\bR^d)$ such that
\begin{equation}\label{modsp norm} \| f \|_{M^{p,q}_{m}} \coloneqq \| V_g f \, m \|_{L^{p,q}} = \left( \int_{\bR^d} \left( \int_{\bR^d} |V_g f(x,\xi) m(x,\xi)|^p dx \right)^{q/p} d\xi \right)^{1/q} < \infty, \end{equation}
with obvious modifications in the case where $p=  \infty$ or $q=\infty$. When $s = 0$ we simply write $M^{p,q}$. Here we used the notation $L^{p,q}(\bR^d \times \bR^d)$ for the mixed-norm Lebesgue spaces, with (quasi-)norm $\|F(x,\xi)\|_{L^{p,q}}= \|\|F(x,\xi)\|_{L^p_x}\|_{L^q_\xi}$.

It turns out that modulation spaces are quasi-Banach spaces (Banach spaces if $p,q\geq 1$) whose definition is independent of the choice of the window function $g$ - in the sense that different choices of the window provide equivalent norms; see \cite{CR,galperin,gro book} for proofs and further details. Here we observe that they have a number of relations with standard function spaces of harmonic analysis, the most notable being that $M^{2,2}(\bR^d) = L^2(\bR^d)$, and the so-called Shubin-Sobolev, alias Hermite-Sobolev, spaces $Q^s$, $s\in\mathbb{R}$ \cite{shubin}, which can be defined \cite[Thm.\ 2.1]{gramchev} as the space of $f\in\mathcal{S}'(\bR^d)$ such that
\begin{equation}\label{eq char shusob}
\|f\|^2_{Q^s}\coloneqq \|H^{s/2}f\|^2_{L^2}=\sum_{k=0}^{+\infty} ||P_k f||^2_{L^2}(2k+d)^{s}<\infty.
\end{equation}
It is well known that, for every $0<p,q\leq\infty$, if $s$ is large enough,
\begin{equation}\label{eq emb shusob}
Q^s\hookrightarrow M^{p,q}\hookrightarrow M^\infty\hookrightarrow Q^{-s}\end{equation} (this follows easily from the characterisation $Q^s=M^{2,2}_{v_s}$ \cite[Lem.\ 4.4.19]{CR}, H\"older's inequality and the inclusion relations \cite[Thm.\ 2.4.17]{CR}).

We recall that reversing the order of integration in \eqref{modsp norm}, namely
\begin{equation}\label{was norm} \| f \|_{W^{p,q}} \coloneqq \left( \int_{\bR^d} \left( \int_{\bR^d} |V_g f(x,\xi)|^p d\xi \right)^{q/p} dx \right)^{1/q}, \end{equation} (here we take $m=1$ for simplicity) gives rise to a norm that characterizes the so-called Wiener amalgam spaces $W^{p,q}$. In fact, they are strictly related to modulation spaces via the Fourier transform, since $W^{p,q} = \mathcal{F} M^{p,q}$.  We stress that such spaces can be equivalently characterized by decomposition methods as Wiener amalgams with local component $\mathcal{F} L^p$ and global component $L^q$, that is $W^{p,q} = W(\mathcal{F} L^p,L^q)(\bR^d)$ - see \cite{CR,F0} for further details.

\begin{example}\label{ex}
\begin{itemize}
\item[(i)]
Let $f_\alpha(x)=|x|^{-\alpha}$, $0<\alpha<d$ (cf.\ Remark \ref{fe}). Let us show that the function $f_\alpha(x)$ is in $M^{p,q}(\mathbb R^d)$ for  $p>d/\alpha$ and $q>d/(d-\alpha)$. We just sketch the proof, leaving the details to the interested reader. \par One can estimate separately the STFT of $\chi f_\alpha$ and $(1-\chi) f_\alpha$, where $\chi$ is smooth with compact support in $\mathbb{R}^d$, $\chi=1$ in a neighborhood of the origin, taking a window $g$ with compact support. Then $V_g (\chi f)(x,\xi)=0$ if $x$ is large enough and for $x$ in a compact subset one uses $|V_g (\chi f)(x,\xi)|\lesssim |\hat{f_\alpha}|\ast_\xi |\hat{\chi}| \ast_\xi |\hat{\overline{g}}|(\xi)=:F(\xi)$. Indeed, since $\hat{f_\alpha}=c_\alpha f_{d-\alpha}$ for some $c_\alpha\in\mathbb{R}$, by invoking the Hardy-Littlewood-Sobolev inequality we have $F\in L^q(\mathbb R^d)$ for $q> \frac{d}{d-\alpha}$. On the other hand,  $(1-\chi)f_\alpha$ can be estimated using the embedding $L^p_k \hookrightarrow M^{p,q}$ which holds for $1 < p < \infty$ and $k \in \mathbb{N}$ large enough, where $L^p_k$ denotes the space of $L^p$ functions with $k$ distribution derivatives in $L^p$ \cite{kob}.
\item[(ii)]
Since, as already observed, $\hat{f_\alpha}=c_\alpha f_{d-\alpha}$, we have $f_\alpha\in W^{q,p}(\mathbb R^d)$ for the same range of exponents $p,q$.
\item[(iii)] Note that $\cos |x| \in M^{\infty,1}$ \cite[Cor.\ 15]{bgor unimod} and by the algebra property (Proposition \ref{gap}) we have $M^{p,q}\cdot M^{\infty,1}\subset M^{p,q}$, $1\le p,q \le \infty$. Thus, if $f\in M^{p,q}$ then $ f(x) \cos |x|\in M^{p,q}$. Similarly, using that $\cos |x|^2 \in W^{1,\infty}$ \cite[Thm.\ 14]{bgor unimod} and  $W^{q,p}\cdot W^{1,\infty}\subset W^{q,p}$, $1 \le p,q \le \infty$, we see that if $f$ belongs to $W^{q,p}$ then $ f(x) \cos |x|^2$ belongs to $W^{q,p}$ too.
\item[(iv)] Constant functions are in $M^{\infty,1}\subset W^{1,\infty}$.
\item[(v)] If $\chi$ is any smooth function with compact support in $\bR^d$ and $\Lambda\subset\mathbb{R}^d$ is any lattice in $\bR^d$, $q>d/(d-\alpha)$, then $\sum_{\mu \in\Lambda} f_\alpha(x-\mu)\chi(x-\mu)$ belongs to $W^{q,\infty}$, as one verifies easily.
\end{itemize}

\end{example}

Modulation spaces can be used both as symbol classes as well as environment where to study boundedness of pseudodifferential operators \cite{CR}, i.e.\ operators formally given by
\[
a^\mathrm{w} f(x)=(2\pi)^{-d}\int_{\bR^{2d}} e^{i\xi\cdot (x-y)} a\Big(\frac{x+y}{2},\xi\Big)f(y)\, dy d\xi
\]
where $a(x,\xi)$ is a function in phase space - the Weyl symbol of the operator $a^\mathrm{w}$. In fact the above integral, suitably interpreted in a weak sense, gives rise to a continuous operator $a^\mathrm{w}:\mathcal{S}(\bR^d)\to \mathcal{S}'(\bR^d)$ for any distribution symbol $a \in \mathcal{S}'(\bR^{2d})$. However, in the following we will only consider smooth symbols satisfying some growth conditions at infinity. Indeed, the relevant symbol classes in this paper are given by the so-called Shubin classes \cite{helffer,NR,shubin}: for $s\in\mathbb{R}$ we define $\Gamma^s$ as the space of functions $a\in C^\infty(\bR^{2d})$ such that
\[
|\partial^\alpha a(z)|\leq C_\alpha (1+|z|)^{s-|\alpha|}\qquad z\in\bR^{2d}
\]
for every $\alpha\in \mathbb{N}^{2d}$. This space becomes a Fr\'echet space when endowed with the obvious seminorms. \par

We state a generalized version of the Calder\'on-Vaillancourt theorem that will be used below.

\begin{theorem}\label{continuity}
Let $a\in \Gamma^{-s}$, $s\in\mathbb{R}$, and $0< p,q\leq \infty$. Then $a^\mathrm{w}: M^{p,q}\to M^{p,q}_{v_s}$ continuously, with operator norm depending only on a finite number of seminorms of $a$ in $\Gamma^{-s}$.
\end{theorem}
\begin{proof}
The desired continuity result follows from \cite[Thm.\ 3.1]{toft} specialised to the case of weights $\omega_1(x,\xi)= 1$,  $\omega_2(x,\xi)=(1+|x|+|\xi|)^s$ ($x,\xi\in\bR^d$), $\omega_0(z,w)= (1+|z|+|w|)^s$ ($z,w\in\bR^{2d}$), which gives the continuity of $a^\mathrm{w}: M^{p,q}\to M^{p,q}_{v_s}$ if $a\in M^{\infty,r}_{\omega_0}(\bR^{2d})$ for $r\leq\min\{1,p,q\}$. \par
On the other hand, if a symbol $a$ satisfies the estimates
\[
|\partial^\alpha a(z)|  \leq C_\alpha (1+|z|)^{-s}\qquad z\in\bR^{2d}
\]
(no additional decay is needed for the derivatives), then it belongs to $M^{\infty,r}_{\omega_0}(\bR^{2d})$ for any $r>0$. In fact, writing $e^{-iw\cdot y} =(1+|w|^2)^{-N}(1-\Delta_y)^N e^{-iw\cdot y}$ $(z,w,y\in\bR^{2d})$ in the formula for the short-time Fourier transform of $a$, and repeated integration by parts yield
\[
|V_g a(z,w)|\leq C_N (1+|w|)^{-2N}(1+|z|)^{-s}\qquad z,w\in\bR^{2d}
\]
for every $N\in\mathbb{N}$ ($g\in\mathcal{S}(\bR^{2d})$), which easily gives the claim.
\end{proof}

To conclude this section we recall a result about complex powers of pseudodifferential operators.

\begin{proposition}\label{prop shubin}
	Let $\beta >0$. The fractional Hermite operator $H^{\beta} = (-\Delta+|x|^2)^{\beta}$ is a pseudodifferential operator with Weyl symbol $a_{\beta} \in \Gamma^{2\beta}$. More precisely we have
	\begin{equation}\label{eq prop shubin}
	a_\beta(x,\xi) = (|x|^2+|\xi|^2)^\beta + r(x,\xi), \quad |x|+|\xi|\ge 1,
	\end{equation}  where $r \in \Gamma^{2\beta-2}$.
\end{proposition}
\begin{proof}
	The result follows from the machinery of complex powers applied to the operator $H$ (cf.\ \cite[Thm.\ 1.11.1, Thm.\ 1.11.2]{helffer}). Indeed $H$ is a positive operator with Weyl symbol $|x|^2+|\xi|^2$, which is positive globally elliptic - in the sense that $|x|^2+|\xi|^2 \ge C(1+|x|+|\xi|)^2$ for $|x|+|\xi|$ sufficiently large and some $C>0$. The desired result then follows e.g.\ from \cite[Thm.\ 4.3.6]{NR} specialized to the symbol class $S(M,\Phi,\Psi)=\Gamma^m$, namely with $M(x,\xi)=(1+|x|+|\xi|)^{m}$, $\Phi(x,\xi)=\Psi(x,\xi)=1+|x|+|\xi|$; the so-called Planck function $h(x,\xi)=\Phi(x,\xi)^{-1}\Psi(x,\xi)^{-1}=(1+|x|+|\xi|)^{-2}$, which gives the gain in the asymptotic expansions within this symbol class, is responsible for the gain in decay of the remainder $r(x,\xi)$ compared with the ``principal symbol'' $(|x|^2+|\xi|^2)^\beta$.
\end{proof}

\section{Proof of the main result}
\begin{proof}[Proof of Theorem \ref{mainthm}]
We prove the desired estimate separately in the regimes $t\geq 1$ and $0< t\leq 1$. \par\medskip
{\bf Case $t\geq 1$.}
It is sufficient to prove the following estimates:
\begin{equation}\label{uno}
\|P_k\|_{M^{p_1,q_1}\to M^{p_2,q_2}}\leq C_0(2k+d)^{s},\qquad k\in\mathbb{N}
\end{equation}
for some $s\geq 0$ and $C_0>0$, and
\begin{equation}\label{due}
\sum_{k=0}^{+\infty} e^{-t(2k+d)^\beta} (2k+d)^{s} \leq C_1 e^{-td^\beta},\qquad t\geq 1
\end{equation} for some $C_1>0$.
Let us prove \eqref{uno}.
As a consequence of the characterization \eqref{eq char shusob} and the embeddings in \eqref{eq emb shusob}, for $s$ large enough we have
\[
\|P_k\|_{M^{p_1,q_1}\to M^{p_2,q_2}}\leq C_0 \|P_k\|_{Q^{-s}\to Q^{s}}=C_0(2k+d)^{s}.
 \]
Let us now prove \eqref{due}. Since the sequence $\mathbb{N}\ni k\mapsto e^{-t(2k+d)^\beta}(2k+d)^{s}$ is decreasing for, say, $k\geq k_0$, we estimate separately
\[
\sum_{k=0}^{k_0} e^{-t(2k+d)^\beta} (2k+d)^{s} \leq C_2 e^{-td^\beta}
\]
for some $C_2>0$, and
\begin{align*}
\sum_{k=k_0+1}^{+\infty} e^{-t(2k+d)^\beta}(2k+d)^{s}&\leq \int_{k_0}^{+\infty} e^{-t(2x+d)^\beta}(2x+d)^{s}\, dx\\
&\leq  \int_{0}^{+\infty} e^{-t(2x+d)^\beta}(2x+d)^{s}\, dx\\
 &=\frac{e^{-td^\beta}}{2\beta}\int_0^{+\infty} e^{-ty} (y+d^\beta)^{-1+1/\beta+s/\beta}\, dy,
\end{align*}
where we applied the change of variable $(2x+d)^\beta=y+d^\beta$. The latter integral is decreasing in $t$, so that its value for $t\geq 1$ is not larger than that corresponding to $t=1$, and the proof is concluded.

\par\bigskip
{\bf Case $0< t\leq 1$.}
First of all we observe that, by the already mentioned inclusion relations of modulation spaces, we can limit ourselves to prove the desired result with $p_2$ replaced by $\min\{p_1,p_2\}$ and $q_2$ replaced by $\min\{q_1,q_2\}$. Hence from now on $p_2\leq p_1$ and $q_2\leq q_1$. \par
Recall from Proposition \ref{prop shubin} that $H^\beta$ is a pseudodifferential operator with a real-valued Weyl symbol $a_\beta \in \Gamma^{2\beta}$. Moreover, the machinery of the heat kernel of pseudodifferential operators applies (see e.g.\ \cite[Thm.\ 4.5.1]{NR}, the global ellipticity assumption being satisfied in view of the structure \eqref{eq prop shubin} of $a_\beta$), and the associated heat semigroup $e^{-tH^\beta}$ is therefore a pseudodifferential operator with Weyl symbol $b_t(x,\xi)$, depending on the parameter $t$, such that, for every $N\geq 0$, the symbol $t^Nb_t$ belongs to a bounded subset of $\Gamma^{-2\beta N}$ when $t$ stays in any compact subset of $[0,+\infty)$.\par Now, fix $N\in\mathbb{N}$ such that $2\beta N> 2d$; we can apply Theorem \ref{continuity} to the symbols $b_t$ and $t^N b_t$ and we obtain, with $m(x,\xi)=v_{2\beta}(x,\xi)=(1+|x|+|\xi|)^{2\beta}$,
\[
\| (1+t^Nm^N) V_g(e^{-tH^\beta} f)\|_{L^{p_1,q_1}}\leq C \|f\|_{M^{p_1,q_1}}
\]
for a constant $C$ independent of $t\in (0,1]$. Hence it is sufficient to prove that
\[
\|F\|_{L^{p_2,q_2}}\leq C(t) \| (1+t^N m^N) F\|_{L^{p_1,q_1}}
\]
with $C(t)$ as in the statement and for every measurable function $F(x,\xi)$. This follows by H\"older's inequality since, under our assumption, $1/p_2=1/p_1+1/\tilde{p}$ and  $1/q_2=1/q_1+1/\tilde{q}$ and
\[
\|(1+t^N m(x,\xi)^N)^{-1}\|_{L^{\tilde{p},\tilde{q}}}\leq \|(1+t^N(|x|+|\xi|)^{2\beta N})^{-1}\|_{L^{\tilde{p},\tilde{q}}}=C(t)
\]
as one sees by a linear change of variable. This proves the desired estimate in the regime $t\in (0,1]$.
  \end{proof}
\begin{remark}
We highlight that the estimate \eqref{uno} for $p_1=p_2=q_1=q_2$ can be improved as follows:
\begin{equation}\label{unobis}
\|P_k\|_{M^{p,p}\to M^{p,p}}\leq 1,\qquad k\in\mathbb{N},\ p\geq 1,
\end{equation}
for a suitable choice of the window implicit in the definition of the $M^{p,p}$-norm.
In fact, for $k\in\mathbb{N}$ we can write
\begin{align*}
P_k&=(2\pi)^{-1}\sum_{\ell =0}^{+\infty}  \left( \int_0^{2\pi} e^{-i\theta(2\ell+d)} e^{i\theta(2k+d)}\,d\theta \right)  \, P_\ell\\
&= (2\pi)^{-1} \int_0^{2\pi} e^{-i\theta H} e^{i\theta(2k+d)}\, d\theta,
\end{align*}
where the interchange of the sum and the integral is justified by Fubini Theorem, after interpreting the above equalities in weak sense and using that for $f,g\in\mathcal{S}(\bR^d)$ one has
\[
\sum_{\ell=0}^{+\infty} |\langle P_\ell  f,g\rangle|= \sum_{\ell=0}^{+\infty} |\langle P_\ell  f,P_\ell g\rangle|\leq \|f\|_{L^2} \|g\|_{L^2}.
\]
The estimate \eqref{unobis} then follows at once from the fact that
\[
\|e^{-i\theta H}\|_{M^{p,p}\to M^{p,p}}=1,\qquad \theta\in\mathbb{R},
\]
for a suitable choice of the window implicit in the definition of the $M^{p,p}$-norm; see the proof of \cite[Thm.\ 1.7]{bhimani2}. \par
This argument fails for $M^{p,q}$ with $p\not=q$, because in that case $e^{-i\theta H}$ is not bounded on $M^{p,q}$ (except for special values of $\theta$).\par
\end{remark}
\section{The nonlinear heat equation for $H^\beta$}
In this section we prove Theorem \ref{mt2}. We recall the following results \cite{F1,kb,wang dec}.

\begin{proposition}[Algebra property]\label{gap}  Let $m\in \mathbb N$, $m\geq1$.  Assume that $\sum_{i=1}^{m} \frac{1}{p_i}= \frac{1}{p_0}$, $\sum_{i=1}^{m} \frac{1}{q_i}= m-1+ \frac{1}{q_0}$ with $0<p_i\leq \infty, 1\leq q_i \leq \infty$  for $1\leq i \leq m.$ Then, for some $C>0$,
\[  \left \| \prod_{i=1}^{m} f_i \right\|_{M^{p_0, q_0}} \leq C  \prod_{i=1}^{m} \|f_i\|_{M^{p_i, q_i}}.\]
\end{proposition}
\begin{lemma} \label{Algebra}
		The multi-linear estimates
		\[
		\||f|^{2k}f\|_{M^{p,r}}\leq C \|f\|^{2k+1}_{M^{p,q}}		\]
		hold for $1\leq p,q,r \leq \infty,  k\in \mathbb N$, $\frac{2k+1}{q}=\frac{1}{r}+2k$.
	\end{lemma}
\begin{proof} From Proposition \ref{gap} we have
	\[ \||f|^{2k}f\|_{M^{p_0,r}} \leq C \|f\|_{M^{p,q}}^{2k+1}, \] with $\frac{2k+1}{p} = \frac{1}{p_0}$ and the conclusion follows from the embedding $M^{p_0,r} \hookrightarrow M^{p,r}$ ($p_0\leq p$).
\end{proof}

\begin{proof}[\textbf{Proof of Theorem \ref{mt2}}]
		{\it (a)} We study \eqref{HE} directly in integral form (Duhamel's principle), namely
				\begin{eqnarray} \label{HE1}
		u(t)=S(t)u_0+ \lambda \int_0^t S(t-\tau) \, [|u(\tau)|^{2k}u(\tau)] \, d \tau=:\mathcal{J}(u)
		\end{eqnarray}
	where $S(t)=e^{-t H^{\beta}}.$\par Let $p,q$ be as in the statement.
By Theorem \ref{mainthm} we have, for $t\geq0$,
\begin{equation}\label{14bis}
\|S(t) f\|_{M^{p,q}}\leq C_1 \|f\|_{M^{p,q}}
\end{equation}
for some $C_1>0$ and, for $r\geq q$,
		\begin{equation*}\label{eq2}
		\|S(t) f\|_{M^{p,q}}\leq C(t)\|f\|_{M^{p,r}},
		\end{equation*}
		where
		\begin{equation*}
		C(t)=\begin{cases}
		C_0 e^{-td^\beta} & t\geq 1\\
		C_0 t^{-\sigma}& 0<t\leq 1
		\end{cases}
		\end{equation*}
		for some $C_0>0,$ with $\sigma=\frac{d}{2\beta}\left(\frac{1}{q}-\frac{1}{r}\right).$  Since $1\leq q \leq \frac{2k+1}{2k}$ we can  choose $r\in [1, \infty]$ such that $\frac{2k+1}{q}=\frac{1}{r}+2k$.
		 Since $ \frac{1}{q}+ \frac{\beta}{kd}>1,$ we have $\sigma<1$.
By Minkowski's inequality for integrals and Lemma \ref{Algebra}, we obtain, for some constants $C_2,C_3>0$,
		\begin{eqnarray*}\label{eq3}
		&&\hspace{-.5in}\left\|\int_0^t S(t-\tau) \, [|u(\tau)|^{2k} u(\tau)] \, d \tau \right\|_{ M^{p,q}} \nonumber \\
		&& \leq \int_{0}^t \left\| S(t-\tau) \, [|u(\tau)|^{2k} u(\tau)] \right\|_{M^{p,q}} \, d \tau \nonumber \\
		&& \leq \int_{0}^t \, C(t-\tau) \, \| |u(\tau)|^{2k}u(\tau) \|_{M^{p,r}} \, d \tau \nonumber \\
		&& \leq C_2  \int_0^t C(t-\tau) \|u(\tau)\|_{M^{p,q}}^{2k+1} d\tau \nonumber \\
		&&  \leq  C_2 \|u\|^{2k+1}_{L^{\infty}([0,t], M^{p,q})} \int_0^t C(s) ds \leq C_3 \, \|u\|^{2k+1}_{L^{\infty}([0,+\infty), M^{p,q})}.
		\end{eqnarray*}
Combining this inequality with \eqref{14bis}, we have
		\begin{eqnarray} \label{eq4}
		\|\mathcal{J}(u)\|_{L^{\infty}([0,+\infty), M^{p,q})} &\leq& C_4 \big(  \|u_0\|_{M^{p,q}}+  \|u\|^{2k+1}_{L^{\infty}([0,+\infty), M^{p,q})}\big)
		\end{eqnarray}
		for some constant $C_4>0$.
For $\varepsilon>0,$ put \[
B_{\varepsilon}=\{u \in L^{\infty}([0,+\infty), M^{p,q}): \|u\|_{L^{\infty}([0,\infty), M^{p,q})} \leq \varepsilon \},\]
which is the closed ball of radius $\varepsilon,$ centred at the origin in $L^{\infty}([0,+\infty), M^{p,q}).$ Next, we show that the mapping $\mathcal{J}$ maps $B_{\varepsilon}$ into itself for suitable choice of $\varepsilon$.
		Now, if we assume $\|u_0\|_{M^{p,q}} \leq \frac{\varepsilon}{2 \, C_4}$ then
		from \eqref{eq4} we obtain for $u \in B_{\varepsilon}$
		\begin{eqnarray*}
			\|\mathcal{J}(u)\|_{L^{\infty}([0,+\infty), M^{p,q})} \leq \frac{\varepsilon}{2}+ C_4 \varepsilon^{2k+1}.
		\end{eqnarray*}
		Since $k>0$, we can choose $\varepsilon$ such that $\varepsilon^{2k} \leq \frac{1}{2 \, C_4}$ and as a consequence we have
		\begin{eqnarray*}
			\|\mathcal{J}(u)\|_{L^{\infty}([0,+\infty), M^{p,q})} \leq \frac{\varepsilon}{2} + \frac{\varepsilon}{2}=\varepsilon,
		\end{eqnarray*}
		that is, $\mathcal{J}(u) \in B_{\varepsilon}.$
Noticing the identity
\[|u|^{2k}u-|v|^{2k}v= (u-v) |u|^{2k} + v(|u|^{2k}-|v|^{2k}) \]
	 and  exploiting similar arguments as before, we obtain
		\begin{eqnarray*}
			\|\mathcal{J}(u)- \mathcal{J}(v)\|_{L^{\infty}([0,+\infty), M^{p,q})} \leq  \frac{1}{2} \, \| u-v\|_{L^{\infty}([0,+\infty), M^{p,q})},
		\end{eqnarray*}
possibly by taking $\varepsilon$ smaller. Therefore, using Banach's contraction principle, we conclude that $\mathcal{J}$ has a unique fixed point in $B_{\varepsilon}$ which is the solution of \eqref{HE1}.\par\bigskip

{\it (b)} Let us now show that when $p<\infty$ ($q<\infty$ because of the assumption $q\leq \frac{2k+1}{2k}$) the unique solution in $L^\infty([0,+\infty),M^{p,q})$ in fact is continuous in $t$, i.e.\ belongs to $C([0,+\infty),M^{p,q}$).
It is sufficient to repeat the above contraction argument with the space $C([0,+\infty),M^{p,q})$ in place of  $L^\infty([0,+\infty),M^{p,q})$, provided that the semigroup $S(t)$ is strongly continuous on $M^{p,q}$. To this end, observe that by \eqref{14bis} it is sufficient to prove that the map $t\mapsto S(t)f$ is continuous with values in $M^{p,q}$ for every $f$ in some dense subset of $M^{p,q}$. \par
Then, let us take $f\in \mathcal{S}(\mathbb{R}^d)$. We know that the semigroup $S(t)$ is strongly continuous on $L^2$. Hence for every $k\in\mathbb{N}$, the map $t\mapsto S(t) H^k  f=H^k S(t) f$ is continuous with values in $L^2$. But the seminorms $p_k(f)\coloneqq \|H^k f\|_{L^2}$, $k\in\mathbb{N}$, define an equivalent family of seminorms for $\mathcal{S}(\mathbb{R}^d)$, see  \cite{NR}. Hence the map $t\mapsto S(t) f$ is continuous with values in $\mathcal{S}(\mathbb{R}^d)$ and a fortiori when regarded as a map valued in $M^{p,q}$.
\par\bigskip

{\it (c)} We	shall now  prove the desired decay of the solution. By Theorem \ref{mainthm} we have, for $t\geq0$,
		\begin{eqnarray*}
			\|S(t)f \|_{M^{p,q}}\leq C_0 \, e^{-td^\beta} \|f\|_{M^{p,q}}
		\end{eqnarray*}
		for some $C_0>0$.
		Thus, we have
		\begin{eqnarray} \label{15bis}
		\|S(t)f\|_{X} \leq C_0 \,  \|f\|_{M^{p,q}}.
		\end{eqnarray}
		We already know from part {\it (a)} that
		\begin{equation*}\label{eq6}
		 \hspace{-.5in} e^{td^{\beta}}\left\|\int_0^t S(t-\tau) \, [|u(\tau)|^{2k} u(\tau)] \, d \tau \right\|_{M^{p,q}} \nonumber
		 \leq  C_1 e^{td^{\beta}} \int_0^tC(t-\tau) \| u(\tau)\|^{2k+1}_{M^{p,q}} d\tau
				\end{equation*}
for some $C_1>0$. To control the above integral we divide it  into two parts.  For $t\geq 0$		we		let $E_1= \{ \tau\in [0,t]: t-\tau <1\}$ and $E_2= \{ \tau\in [0,t]: t-\tau \geq 1\}$. Note that
		\[e^{\tau d^\beta} \|u(\tau)\|_{M^{p,q}}^{2k+1}= e^{-2k\tau d^\beta} (e^{\tau d^\beta}\|u(\tau)\|_{M^{p,q}})^{2k+1}.\] Hence, for some $C_2>0$,
\begin{eqnarray*}
e^{td^{\beta}}\int_{E_1} C(t-\tau) \| u(\tau)\|^{2k+1}_{M^{p,q}} d\tau & \leq & C_2 \, e^{td^{\beta}} \int_{E_1} (t-\tau)^{-\sigma} e^{-\tau d^{\beta}} e^{-2k \tau d^{\beta}} \|u\|_{X}^{2k+1} d\tau\\
& \leq & C_2 \, \|u\|_{X}^{2k+1} e^{d^{\beta}} \int_0^1 s^{-\sigma}  ds
\end{eqnarray*}
and, since $k>0$,
\begin{eqnarray*}
e^{td^{\beta}}\int_{E_2} C(t-\tau) \| u(\tau)\|^{2k+1}_{M^{p,q}} d\tau & \leq & C_2 \, e^{t d^{\beta}} \int_{E_2} e^{- (t- \tau) d^{\beta}} \|u(\tau)\|_{M^{p,q}}^{2k+1} d\tau\\
& \leq & C_2 \, \|u\|_{X}^{2k+1} \int_{0}^{+\infty} e^{-2k \tau d^\beta} d\tau.
\end{eqnarray*}
Combining these inequalities with \eqref{15bis} yields
		\begin{equation*} \label{eq7}
		\|\mathcal{J}(u)\|_{X} \leq C_3 ( \|u_0\|_{M^{p,q}}+  \|u\|^{2k+1}_{X})
		\end{equation*}
for some constant $C_3>0$. Now, repeating  similar arguments as before gives the desired result.	\par\bigskip
{\it (d)} The proof of the global well-posedness in $W^{q,p}(\mathbb R^d)$ goes as that above for the modulation spaces $M^{p,q}$. We can replace indeed $M^{p,q}$ with $W^{q,p}$ everywhere, using the algebra properties of $W^{q,p}$, analogous to Lemma \ref{Algebra}, that are \[
 \||f|^{2k}f\|_{W^{r,p}}\leq C \|f\|^{2k+1}_{W^{q,p}}
		\]
		for $1\leq p,q,r \leq \infty,  k\in \mathbb N$, $\frac{2k+1}{q}=\frac{1}{r}+2k$,
		which are in turn a consequence of the {\it convolution} properties for modulation spaces \cite{toft contp} (the Fourier transform turns convolution into pointwise multiplication and modulation spaces into Wiener amalgam spaces).

\end{proof}

\section*{Acknowledgments} \noindent D.\  G.\ B.\ is thankful to Henri Lebesgue Centre (IRMAR, Univ.\ Rennes 1) for the financial support and the excellent research facilities. D.\ G.\ B.\ is also thankful to DST-INSPIRE and TIFR CAM for the academic leave. \\ R.\ M.\ is thankful to Indian Institute of Science (C.V.\ Raman PDF, R(IA)CVR-PDF/2020/224) for the financial support and the excellent research facilities. \\ S.\ T.\ is supported by J.\ C.\ Bose Fellowship from D.\ S.\ T.,
Government of India. \\ F.\ N.\ and S.\ I.\ T.\ are members of the Gruppo Nazionale per l'Analisi Matematica, la Probabilit\`a e le loro Applicazioni (GNAMPA) of the Istituto Nazionale di Alta Matematica (INdAM). The present research was partially supported by MIUR grant “Dipartimenti di Eccellenza” 2018–2022, CUP: E11G18000350001, DISMA, Politecnico di Torino.


\section{Introduction}\label{sec_intro}

In the current paradigm of hierarchical structure formation, dark matter halos form through the gravitational collapse of the initial density perturbation and evolve via mass accretion and mergers with other halos \citep{White1978}. Since galaxies form and evolve in the potential well of dark matter halos, it is expected that the orientations of galaxies are related to those of their host halos, as well as the large-scale environment. Consequently, understanding the relation between the shapes of galaxies and halos is one of the most important parts of galaxy formation and evolution. Meanwhile, the alignment between galaxies and halos can be used to constrain the systematic errors in the weak gravitational lensing measurements \citep{Codis2015, Joa2015, Kilbinger2015, Kirk2015, Fortuna2020}.

In the past decades, numerous studies have characterized various types of alignments among galaxies, halos and the large-scale environment. The large-scale environment is commonly characterized by the tidal field to define the orientations of filaments or sheets in the cosmic web, leading to studies of the alignment of galaxies or halos with respect to their large-scale environment \citep{Zhang2009, Zhang2013, Zhang2015, Tempel2015, Hirv2017, Xia2017, Gane2018, Chen2019, Gane2019,Lee2019, Kraljic2020}. Observationally, the galaxy groups are usually adopted to study the alignment of the shapes of centrals, satellites or galaxy groups with respect to central-satellite position vector \citep{Brain2005, Pere2005, Agu2006, Yang2006, Fal2007, Kang2007, Siverd2009, Agu2010, Hao2011, LiZhigang2013, Chisari2014, Singh2015, Huang2018, Georgiou2019}, as well as the alignment between the major axes of the galaxies and those of their host groups \citep{WangY2008, Nie2010, Huang2016, West2017,Ume2018}. The observational evidence shows that the major axes of central galaxies are well aligned with their host clusters, and the alignment signals are stronger for redder and more luminous galaxies.
For example, \citet{Huang2016} measured the projected alignments between the major axes of the central galaxies and those of their host clusters using three different
galaxy shape measurement methods. The average alignment angle based on isophote shape measurement is about $33$ degrees in the two-dimensional case.

Besides the few observational measurements using galaxy clusters, a vast majority of studies in literature rely on the cosmological hydrodynamical simulations to determine the alignment between galaxies and their host halos \citep{Bett2010, Okabe2018, Okabe2020, Brain2019, Bho2020, Soussana2020, Tang2020, Tenneti2020}. There is general agreement in literature that the major axes of central galaxies are well aligned with their halos, and the alignment signals are stronger galaxies in more massive halos. However, the average alignment angle measured in the simulations spans a wide range from $30^\circ$ to $50^\circ$ in the three-dimensional case \citep[see e.g.,][]{Tenneti2014,Vell2015, Shao2016, Chisari2017}. The wide range of variations on the alignment signals from cosmological simulations may be caused by the adoption of different hydro subgrid physics modeling, different estimators in the quantification of position angle, and different
numerical techniques with or without adaptive mesh refinement \citep{Springel2010}.

The above two methods both have pros and cons. In observation, the orientations of the galaxy clusters determined from the satellite galaxy distribution may not be the exactly same as the halos, which would affect the accuracy of the alignment angle. On the other hand, the alignment angle measured in the hydrodynamical simulations apparently strongly correlates with the complex baryonic physics in the galaxy formation and evolution models. For example, \citet{Tenneti2017} found that the orientation of the stellar distribution can be affected by the angular momentum of the galactic wind used in the hydrodynamical simulation. In this paper, we propose a different method of measuring the alignment angle by connecting the observed galaxies to the dark matter
halos in the constrained $N$-body simulation of ELUCID \citep{WangHY2014}, which reasonably reproduces the large-scale environment of the observed galaxies.

The initial condition of ELUCID is extracted from the
density field of galaxies in the SDSS DR7 Main Sample \citep{WangHY2014}. The mass and positions of the halos at $z\sim0$ in the ELUCID simulation are consistent with those in the galaxy group catalogs, especially for halos with mass larger than $10^{13}\,h^{-1}{\rm M}_\odot$. Using a novel neighborhood abundance matching method, \citet{Yang2018} built up the connection between galaxies in SDSS DR7 to the subhalos in the corresponding observed region in the ELUCID simulation. Based on their galaxy-subhalo connection, we are able to make a one-to-one comparison between the projected major axes of galaxies and those of the subhalos.

The paper is organized as follows. In Section~\ref{sec_data}, we present the observational data from SDSS DR7 and the subhalo catalog from the ELUCID simulation. In Section~\ref{sec_method}, we describe the novel neighborhood abundance matching method which links galaxies in SDSS DR7 to dark matter subhalos in the ELUCID simulation and present the method to calculate the shapes of dark matter subhalos in simulation and groups in observation. In addition, we describe the methodology to quantify the various alignment signals. In Section~\ref{sec_result},  we present the alignments of galaxies with respect to their host groups in observation and their corresponding subhalos in simulation. In Section~\ref{sec_summary}, we summarize our main results.

\section{Data}\label{sec_data}

\subsection{Galaxies in the SDSS DR7}

Based on galaxy sample of the New York
University Value-Added Galaxy Catalog \citep[NYU-VAGC;][]{Blanton2005}, constructed from the SDSS DR7 \citep{Abazajian2009}, we collect a total of $639,359$ galaxies in the redshift range of $0.01 \le z \le 0.2$. We use $472,416$ groups identified from these galaxies with the adaptive halo-based group finder as in \citet{Yang2005} and \citet{Yang2007}. Most of groups only have a single member galaxy, and there are $68,170$ groups with at least two members. We refer the readers to \citet{Yang2007} for details.

Following \citet{Yang2008}, we separate the galaxies into red and blue subsamples according to the $g-r$ color cut in the absolute magnitudes,
\begin{equation}
\centering
g-r = 1.022 - 0.0651x - 0.00311x^2,
\label{eqn:color}
\end{equation}
where $x=M_r - 5\log h + 23.0$ and $M_r$ is the absolute r-band magnitude k+e corrected to $z=0.1$. Since the initial condition of the ELUCID simulation is constructed from the distribution of the galaxies in the continuous Northern Galactic Cap (NGC), in this paper we only select galaxies in the NGC region of the range $99^{\circ} < \alpha < 283^\circ$ and $-7^{\circ}<\delta < 75^{\circ}$ , where $\alpha$ and $\delta$ are the right ascension and declination, respectively. This results in a final sample of $396,069$ galaxies in the NGC.

\subsection{Subhalos in the ELUCID simulation}

The halo catalog used in this study is from the ELUCID simulation  \citep{WangHuiyuan2016}, which is a dark matter only, constrained simulation carried out in the Center for High Performance Computing of Shanghai Jiao Tong University using L-GADGET2 code, a memory-optimized  version of GADGET2 \citep{Springel2005}. The simulation evolves $3072^3$ dark matter particles in a periodic box of $500\,h^{-1}{\rm Mpc}$ on a side from redshift $z=100$ to the present epoch. The particle mass and softening length are $3.0875 \times 10^8 \,h^{-1}{\rm M}_\odot$ and $3.5 \,h^{-1}{\rm kpc}$, respectively. The cosmological parameters adopted in the simulation are  $\Omega_{\rm m} = 0.258$, $\Omega_{\rm b} = 0.044$, $\Omega_\Lambda = 0.742$, $h =0.72$, $n_s = 0.963$ and $\sigma_8 = 0.796$, where $\Omega_{\rm m}$ is the matter density, $\Omega_{\rm b}$ is the baryon density, $\Omega_\Lambda$ is the dark energy density, $h$ is the normalized Hubble constant today, $n_s$ is the power law index of the initial power spectrum, and $\sigma_8$ is the amplitude of matter fluctuations on $8\,h^{-1}{\rm Mpc}$ scales today.

Dark matter halos are identified using the standard friends-of-friends (FOF) algorithm \citep{Davis1985} with a linking length of $b=0.2$ times the mean particle separation and containing at least $20$ particles. Since the FOF algorithm may accidentally identify two independent structures linking with particle bridge as one structure, it is necessary to detect substructures inside a larger FOF halo. In this paper, a ${\tt subhalo}$ is defined as the locally overdense, self-bound substructure within a larger parent FOF halo. Using the algorithm SUBFIND designed by \citet{Springel2001}, we decomposed a given FOF halo into a set of disjoint self-bound subhalos. The center of the subhalo is defined as the position of the most bound particle, and the velocity is  defined as the mean velocities of the particles in the subhalo. Within these subhalos, the most massive one is regarded as the main subhalo. Obviously, compared with FOF halos, self-bound subhalos are more suitable to link central or satellite galaxies in observation.

\section{Methods}\label{sec_method}

\subsection{Matching between galaxies and subhalos}

Since the initial condition of the ELUCID simulation is constrained by the mass density field extracted from the galaxy distribution in observation, it is expected that the subhalo distribution at present is tightly correlated with the spatial distribution of galaxies in SDSS DR7.

\citet{Yang2018} proposed a novel neighborhood abundance matching method to link galaxies in SDSS DR7 to dark matter subhalos in the ELUCID simulation, according to the likelihood of the subhalo to be linked to the candidate galaxy \begin{equation}\label{eqn:match} P = M_{\rm sh} \exp \left( - \frac {r_p^2} {2 r_{\rm off}^2} \right) \exp \left( -\frac {\pi^2} {2 v_{\rm off}^2}  \right), \end{equation} where $r_p$ and $\pi$ are the separations between the galaxy and the subhalo in the perpendicular and parallel to the line-of-sight direction, $M_{\rm sh}$ is the mass of the subhalo under consideration, $r_{\rm off}$ and $v_{\rm off}$ are two free parameters. Note that $r_{\rm off} = \infty$ and $v_{\rm off} = \infty$ correspond to the traditional abundance matching method. Here, the parameters are set to be $r_{\rm off} = 5 \,h^{-1}{\rm Mpc}$ and $v_{\rm off} = 1000 ~{\rm km/s}$.  This results in a total of $396,069$ galaxy-subhalo pairs in the catalog.

In what follows, we use the galaxy-subhalo pairs by matching separately for central and satellite galaxies. Central galaxies are linked to the main subhalos in the host FOF halos and satellite galaxies to the other subhalos. The matching criteria can give better constraint for the luminosity (stellar mass)-subhalo mass relation. This results in a total of $396,069$ galaxy-subhalo pairs in the matching catalog. Using the galaxy-subhalo pairs, \citet{Yang2018} reproduced the satellite fraction, the conditional luminosity function(the conditional stellar mass function), and the biases of the galaxies.

In this paper, only galaxies in groups with at least two member galaxies are selected to measure the alignment signals. After this selection, there are $43,316$ central galaxies and $118,095$ satellite galaxies.

\subsection{Shape and alignment definition}

The shape of a subhalo containing $N$ dark matter particles is calculated by the simple inertia tensor \begin{equation}\label{eqn:it} I_{\alpha \beta} = m  \sum\limits_{i=1}^N x_{i,\alpha} x_{i,\beta}, \end{equation} where $m$ is the particle mass, $\alpha$ and $\beta$ are the inertia tensor indices with values of $1$, $2$ or $3$, and $x_{i,\alpha}$ is the position of the particle $i$ with respect to the center of the subhalo, which is defined as the position of the particle with the minimum potential. The axis lengths $a$, $b$ and $c$ ($a\geq b \geq c$) of the ellipsoidal subhalo are proportional to the square roots of the eigenvalues $\lambda_1$, $\lambda_2$ and $\lambda_3$ ($\lambda_1 \geq \lambda_2 \geq \lambda_3$) of the inertia tensor by $a=\sqrt{\lambda_1}$, $b=\sqrt{\lambda_2}$ and $c=\sqrt{\lambda_3}$ . The corresponding eigenvectors denote the directions of the major, middle and minor axes of the subhalos.

In order to calculate the projected orientation on the sky of the halos, the Cartesian coordinates in the simulation are transformed into the redshift $z$ and sky coordinates $\alpha$ and $\delta$, where $\alpha$ and $\delta$ are the right ascension and declination, respectively. For a subhalo at the location $\boldsymbol x$ with the three-dimensional direction of the major axis $\Delta{\boldsymbol x}$, the projected direction $\theta_{\rm H}$ on the sky can be calculated by \begin{equation}\label{eqn:project} \theta_{\rm H} = \arctan \left( \frac {\Delta \alpha \cos \delta} {\Delta \delta} \right), \end{equation} where $\Delta\alpha$ and $\Delta\delta$ are the right ascension and declination differences between the locations $\boldsymbol x$ and $\Delta \boldsymbol x$.

In observation, for each group with more than one member, the group shape is calculated using the projected two-dimensional case of Equation~\ref{eqn:it}, where $N$ is the number of the galaxies in the group, $x_{i,1} = \Delta \alpha \cos \delta$ and  $x_{i,2}= \Delta \delta $ are the $i$-th projected coordinates on the sky with the origin at the position of the central galaxy, where $\Delta\alpha$ and $\Delta\delta$ are the right ascension and declination differences between the locations of satellite galaxies and central galaxies. For each galaxy, the orientation angle $\theta_{\rm G}$ of the major axis of the galaxy on the sky is specified by the $25$ mag~arcsec$^{-2}$ isophote in the r-band.

In order to quantify the alignment signal between  the orientation of the major axis of a galaxy in observation and its corresponding subhalo in the ELUCID simulation, we calculate the normalized probability distribution of the angle $\theta$ between the two orientations as
\begin{equation}
P(\theta) = N(\theta)/\langle N_{\rm R}(\theta)\rangle,
\end{equation}
where $N(\theta)$ is the number of the galaxy-subhalo or galaxy-group pairs in each $\theta$ bin, and $\langle N_R(\theta)\rangle$ is the average number of such pairs obtained from $100$ random samples, in which the orientation of galaxies are kept fixed, but the other orientations (subhalos or groups) are randomized. The standard deviation of $P_{\rm R}(\theta) = N_{\rm R}(\theta)/\langle N_R(\theta)\rangle$ calculated from the random samples is used to assess the  significance of the deviation of $P(\theta)$ from unity. We note that $P(\theta)=1$ is the case without any alignment. Since the significance is quantified with respect to the null hypothesis, throughout the paper the error bars are plotted on top of $P(\theta)=1$ line. The angle $\theta$ is constrained in the range of $0^{\circ}\leq \theta \leq 90^{\circ}$. For two parallel orientations, $\theta = 0$, while for two perpendicular orientations, $\theta = 90$. In addition, we calculate the average angle $\langle \theta \rangle$ and $\theta_\sigma$, which is the standard deviation of $\langle \theta_{\rm R} \rangle$ of the $100$ random samples. In the absence of the alignment, $\langle \theta \rangle = 45^{\circ}$, however $\langle \theta \rangle = 45^{\circ}$ does not necessarily mean a isotropic distribution.

\section{Results}\label{sec_result}

\subsection{galaxy-group alignment in observation}

\begin{figure}
\includegraphics[width=0.5\textwidth]{gg_cen.eps}
\caption{The normalized probability distribution of the angles between
the major axes of central galaxies and the major axes of their host groups
in observation. The black solid line shows the alignment signal for a
total of $43,316$ central galaxies in the groups with member galaxies
larger than $2$. The horizontal line shows an isotropic distribution of
alignment angles, while the error bars are obtained from $100$ random
realizations. The average angle and its error are also indicated.}
\label{fig:gg_cen}
\end{figure}

\begin{figure}
\includegraphics[width=0.5\textwidth]{gg_sat.eps}
\caption{Similar to Figure~\ref{fig:gg_cen}, but for the alignments
between satellite galaxies and their host groups. }
\label{fig:gg_sat}
\end{figure}

\begin{figure}
\includegraphics[width=0.5\textwidth]{gg_cen_member.eps}
\caption{Galaxy-group alignment for central galaxies in groups with
member galaxies $N \geq 2$, $N \geq 4$, $N \geq 4$, and $N \geq 8$,
respectively. }
\label{fig:gg_member}
\end{figure}

\begin{figure}
\includegraphics[width=0.5\textwidth]{gg_rb.eps}
\caption{Galaxy-group alignment for red and blue central galaxies. The
error bars on top of the horizontal line are also shown for red and
blue subsamples, respectively. }
\label{fig:gg_rb}
\end{figure}

In this subsection, we measure the alignment signal between the major axes of the galaxies in the SDSS DR7 group catalog and those of their host groups. For each group with at least two members, we calculate the shape of the group using the projected two-dimensional case of Equation~\ref{eqn:it} based on the angular positions of satellite galaxies with respect to their central galaxies. There are a total of $43,316$ groups with more than one member for $396,069$ galaxies in the continuous region.  Figure~\ref{fig:gg_cen} shows the probability distribution $P(\theta)$, with the black line for all centrals. With an average alignment angle of $\langle \theta \rangle = 42\degree86\pm 0\degree13$ for a total of $43,316$ central galaxies, there is a clear trend that the major axes of central galaxies are preferentially aligned parallel to the major axes of their host groups. The error bars of the measurements are shown on top of the horizontal dotted line taken from $100$ realizations, in which the major axes of the groups have been randomized.

In observation, based on the redMapper cluster catalog from SDSS DR8, \citet{Huang2016} also found that the shapes of central galaxies are aligned with the shapes of their parent clusters traced by the satellite location in the clusters at redshifts $0.1<z<0.35$. Using a sample of $65$ distant galaxy clusters from HST observation at redshifts $0.19<z<1.8$, \citet{West2017} reported the similar alignment signals, which are also confirmed in cosmological hydrodynamical simulations \citep{Shao2016,Tenneti2020}.

To study the mass dependence of the alignment signal, the central galaxies in the galaxy-group pairs are then separated into three different halo mass subsamples with $\log (M_{\rm h}/\,h^{-1}{\rm M}_\odot)$ in the ranges of $(11,12)$, $(12,13)$, and $(13,\infty)$. Here for each central galaxy, the halo mass is taken from its matched main subhalo in the ELUCID simulation. The alignment signals for the subsamples are shown using different types of lines with different colors in Figure~\ref{fig:gg_cen}. There is a clear indication that the alignment is stronger for more massive halos. Similar trend are also found for the galaxy stellar mass.

We also study the alignment signals of satellite galaxies with respect to their host groups, as shown in Figure~\ref{fig:gg_sat}. Similar trends of alignment angle distribution, as well as the halo mass dependence, are found for the satellite galaxies. The alignment signals of satellite galaxies with $\langle \theta \rangle = 44\degree39\pm 0\degree07$ are weaker than those of central galaxies with $\langle \theta \rangle = 42\degree86 \pm 0\degree13$. In the following analysis, we mainly focus on the  central galaxies.

To investigate the effect of group richness, all groups are separated into different richness subsamples, as shown in Figure~\ref{fig:gg_member}. Central galaxies in richer groups have higher chance to be aligned with their host groups, confirming the result of \citet{Nie2010} using the galaxy cluster catalog from SDSS DR6.

Finally, we show the dependence of the alignment signal on the galaxy color in Figure~\ref{fig:gg_rb}, where we separate the galaxies into red and blue subsamples according to Equation~(\ref{eqn:color}). The alignment signal for the red centrals with $\langle \theta \rangle = 41\degree93\pm 0\degree16$ is much stronger than that of the blue centrals with $\langle \theta \rangle = 44\degree10\pm 0\degree18$, similar to the finding of \citet{Agu2010}.

\subsection{Alignment between galaxies in observation and subhalos in simulation}

\begin{figure}
\includegraphics[width=0.5\textwidth]{spin_axis_3d.eps}
\caption{The probability distribution of the angle $\cos\theta$
between the angular momentum and the  major (red),
middle (green), and minor (blue) axes of the subhalos in
the ELUCID simulation in three-dimensional case.}
\label{fig:ss_3d}
\end{figure}

\begin{figure}
\includegraphics[width=0.5\textwidth]{spin_axis_2d.eps}
\caption{The probability distribution of the angle $\theta$
between the projected angular momentum and the projected major (red),
middle (green), and minor (blue) axes in two-dimensional case.}
\label{fig:ss_2d}
\end{figure}

\begin{figure}
\includegraphics[width=0.5\textwidth]{sa.eps}
\caption{The average alignment angle $\langle \theta \rangle$
as a function of the inner particle fraction of the subhalos in  the
ELUCID simulation for spin-major (red), spin-middle (green), and
spin-minor (blue) alignments.}
\label{fig:sa}
\end{figure}

In this section, we investigate the alignment signal between the major axes of the SDSS galaxies and those of matching subhalos in the ELUCID simulation \citep{Yang2018}. As the satellite galaxy-group alignment is very week, in what follows, we only focus on the central-main subhalo pairs.

\subsubsection{spin-shape alignment in simulation}

In order to check the reliability of the shape calculation for subhalos in simulation, we first examine the spin-shape alignment within the subhalos themselves in the ELUCID simulation. The subhalo spin vector is defined as \begin{equation}\label{eqn:spin} \boldsymbol J = m \sum \limits_{i=1}^N \boldsymbol r_i \times \boldsymbol v_i, \end{equation} where $m$ is the dark matter particle mass, $\boldsymbol r_i$ is the position vector of the $i$-th particle with respect to the subhalo center of mass, and $\boldsymbol v_i$ is the velocity of the $i$-th particle relative to the bulk velocity of the subhalo.

Previous studies mainly focus on the spin-shape alignment in the three-dimensional case. For comparison,  Figure~\ref{fig:ss_3d} shows probability distribution of $\cos\theta$, where $\theta$ is the three-dimensional alignment angle between the major axes of subhalo shape and the subhalo spin vector. Note that two parallel orientations correspond to $\cos\theta = 1$, while two perpendicular orientations correspond to $\cos\theta =0$. As shown in Figure~\ref{fig:ss_3d}, the spin direction is preferentially parallel to the subhalo minor axis, which is in good agreement with the previous results \citep{Bett2007, Zhang2009, Chisari2017, Gane2018}.

The spin direction of the subhalo in three-dimensional case is calculated using Equation~\ref{eqn:spin}, and then projected on the sky using Equation~\ref{eqn:project}, similar to the projected shape orientations \citep{Zhang2013, Zhang2015}. Figure~\ref{fig:ss_2d} shows the alignment signals of Figure~\ref{fig:ss_3d} in the projected two-dimensional case. The spin direction has a strong tendency to be parallel to the minor axis with $\langle \theta \rangle = 36\degree99\pm 0\degree06$, and much weaker chance to be perpendicular to the major axes with $\langle \theta \rangle = 50\degree94\pm 0\degree06$.

The spin and shape orientations have been proven to vary as a function of distance from the center \citep{Bailin2005, Bett2010}. For each subhalo, we also calculate the spin and shape orientations for the inner subhalo regions consisting of the inner $1/8$, $1/4$, and $1/2$ of the subhalo particles. Figure~\ref{fig:sa} shows the average angle $\langle \theta \rangle$ of the spin-shape alignment in two-dimensional case as a function of the inner particle fraction of the subhalos in the ELUCID simulation.  As indicated in Figure~\ref{fig:sa}, the spin-shape alignment is dependent on the radial extent of the subhalo. Both the spin-minor and spin-major alignment are weaker in the inner part, compared to those of the entire subhalos. In the following analysis, we will mainly use the subhalo shape calculated from the entire subhalo particles.

\subsubsection{galaxy-subhalo alignment}

\citet{WangHY2012} generated the initial condition of the constrained ELUCID simulation based on the galaxy density field using the groups with halo masses larger than $10^{12}\,h^{-1}{\rm M}_\odot$. \citet{Tweed2017} compared the halos in reconstructed and original simulations and found that the reconstruction techniques are reliable for the most massive halos with masses larger than $5\times10^{13} \,h^{-1}{\rm M}_\odot$, where more than half of the halo particles in the original simulation are matched in the reconstructed one.

In the neighborhood abundance matching method developed by \citet{Yang2018}, for each galaxy, the likelihood of the corresponding subhalo calculated by Equation~\ref{eqn:match} is proportional to the subhalo mass in a certain small volume. In this method, \citet{Yang2018} started from the most massive galaxy to search for its corresponding subhalo. Therefore, the most massive subhalo in the searching volume is mostly linked to the most massive galaxy, leading to the fact that matching pairs of central galaxies are more reliable than those of satellite galaxies. In the following analysis, we mainly focus on the alignment signals of central galaxies.  In total, there are $43,316$ central galaxies in groups with at least two members.

\begin{figure}
\includegraphics[width=0.5\textwidth]{gh_cen.eps}
\caption{The probability distribution of the angles between the major
axes of the galaxies in observation and the projected major axes of
their embedding main subhalos from the ELUCID simulation.
The black solid
line shows the alignment signal for a total of $43,316$ central galaxies
in the galaxy-halo
pairs catalog. The horizontal line shows an isotropic distribution of
alignment angles, while the error bars are obtained from $100$ random
realizations for total galaxies. The average angle and its error are
also indicated.}
\label{fig:gh}
\end{figure}

Figure~\ref{fig:gh} shows the probability distribution of the angles between the major axes of the central galaxies from the SDSS DR7 and the projected major axes of their host subhalos in the ELUCID simulation. With an average alignment angle of $\langle \theta \rangle = 44\degree55 \pm 0\degree13 $ for $43,316$ central galaxies, there is a clear tendency that the major axes of central galaxies are preferentially parallel to those of their host subhalos. However, the alignment signals are weaker than those of galaxy-group alignment. Here the shapes of the subhalos in simulation are computed from the distribution of the entire dark matter particles in the subhalos. We have checked the alignment signals only using the inner $1/8$ particles in the subhalos, and the alignment signals are still weak between galaxies and the inner region of the subhalos.

Besides, we have checked the alignment signals using the galaxy-subhalo pairs with smaller separation $r_p$ and $\pi$ in Equation~\ref{eqn:match}. The alignment signals are slightly
stronger for the galaxy-subhalo pairs with smaller separation.  For $50\%$ cleaner galaxy-subhalo pairs with $r_p < 2.8 \,h^{-1}{\rm Mpc}$ and $\pi < 3.2 \,h^{-1}{\rm Mpc}$, the average angle of
galaxy-subhalo alignment is $44\degree15 \pm 0\degree25$.

We also study the mass dependence of the alignment signals using the subsamples separated by the embedding subhalo mass with $\log (M_{\rm h} /\,h^{-1}{\rm M}_\odot)$ in the ranges of $(11,12)$, $(12,13)$, and $(13,\infty)$. In Figure~\ref{fig:gh}, different types of lines with blue, green and red colors show the alignment signals for galaxies in different subhalo mass ranges. There is a clear indication that galaxies in more massive subhalos have stronger alignment signals. Based on the previous studies \citep{WangHY2012, Tweed2017, Yang2018}, the cross-identification method is more reliable for more massive subhalos. Besides, the shape measurements are more accurate for subhalos containing more dark matter particles. Obviously, the results of galaxies in subhalos with $\log (M_{\rm h}/ \,h^{-1}{\rm M}_\odot) \geq 13 $ are more reliable. As shown by the red line in Figure~\ref{fig:gh}, there is a clear and significant alignment signal with $\langle \theta \rangle = 43\degree66\pm 0\degree30$ for galaxies in subhalos of $\log (M_{\rm h}/ \,h^{-1}{\rm M}_\odot) \geq 13$.

\subsection{The impact of survey selection effects on the alignment signals}

\begin{figure}
\includegraphics[width=0.5\textwidth]{sub_fof_dm.eps}
\caption {The probability distribution of the angles between the projected major axes of the main subhalos and the projected major axes of FOF halos in simulation.
The black solid line shows the result for a total of $43,316$ FOF halos. Different colors correspond to the subsamples with different main subhalo mass.}\label{fig:sub_fof_dm}
\end{figure}

There are a total of $43,316$ FOF halos in the ELUCID simulation corresponding to the observed groups with at least two members.
In this section, based on the subhalo and FOF halo catalogs from the ELUCID simulation, we investigate three types of alignment signals, which include the alignment
between main subhalo shapes and the host FOF halos shapes traced by all dark matter particles, the alignment between main subhalo shapes and the satellite subhalo distribution systems,
and the alignment between main subhalo shapes and the SDSS matched satellite subhalo distribution systems.

\subsubsection{Signals for dark matter particles }

Based on the subhalo and FOF halo catalogs from the ELUCID simulation, we first calculate the alignments between the main subhalos and their parent FOF halos,
where the shapes of the FOF halos are calculated by the position distribution of all the dark matter particles in the FOF halos. The sample is separated into four subsamples according to
the main subhalo mass, resulting in $15,288$, $19,379$, $5,765$, and $471$ subhalos with the mass $\log (M_{\rm h} /\,h^{-1}{\rm M}_\odot)$ in the ranges of $(11,12)$, $(12,13)$, $(13,14)$ and $(14,\infty)$.

Figure~\ref{fig:sub_fof_dm} shows the alignments of the major axes of the main subhalos with respect to the FOF halo shapes traced by dark matter particles. With the average alignment angle $\langle \theta \rangle = 18\degree7\pm 0\degree13$, in dark matter only simulation, the major axes of the main subhalos are strongly correlated with the major axes of their parent
FOF halos. Besides, there is a mass dependence that the alignment signals are stronger for more massive halos.

\begin{figure}
\includegraphics[width=0.5\textwidth]{sub_fof.eps}
\caption{The probability distribution of the angles between the projected major axes of the main subhalos and the distribution of all the subhalos in FOF halos in simulation. The black solid line shows the result for a total of $43,316$ halos. Different colors correspond to the subsamples with different halo mass.}\label{fig:sub_fof}
\end{figure}

\subsubsection{Signals for all the subhalos }

To fully understand the alignment signals in observation and the ELUCID simulation, we then calculate the alignments between the main subhalos and their subhalo systems, where the
shapes of the subhalo systems are calculated by the position distribution of all the subhalos in the FOF halos.

Figure~\ref{fig:sub_fof} shows the alignments of the major axes of the main subhalos with respect to the distribution of all the subhalos in the FOF halos. There is a mass dependence that the alignment signals are stronger for more massive halos. For halo mass larger than $10^{14}\,h^{-1}{\rm M}_\odot$, the average alignment angle is $\langle \theta \rangle = 11\degree55\pm 1\degree05$.

Such kind of mass dependence of the alignment signal is also found in the hydrodynamical simulations \citep{Vell2015, Tenneti2020}. From the Horizon-AGN simulation,  \citet{Okabe2020} claimed that the major axes of central galaxies are aligned with their subhalo systems with the average angle $\sim 20$ degrees in the projected plane for $40$ cluster-sized halos with masses larger that $5\times 10^{13}\,h^{-1}{\rm M}_\odot$. As shown in Figure~\ref{fig:sub_fof}, the two-dimensional average alignment angle is $\langle \theta \rangle = 17\degree70\pm 0\degree32$ for halo masses of $ 13\leq \log (M_{\rm h}/\,h^{-1}{\rm M}_\odot)\leq 14$, which are in agreement with the results of \citet{Okabe2020}. This agreement is expected in the fact that the galaxy-halo alignment in hydrodynamical simulation is determined, to first order, by the alignment between the inner dark matter distribution and the entire dark matter in the FOF halos. \citet{Vell2015} has claimed that, for stars and dark matter enclosed in sphere of the same radius, the orientation of the stellar distribution follows that of dark matter (see the right panel of their Figure~$8$), while dark matter itself changes orientation with increasing radii, resulting in the misalignment between galaxies and halos in hydrodynamical simulations.

\begin{figure}
\includegraphics[width=0.5\textwidth]{sh_cen.eps}
\caption{The probability distribution of the angles between the projected major axes of the main subhalos and their hosts calculated using the subhalos linked with satellites
in the groups from SDSS DR7.}\label{fig:sh_cen}
\end{figure}

\begin{figure}
\includegraphics[width=0.5\textwidth]{sh_member.eps}
\caption{Similar to Figure~\ref{fig:sh_cen}, but for the host system richness $N \geq 2$, $N \geq 4$, $N \geq 4$, and $N \geq 8$, respectively. }\label{fig:sh_member}
\end{figure}

\subsubsection{Signal for the observed subhalos}

In simulation, the alignments between main subhalos and their host halos are commonly much stronger than the alignments between galaxies and their host groups in observation, which may be due to the fact that the satellites subhalos resolved in N-body or hydrodynamical simulations are much more complete than the observation. To explore the impact of the survey selection effects, we use only the subhalos that are associated with the satellite galaxies in SDSS observation to measure the shape and orientation of subhalo systems.

In observation, the shape of the group is calculated by the distribution of the satellites in the group. For consistency, here the shape of the subhalo system in simulation is also calculated using the position distribution of the subhalos that are linked with the satellites in the observed group, which makes sure that the shape calculation of the subhalo system in simulation has used the same number $N$ in Equation~\ref{eqn:it} as that of the group in observation.

Figure~\ref{fig:sh_cen} shows the distribution of the angles between the projected major axes of
the main subhalos and their subhalo systems, which are calculated only using the subhalos linked with
satellites in the observed group.
Similar to the galaxy-group alignments, there is a clear indication that the major axes of the main subhalos are preferentially aligned to the major axes of their subhalo systems in simulation. The alignments between the main subhalos and the SDSS matched satellite subhalo systems in
simulation with $\langle \theta \rangle = 42\degree08\pm 0\degree13$ are slightly larger than the galaxy-group alignments in observation with
$\langle \theta \rangle = 42\degree86\pm 0\degree13$. As shown by the red lines in Figure~\ref{fig:sh_cen},
there is a clearly significant alignment signal with $\langle \theta \rangle = 37_\cdot^\circ 91\pm 0\degree30$ for $6,236$ halos with mass larger than $10^{13}\,h^{-1}{\rm M}_\odot$, again somewhat larger than the signals obtained from the observations.

We also investigate the dependence of the alignment signals on the richness. The host systems are separated into different subsamples according to the richness at least $2$, $4$, $6$, and $8$. Figure~\ref{fig:sh_member} shows the alignment signal as a function of the richness. Similar to the tendency of the galaxy-group alignment, the alignment signals are stronger in richer systems. For host systems with the richness $N \geq 8$,  there is a clearly significant alignment with $\langle \theta \rangle = 39\degree40\pm 0\degree51$, which is in good agreement with the galaxy-group alignment with $\langle \theta \rangle = 39\degree98\pm 0\degree51$ with $N \geq 8$.

The results demonstrate that using only a small fraction of galaxies to trace the shape of host halo will significantly reduce the alignment signals between the orientations of the main subhalo and FOF host halo. Given that the ELUCID simulation can only reproduce the most massive halos, this results in a further
reduction of the alignment signals.

Taken the results shown in Figure~\ref{fig:sh_cen} as a benchmark for the case of central (main subhalo) - group shape (as traced by subhalos) alignment in the situation without baryonic effect from a dark matter only simulation, the lower alignment signals shown in Figure 1 highlight the importance of taking into account the baryon effect on the intrinsic alignment signals, as discussed in the works of \citet{Tenneti2017} and \citet{Soussana2020}.

\section{Summary}\label{sec_summary}

In this paper, we have studied the orientation of galaxies relative to their host groups in observation and their corresponding subhalos in the ELUCID simulation.

Observationally, from the group catalog, we select $43,316$ groups with member galaxies larger than $2$.  The orientations of the galaxies are defined by the position angles of the major axes specified by the $25$ mag~arcsec$^{-2}$ isophote in the r-band. The shapes of the groups are chracterized by the inertia tensor of member galaxies in the groups. Based on the $43,316$ groups from SDSS DR7, we have investigated the alignment between the major axes of galaxies and their host groups. We have found that the major axes of central and satellite galaxies have a tendency to be aligned with the major axes of their host  groups. The galaxy-group  alignment signals of satellite galaxies are weaker than those of central galaxies. There is a mass and richness dependence that the alignment signals are stronger for galaxies in massive halos and in richer groups.   In addition, the galaxy-group alignment is found to depend on the galaxy colors. Red centrals show a stronger alignment than that of blue centrals.

For 43,316 central galaxies in groups with member galaxies larger than $2$, we have matched 43,316 main subhalos from the ELUCID simulation using a novel neighborhood abundance matching method \citep{Yang2018}. From the ELUCID simulation,  we have calculated the shapes of the halos by the inertia tensor of thousands of dark matter particles within the halos.  The shapes of the halos are then projected on the sky using Equation~\ref{eqn:project}, in order to compare with the shapes of the galaxies in observation. Using 43,316 main subhalos matched to central galaxies in observation, we have examined the alignment between the major axes of galaxies and their corresponding subhalos. We find that central galaxies are  preferentially parallel to the major axes of their corresponding subhalos. Galaxies in more massive subhalos have stronger galaxy-subhalo alignment signals.

For $43,316$ main subhalos matched to central galaxies in observation, we have calculated the alignments between main subhalos and their host systems in simulation. The shapes of the host systems in simulation are calculated using the positions of the subhalos matched to galaxies in observation. Totally, the alignments between main subhalos and SDSS matched subhalo systems
in simulation are slightly stronger than galaxy-group alignments in observation. Similar to the galaxy-group alignments in observation, the projected major axes of the main subhalos are
aligned to the major axes of their subhalo systems in simulation.

Totally, the major axes of central galaxies, groups and halos are preferentially parallel each other. The alignment signals are stronger for galaxies in more massive halos.  Especially for $6,236$ central galaxies with their corresponding halo mass larger than $10^{13}\,h^{-1}{\rm M}_\odot$, there are clearly significant alignment signals with the average angles $\langle \theta \rangle = 40\degree70\pm 0\degree30$, $\langle \theta \rangle = 43\degree66\pm 0\degree30$, and $\langle \theta \rangle = 37\degree91\pm 0\degree30$  of the galaxy-group, galaxy-subhalo, and main subhalo-subhalo system alignments, respectively.

In addition, we have examined the spin-shape alignments within the subhalos themselves in the ELUCID simulation. The spin vectors of the subhalos are calculated by the angular momenta of the subhalos. The major (minor) axes of the subhalos are preferentially perpendicular (parallel) to the directions of angular momenta, which is in good agreement with the previous studies \citep{Bett2007, Zhang2009, Chisari2017, Gane2018}.  For each subhalo, we also calculate the spin-shape alignment within the inner regions consisting of the inner $1/8$, $1/4$, and $1/2$ of the subhalo particles. The spin-shape alignment signals are found to be weaker in the inner part of the subhalos.

\section*{Acknowledgements}

We thank the anonymous referee for the helpful comments that significantly improve the presentation of this paper.
This work is supported by the national science foundation of China
(Nos. 11833005,  11890692, 11621303), 111 project No. B20019 and
Shanghai Natural Science Foundation, grant No. 15ZR1446700.
We also thank the support of the Key Laboratory for Particle
Physics, Astrophysics and Cosmology, Ministry of Education.

This work is also supported by the High Performance Computing Resource
in the Core Facility for Advanced Research Computing at Shanghai
Astronomical Observatory.

Funding for the Sloan Digital Sky Survey IV has been provided by the
Alfred P. Sloan Foundation, the U.S. Department of Energy Office of
Science, and the Participating Institutions. SDSS acknowledges support
and resources from the Center for High-Performance Computing at the
University of Utah. The SDSS web site is www.sdss.org.

SDSS is managed by the Astrophysical Research Consortium for the
Participating Institutions of the SDSS Collaboration including the
Brazilian Participation Group, the Carnegie Institution for Science,
Carnegie Mellon University, the Chilean Participation Group, the
French Participation Group, Harvard-Smithsonian Center for
Astrophysics, Instituto de Astrof{\'i}sica de Canarias, The Johns
Hopkins University, Kavli Institute for the Physics and Mathematics of
the Universe (IPMU) / University of Tokyo, Lawrence Berkeley National
Laboratory, Leibniz Institut f{\"u}r Astrophysik Potsdam (AIP),
Max-Planck-Institut f{\"u}r Astronomie (MPIA Heidelberg),
Max-Planck-Institut f{\"u}r Astrophysik (MPA Garching),
Max-Planck-Institut f{\"u}r Extraterrestrische Physik (MPE), National
Astronomical Observatories of China, New Mexico State University, New
York University, University of Notre Dame, Observat{\'o}rio Nacional /
MCTI, The Ohio State University, Pennsylvania State University,
Shanghai Astronomical Observatory, United Kingdom Participation Group,
Universidad Nacional Aut{\'o}noma de M{\'e}xico, University of
Arizona, University of Colorado Boulder, University of Oxford,
University of Portsmouth, University of Utah, University of Virginia,
University of Washington, University of Wisconsin, Vanderbilt
University, and Yale University.

\section*{Data availability}
The data underlying this article will be shared on reasonable request to the corresponding author.

\bibliographystyle{mnras}