\section{Introduction}
Astronomical data from the supernova type Ia (SNeIa)
\cite{1,1b,1c,1d,1e}, cosmic microwave background (CMB)
\cite{2,2b,2c}, large scale structure (LSS) \cite{3,3b,3c}, baryon
acoustic oscillations (BAO) \cite{4} and weak lensing \cite{5}
indicate that expansion of the universe is speeding up rather than
decelerating. It is believed that the present accelerated expansion
is driven by gravitationally repulsive dominant energy component
known as ``dark energy'' (DE). Although the nature of DE remains a
mystery, various models of DE have been proposed in the literature
(for review see \cite{Padmanabhan,Padmanabhanb,Padmanabhanc}).
Theoretically, the simplest candidate for such a component is a
small positive cosmological constant, but it suffers the
difficulties associated with the fine tuning and the cosmic
coincidence problems \cite{Weinberg}.

More recently, a new DE model called ghost DE (GDE) has been
motivated from the Veneziano ghost of choromodynamics (QCD)
\cite{Urban,Urbanb,Urbanc}. The Veneziano ghost is required to exist
for the resolution of the $U(1)_A$ problem in low energy effective
theory of QCD. The ghosts make no contribution in the flat Minkowski
space, but once they are in the curved space or time-dependent
background such as our Friedmann-Robertson-Walker (FRW) universe,
the cancelation of their contribution to the vacuum energy is
off-set, leaving a small energy density $\sim\Lambda^3_{\rm QCD}H$,
where $H$ is the Hubble parameter and $\Lambda_{\rm QCD}\sim100$ MeV
is the QCD mass scale \cite{Witten,Wittenb,Wittenc,Wittend,Wittene}.
This small contribution can play an important role in the
evolutionary behavior of the universe. For instance, taking
$H\sim10^{-33}$eV at the present, $\Lambda^3_{\rm QCD}H$ gives the
right order of observed magnitude of the DE density. This
coincidence is remarkable and implies that the GDE model gets rid of
fine tuning problem \cite{Urban,Urbanb,Urbanc}. In addition, the
appearance of the QCD scale could be relevant for a solution to the
cosmic coincidence problem, as it may be the scale at which dark
matter (DM) forms \cite{Forbes}. The other advantage of the GDE with
respect to other DE models include the fact that it can be
completely explained within the standard model and general
relativity (GR), without recourse to any new field, new degree(s) of
freedom, new symmetries or modifications of GR. It is worth to note
that the GDE model does not violate unitarity, causality, gauge
invariance and other important features of renormalizable quantum
field theory, as advocated in
\cite{Zhitnitsky1,Zhitnitsky1b,Zhitnitsky1c}.

The ordinary GDE density is given by \cite{Urban,Urbanb,Urbanc}
\begin{equation}
\rho_{D}=\alpha H,\label{GDE}
\end{equation}
where $\alpha$ is a constant with dimension $[{\rm energy}]^3$, and
roughly of order of $\Lambda_{\rm QCD}^3$. This new kind of DE model
has got a lot of enthusiasm recently in the literature
\cite{Cai,sheikh,sheikhb,sheikhc,sheikhd,sheikhe,sheikhf,Sheykhi1,Sheykhi2,Sheykhi3,Rozas-Fernandez}.
Although the ordinary GDE model is consistent with the observational
data, it suffers from the difficulty to describe the early evolution
of the universe. This motivated Cai et al. \cite{Cai2} to introduce
the modified GDE (MGDE) density as
\begin{equation}\label{rho}
    \rho_{D}=\alpha H+\beta H^{2},
\end{equation}
where the constant $\alpha$ is same as that defined in the ordinary
GDE and $\beta$ is another constant with dimension $[{\rm
energy}]^2$. Note that the Veneziano ghost field in QCD is
originated from the vacuum energy of quantum fields which is of the
form $H+{\mathcal O}(H^2)$ \cite{Zhitnitsky2}. However, in the
ordinary GDE model, only the leading term $H$ has been considered.
Although the subleading term $H^2$ is too small and cannot drive the
universe to accelerating expansion, Cai et al. \cite{Cai2} showed
that this term can play an important role in the early evolution of
the universe, acting as an early DE. Using the joint analysis of the
astronomical data, Cai et al. \cite{Cai2} found that the subleading
term $H^2$ of the DE density--the early DE--could have a fraction
energy density around $10\%$. Indeed, the term proportional to $H^2$
is related to the difference between the vacuum energies in
Minkowski space and in a FRW universe \cite{Maggiore,Maggioreb}. On
the other hand, the vacuum energy difference from the Veneziano
ghost field introduced in order to solve the so-called $U(1)_A$
problem in QCD has the exact form $\alpha H+\beta H^{2}$, where
$\alpha\sim\Lambda_{\rm QCD}^3\sim (100~{\rm MeV})^3$. For other
motivations to consider this form of DE, see again \cite{Cai2},
where the MGDE model has been fitted with current observational data
including SNeIa, BAO, CMB, big bang nucleosynthesis (BBN), Hubble
parameter and growth rate of matter perturbation in order to get
some constraints on the model parameters. It was found that the MGDE
model, without having the two fundamental cosmological puzzles, like
the $\Lambda$CDM fit the astronomical data very well.

On the other hand, as is well known, an alternative proposal for DE
is scalar field scenarios such as quintessence, tachyon, K-essence
and dilaton (for review see \cite{Copeland} and references therein).
The scalar field models are an effective description of an
underlying theory of DE and can alleviate the fine tuning and cosmic
coincidence problems \cite{Ali}. Scalar fields naturally arise in
particle physics including supersymmetric field theories and
string/M theory. Therefore, scalar field is expected to reveal the
dynamical mechanism and the nature of DE. Although fundamental
theories such as string/M theory do provide a number of possible
candidates for scalar fields, they do not uniquely predict their
field $\phi$ and potential $V(\phi)$ \cite{JPWu}. Therefore it
becomes meaningful to reconstruct both the dynamics and potentials
of the scalar fields from some DE models such as holographic
\cite{Zhang11,Zhang11b,Zhang11c,Zhang11d,Zhang11e,Zhang11f,Zhang11g},
agegraphic \cite{Zhang15,Zhang15b,Zhang15c,Zhang15d} and ordinary
ghost \cite{Sheykhi2,Sheykhi3,Rozas-Fernandez}.

Describing the DE model in a scalar field framework provides a more
fundamental representation of the dark component. This motivates us
to establish different scalar field models according to evolutionary
behavior of the MGDE scenario. To do so, in section 2 we investigate
the MGDE in a spatially non-flat FRW universe. In section 3 we study
the growth of structure formation in our model. In section 4 we
reconstruct both the dynamics $\phi$ and potentials $V(\phi)$ of the
quintessence, tachyon, K-essence and dilaton scalar field models of
DE according to the evolution of MGDE density. Section 5 is devoted
to conclusions.
\section{The Veneziano MGDE and DM}

Within the framework of Einstein gravity, we consider a spatially
non-flat FRW universe filled with MGDE density $\rho_D$ and DM
energy density $\rho_{\rm m}$. Therefore, the first Friedmann
equation reads
\begin{equation}
H^2+\frac{k}{a^2}=\frac{1}{3M_P^2}~ (\rho_{D}+\rho_{\rm
m}),\label{eqfr}
\end{equation}
where the scalar curvature $k=0,1,-1$ denote a flat, closed and open
FRW universe, respectively. Also $M_{P}=(8\pi G)^{-1/2}$ is the
reduced Planck mass.

Using the fractional energy densities
\begin{equation}
\Omega_{\rm m}=\frac{\rho_{\rm m}}{\rho_{\rm cr}}=\frac{\rho_{\rm
m}}{3M_P^2H^2},~~~\Omega_{D}=\frac{\rho_{D}}{\rho_{\rm
cr}}=\frac{\rho_{D}}{3M_P^2H^2},~~~\Omega_{k}=\frac{k}{a^2H^2},
\label{eqomega}
\end{equation}
the Friedmann equation (\ref{eqfr}) can be rewritten as
\begin{equation}
1+\Omega_{k}=\Omega_{D}+\Omega_{\rm m}.\label{eq10}
\end{equation}
Substituting Eq. (\ref{rho}) into the middle relation of Eq.
(\ref{eqomega}) gives
\begin{equation}\label{OmegaDMGDE}
 \Omega_{D}=\frac{\alpha}{3M^{2}_{P}H}+1-\gamma,
\end{equation}
where
\begin{equation}
\gamma=1-\frac{\beta}{3M_{P}^{2}}.
\end{equation}
Using Eq. (\ref{OmegaDMGDE}), the curvature energy density parameter
takes the form
\begin{equation}\label{omegak}
    \Omega_{k}=\frac{\Omega_{k_0}}{(\Omega_{D_0}+\gamma-1)^{2}}\left(\frac{\Omega_{D}+\gamma-1}{a}\right)^2,
\end{equation}
where we take $a_0=1$ for the present time and the subscript ``0''
denotes the present value of a quantity.

Following the observational evidences we extend our study to the
case in which the MGDE has a bulk viscosity property
\cite{Zimdahl,Zimdahlb,Zimdahlc} and also interact with DM
\cite{Bertolami,Bertolamib}. It is well known that in the framework
of FRW metric, the shear viscosity has no contribution in the energy
momentum tensor, and the bulk viscosity behaves like an effective
pressure. Because, the CMB does not indicate significant anisotropy
due to shear viscosity and only bulk viscosity is taken into account
for the fluids in the cosmological context \cite{Jaffe}. In viscous
cosmology, shear and bulk viscosities arise in relation to space
anisotropy and isotropy, respectively. It was also pointed out that
the bulk viscosity can play a significant role in the formation of
the LSS of the universe \cite{Folomeev}. It can also alleviate
several cosmological puzzles like age problem and cosmic coincidence
problem \cite{Brevik,Brevikb,Brevikc,Brevikd}.

The interaction between DE and DM can be detected in the formation
of LSS. It was suggested that the dynamical equilibrium of collapsed
structures such as galaxy clusters would be modified due to the
coupling between DE and DM. The recent observational evidence
provided by the galaxy cluster Abell A586 supports the interaction
between DE and DM \cite{Bertolami8}. The other observational
signatures on the dark sectors' mutual interaction can be found in
the probes of the cosmic expansion history by using the SNeIa, BAO
and CMB shift data \cite{Honorez}.

In the presence of bulk viscosity and interaction, the energy
densities of MGDE and DM do not conserve separately and continuity
equations take the forms
\begin{equation}
\dot{\rho}_{D}+3H(1+\omega_{D})\rho_{D}=9H^{2}\xi-Q,\label{eqpol}
\end{equation}
\begin{equation}
\dot{\rho}_{\rm m}+3H\rho_{\rm m}=Q,\label{eqCDM}
\end{equation}
where $\omega_{D}=p_{D}/\rho_{D}$ is the equation of state (EoS)
parameter of MGDE. Also $\xi=\epsilon H^{-1}\rho_{D}$ is the bulk
viscosity coefficient with the viscosity constant $\epsilon$
\cite{Brevik,Brevikb,Brevikc,Brevikd} and
$Q=3b^{2}H(\rho_{D}+\rho_{\rm m})$ is the interaction term with the
coupling constant $b^2$ \cite{Kim06}. This expression for the
interaction term $Q$ was first introduced in the study of the
suitable coupling between a quintessence scalar field and a
pressureless cold DM field \cite{Amendola1}.

Taking time derivative of Eq. (\ref{eqfr}) and using Eqs.
(\ref{eqomega}), (\ref{eq10}), (\ref{eqpol}) and (\ref{eqCDM}) one
can get
\begin{equation}
\frac{\dot{H}}{H^2}=-\frac{3}{2}\left[1+\frac{\Omega_k}{3}+(\omega_{D}-3\epsilon)\Omega_{D}\right].\label{HdotH2}
\end{equation}
Taking time derivative of Eq. (\ref{rho}), using (\ref{HdotH2}), and
substituting the obtained result into Eq. (\ref{eqpol}) gives the
EoS parameter of the interacting viscous MGDE as
\begin{equation}\label{omega}
    \omega_{D}=\frac{1}{\Omega_{D}-\gamma-1}
    \left[1-\frac{\Omega_{k}}{3}+\frac{\gamma-1}{\Omega_{D}}\left(1+\frac{\Omega_{k}}{3}\right)+2b^{2}
    \left(\frac{1+\Omega_{k}}{\Omega_{D}}\right)\right]+3\epsilon.
\end{equation}
For completeness, we give the deceleration parameter
\begin{equation}
q=-1-\frac{\dot{H}}{H^2},
\end{equation}
which combined with the dimensionless density parameter $\Omega_D$
and the EoS parameter $\omega_D$ form a set of useful parameters for
the description of the astrophysical observations. Replacing Eq.
(\ref{omega}) into (\ref{HdotH2}) gives
\begin{equation}
q=\frac{\gamma+2(\Omega_{D}-1)-\Omega_k+3b^2(1+\Omega_k)}{\Omega_{D}-\gamma-1}
  .\label{q}
\end{equation}
Taking the derivative of Eq. (\ref{OmegaDMGDE}) with respect to
redshift $z=\frac{1}{a}-1$ and using Eqs. (\ref{HdotH2}) and
(\ref{omega}) gives
\begin{equation}\label{omegaprim}
\frac{{\rm d}\Omega_{D}}{{\rm
d}z}=-\left(\frac{3}{1+z}\right)\left(\frac{\Omega_{D}+\gamma-1}{\Omega_{D}-\gamma-1}\right)
\left[\Omega_{D}-1-\frac{\Omega_{k}}{3}+b^2(1+\Omega_{k})\right].
\end{equation}
Note that Eqs. (\ref{q}) and (\ref{omegaprim}) show that both the
deceleration and MGDE density parameters are independent of
viscosity constant. The viscosity appeared only in the dynamical EoS
parameter of MGDE, Eq. (\ref{omega}).

We solve the differential equation governing the MGDE density
parameter, Eq. (\ref{omegaprim}), numerically. In Fig.
\ref{OmegaD-z-b2}, variation of the MGDE density parameter
$\Omega_{D}$ versus redshift $z=\frac{1}{a}-1$ for different
coupling constants $b^2$ is plotted. Figure \ref{OmegaD-z-b2} shows
that: i) for a given $b^2$, $\Omega_D$ increases during history of
the universe. ii) At early and late times, $\Omega_D$ increases and
decreases, respectively, with increasing $b^2$.

In Fig. \ref{q-z-b2e} we plot the evolutionary behavior of the
deceleration parameter, Eq. (\ref{q}), for different $b^2$. Figure
\ref{q-z-b2e} presents that: i) the universe transitions from a
matter dominated epoch at early times to the de Sitter phase, i.e.
$q=-1$, in the future, as expected. The result for $b^2=0$ is in
agreement with that obtained by \cite{Cai2}. ii) For $b^2=0.0$,
$0.02$ and $0.04$ at $z=0.75$, $0.91$ and $1.10$, respectively, we
have a cosmic deceleration $q>0$ to acceleration $q<0$ transition
which is compatible with the observations \cite{Ishida}. iii) For a
given $z$, although $q$ decreases with increasing $b^2$, in the far
future $q$ becomes independent of $b^2$.

The variation of the EoS parameter of interacting viscous MGDE, Eq.
(\ref{omega}), for different coupling $b^2$ and viscosity $\epsilon$
constants is plotted in Figs. \ref{wD-z-b2}, \ref{wD-z-e} and
\ref{wD-z-b2e}. Figures \ref{wD-z-b2} and \ref{wD-z-e} illustrate
that: i) in the absence of viscosity ($\epsilon=0$), $\omega_D(z)$
varies from the quintessence phases ($\omega_D>-1$) to the phantom
regime ($\omega_D<-1$). ii) In the absence of interaction ($b^2=0$),
$\omega_D(z)$ varies from $\omega_D>-1$ to $\omega_D=-1$, which is
similar to the freezing quintessence model \cite{Caldwell2}. iii)
For a given $z$, $\omega_D$ decreases and increases with increasing
$b^2$ and $\epsilon$, respectively. Note that the result for
$b^2=\epsilon=0$ is in accordance with that obtained by \cite{Cai2}.
Figures \ref{wD-z-b2e} shows that the interaction and viscosity have
opposite effects on $\omega_D$.
\section{The growth of structure formation}
Here we investigate the growth rate of matter in the presence of
interaction between viscous MGDE and DM. Following \cite{book}, the
structure formation takes place in the Newtonian regime. Hence, the
equations of motion containing the Euler and Poisson equations for
DM in the non-relativistic approximation reduce to
\begin{equation}
\frac{\partial {\mathbf v}}{\partial t}+({\mathbf v}.\nabla){
\mathbf v}=-\frac{\nabla p_{\rm m}}{\rho_{\rm
m}}-\nabla\phi,\label{eqm1}
\end{equation}
\begin{equation}
\nabla^2\phi=4\pi G\rho_{\rm m},\label{eqm2}
\end{equation}
where ${\mathbf v}$ and $\phi$ are the velocity of DM and
gravitational potential, respectively. According to the linear
perturbation theory \cite{book}, substituting the perturbations
\begin{equation}
{\mathbf v}\rightarrow{\mathbf v}+\delta{\mathbf v},~~~p_{\rm
m}\rightarrow p_{\rm m}+\delta p_{\rm m},~~~\rho_{\rm
m}\rightarrow\rho_{\rm m}+\delta\rho_{\rm
m},~~~\phi\rightarrow\phi+\delta\phi,
\end{equation}
into Eqs. (\ref{eqm1}) and (\ref{eqm2}) and retaining only first
order corrections to the background variable, one can obtain
\begin{equation}
\delta\dot{H}=-2H\delta H-\frac{4\pi G}{3}\delta\rho_{\rm
m},\label{pert1}
\end{equation}
where we have used $\nabla.{\mathbf v}=3H$ and set $v_s^2=0$ for the
pressureless DM ($p_{\rm m}=0$).

Now we generalize the linear perturbation formalism introduced in
\cite{book} to the case of interaction between viscous MGDE and DM.
Note that Eq. (\ref{pert1}) still holds in the presence of
interaction and viscosity. Because the bulk viscous coefficient
$\epsilon$ and the coupling constant of interaction term $b^2$ only
appear in the continuity Eqs. (\ref{eqpol}) and (\ref{eqCDM}).

From Eqs. (\ref{rho}) and (\ref{eqCDM}), we have
\begin{equation}
\dot{\rho}_{\rm m}+3H\rho_{\rm m}=3b^2H(\alpha H+\beta H^2+\rho_{\rm
m}).\label{rhok1}
\end{equation}
Substituting the perturbations $\rho_{\rm m}\rightarrow\rho_{\rm
m}+\delta\rho_{\rm m}$ and $H\rightarrow H+\delta H$ into the above
relation yields the density perturbation equation for DM as
\begin{equation}
\delta\dot{\rho}_{\rm m}+3H(1-b^2)\delta\rho_{\rm
m}=\left[\frac{\rho_{\rm m}}{H}+3b^2(\alpha H+2\beta
H^2)\right]\delta H.\label{deltarho}
\end{equation}
In terms of the DM density contrast $\delta_{\rm m}=\delta\rho_{\rm
m}/\rho_{\rm m}$, Eq. (\ref{deltarho}) can be rewritten for $\delta
H$ as
\begin{equation}
\delta H=\frac{\dot{\delta}_{\rm m}+\frac{3b^2}{\rho_{\rm m}}(\alpha
H^2+\beta H^3)\delta_{\rm m}}{-3(1-b^2)+\frac{3b^2}{\rho_{\rm
m}}(2\alpha H+3\beta H^2)}.\label{deltaH}
\end{equation}
Taking the time derivative of Eq. (\ref{deltaH}) gives
\begin{equation}
\delta\dot{H}=\frac{{\rm I}}{3\left[1-2b^2\Big(1+\frac{2\alpha
H+3\beta H^2}{\rho_{\rm m}}\Big)\right]},\label{deltaHdot}
\end{equation}
where
\begin{eqnarray}
{\rm I}&=&\left[b^2\left(1+\frac{2\alpha H+3\beta H^2}{\rho_{\rm
m}}\right)-1\right]\ddot{\delta}_{\rm m} \nonumber\\&
+&\frac{b^2}{\rho_{\rm m}}\left[-7\alpha H^2-6\beta H^3+\frac{8\pi
G}{3}\rho_{\rm m}(\alpha+3\beta H)\right]\dot{\delta}_{\rm m}
\nonumber\\& +&\frac{b^2}{\rho_{\rm m}}\left[-3\alpha H^3+4\pi
G\rho_{\rm m}\Big(2\alpha H+3\beta H^2\Big)\right]\delta_{\rm m},
\end{eqnarray}
and from Eq. (\ref{pert1}) we have used $\dot{H}=-H^2-4\pi
G\rho_{\rm m}/3$. We also have neglected the terms higher than
${\mathcal O}(b^2)$. Using the latest observations (golden SNeIa,
the shift parameter of CMB and the BAO) and combining them with the
lookback time data we have that $b^2$ could be as large as 0.2 (see
\cite{Feng}) although a value of $b^2<0.04$ is favored.

Inserting Eqs. (\ref{deltaH}) and (\ref{deltaHdot}) into
(\ref{pert1}) and using $\delta\rho_{\rm m}=\rho_{\rm m}\delta_{\rm
m}$, one can get the evolution equation for the dimensionless DM
density perturbation $\delta_{\rm m}$ as
\begin{eqnarray}
&&~~\left\{1-b^2\left(1+\frac{2\alpha H+3\beta H^2}{\rho_{\rm
m}}\right)\right\}\ddot{\delta}_{\rm m} \nonumber\\&&
+\left\{2H+b^2\left(\frac{3\alpha H^2}{\rho_{\rm m}}-\frac{8\pi
G}{3}(\alpha+3\beta H)-2H\right)\right\}\dot{\delta}_{\rm m}
\nonumber\\&& -\left\{4\pi G\rho_{\rm
m}-b^2\left[3H^3\left(\frac{3\alpha+2\beta H}{\rho_{\rm
m}}\right)+4\pi G\Big(2\alpha H+3\beta H^2+2\rho_{\rm
m}\Big)\right]\right\}\delta_{\rm m}=0.\label{deltarhom}
\end{eqnarray}
Note that the viscosity constant $\epsilon$ does not appear in Eq.
(\ref{deltarhom}). This comes back to the fact that in our model,
the DM has not the viscosity properties (see Eq. \ref{eqCDM}) and
due to choosing a specific form for the bulk viscosity coefficient
$\xi=\epsilon H^{-1}\rho_{D}$ in Eq. (\ref{eqpol}), the viscosity
constant $\epsilon$ does not affect $\delta_{\rm m}$ in our model.

In the absence of interaction, i.e. $b^2=0$, Eq. (\ref{deltarhom})
recovers the well known relation given by \cite{book}
\begin{equation}
\ddot{\delta}_{\rm m}+2H\dot{\delta}_{\rm m}-4\pi G\rho_{\rm
m}\delta_{\rm m}=0.\label{deltamb20}
\end{equation}
Using the following definitions
\begin{equation}
\bar{\ddot{\delta}}_{\rm m}=\ddot{\delta}_{\rm
m}/H_0^2,~~~\bar{\dot{\delta}}_{\rm m}=\dot{\delta}_{\rm
m}/H_0,~~~\bar{H}=H/H_0,~~~\bar{\alpha}=\frac{\alpha}{3M_P^2H_0},\label{dimenlesspara}
\end{equation}
one can rewrite Eq. (\ref{deltarhom}) in dimensionless form as
\begin{equation}
{\mathcal{J}}_1\bar{\ddot{\delta}}_{\rm
m}+{\mathcal{J}}_2\bar{\dot{\delta}}_{\rm
m}+{\mathcal{J}}_3\delta_{\rm m}=0,\label{delta}
\end{equation}
where
\begin{equation}
{\mathcal{J}}_1=1-b^2\left[1+\frac{2\bar{\alpha}/\bar{H}+3(1-\gamma)}{\Omega_{\rm
m}}\right],
\end{equation}
\begin{equation}
{\mathcal{J}}_2=2\bar{H}+b^2\left[\bar{\alpha}\left(\frac{3}{\Omega_{\rm
m}}-1\right)+(3\gamma-5)\bar{H}\right],
\end{equation}
\begin{equation}
{\mathcal{J}}_3=-\frac{3}{2}\Omega_{\rm
m}\bar{H}^2+3b^2\bar{H}\left[\left(\frac{3\bar{\alpha}+2(1-\gamma)\bar{H}}{\Omega_{\rm
m}}\right)+\bar{\alpha}+\frac{3}{2}(1-\gamma)\bar{H}+\Omega_{\rm
m}\bar{H}\right].
\end{equation}
In terms of the growth factor $f$ defined as \cite{book}
\begin{equation}
f(z)=\frac{{\rm d}\ln \delta_{\rm m}}{{\rm d}\ln a}=-(1+z)\frac{{\rm
d}\ln \delta_{\rm m}}{{\rm d}z},\label{f}
\end{equation}
equation (\ref{delta}) can be rewritten as
\begin{equation}
-(1+z)\bar{H}^2\frac{{\rm d}f}{{\rm
d}z}+\bar{H}\left[-(1+z)\frac{{\rm d}\bar{H}}{{\rm
d}z}+\bar{H}f+\left(\frac{{\mathcal{J}}_2}{{\mathcal{J}}_1}\right)\right]f+\left(\frac{{\mathcal{J}}_3}{{\mathcal{J}}_1}\right)=0.\label{fdiff}
\end{equation}
Note that for the matter dominated universe, i.e. $H^2=\rho_{\rm
m}/3M_P^2$, solution of Eq. (\ref{deltamb20}) yields $\delta_{\rm
m}=a$. In this case, Eq. (\ref{f}) gives the growth factor $f=1$.

In general, the differential equation (\ref{fdiff}) has no
analytical solution. Hence, we need to solve it numerically. To do
so, we first need to know the dimensionless Hubble parameter
$\bar{H}(z)$. From Eq. (\ref{OmegaDMGDE}) and last relation in Eq.
(\ref{dimenlesspara}) we obtain
\begin{equation}
\bar{H}(z)=\frac{\bar{\alpha}}{\Omega_{D}(z)+\gamma-1},
\end{equation}
and
\begin{equation}
\bar{\alpha}=\Omega_{D_0}+\gamma-1.
\end{equation}
Now with the help of $\Omega_D(z)$ plotted in Fig. \ref{OmegaD-z-b2}
which has been already obtained by numerically solving Eq.
(\ref{omegaprim}), one can obtain the evolutionary behavior of the
growth factor $f(z)$ of DM. In Fig. \ref{f-z-b2}, variation of the
growth factor $f(z)$ versus redshift for different coupling
constants $b^2$ is plotted. Figure \ref{f-z-b2} shows that: i) for a
given $b^2$, $f$ decreases during history of the universe. ii) For a
given $z$, $f$ increases with increasing $b^2$. The result for
$b^2=0$ is same as that obtained by \cite{Cai2}.

Using Eq. (\ref{f}), the dimensionless DM density perturbation
$\delta_{\rm m}$ takes the form
\begin{equation}
\delta_{\rm m}(z)=\delta_{\rm
m_0}\exp{\left[-\int_0^z\frac{f(z)}{1+z}~{\rm d
}z\right]}.\label{deltamz}
\end{equation}
In Fig. \ref{delta-z-b2}, we plot the evolutionary behavior of
$\delta_{\rm m}(z)$ for different $b^2$. Figure \ref{delta-z-b2}
presents that: i) for a given $b^2$, $\delta_{\rm m}$ increases
during history of the universe. ii) For a given $z$, $\delta_{\rm
m}$ decreases with increasing $b^2$. This shows a suppression of
structure growth relative to the non-interacting case. In the
presence of interaction, there is less DM in the past (see Fig. \ref
{OmegaM-z-b2}), and this leads to a suppression in the growth of
structure. This is in agreement with the result obtained by
\cite{Caldera}. Note that the suppression is specific to the
circumstance that the interacting and non-interacting cases are
normalized to have the same parameters today, i.e. $\Omega_{\rm
m_0}$ and $\delta_{\rm m_0}$.

\section{Correspondence with scalar field models}

Here, our aim is to investigate whether a minimally coupled scalar
field with a specific action/Lagrangian can mimic the dynamics of
the Veneziano MGDE model so that this model can be related to some
fundamental theory (such as string/M theory), as it is for a scalar
field. For this task, it is then meaningful to reconstruct the
dynamics $\phi$ and potential $V(\phi)$ of a scalar field model
possessing some significant features of the underlying theory of DE,
such as the MGDE model. To do so and following the method proposed
by \cite{Starobinsky}, we establish a correspondence between the
MGDE and various scalar field models by identifying their respective
energy densities and equations of state and then reconstruct both
the dynamics and potential of the field.

\subsection{Quintessence MGDE}

 Quintessence is described by an
ordinary time dependent and homogeneous scalar field $\phi$ which is
minimally coupled to gravity, but with a particular potential
$V(\phi)$ that leads to the accelerating universe. The action for
quintessence is given by \cite{Copeland,Ratra,Ratrab}
\begin{equation}\label{1}
S=\int {\rm d}^4x\sqrt{-g}\left[
-\frac{1}{2}g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi-V(\phi)
\right].
\end{equation}
The energy density and pressure of the quintessence scalar field
$\phi$ are as follows
\begin{equation}
\rho_Q=\frac{1}{2}\dot \phi^2+V(\phi),\label{rhoQ}
\end{equation}
\begin{equation}
p_Q=\frac{1}{2}\dot \phi^2-V(\phi).\label{p q}
\end{equation}
The quintessence EoS parameter takes the form
\begin{equation}
\omega_Q=\frac{p_Q}{\rho_Q}=\frac{\dot \phi^2-2V(\phi)}{\dot
\phi^2+2V(\phi)}.\label{omegaQ}
\end{equation}
Identifying Eq. (\ref{omegaQ}) with the EoS parameter of interacting
viscous MGDE (\ref{omega}), $\omega_Q=\omega_D$, and also equating
Eq. (\ref{rhoQ}) with (\ref{rho}), $\rho_Q=\rho_D$, one can get
\begin{equation}
V(\phi)=\frac{1}{2}(1-\omega_D)\rho_D,\label{Vphi-2}
\end{equation}
\begin{equation} \dot
\phi^2=(1+\omega_D)\rho_D.\label{phidot2-2}
\end{equation}
Inserting Eqs. (\ref{rho}) and (\ref{omega}) into the above
equations, one can get the quintessence potential and kinetic energy
as
\begin{eqnarray}
V(\phi)=\frac{\alpha^{2}\Omega_{D}}{6M^{2}_{P}(\Omega_{D}+\gamma-1)^{2}}
\left[\frac{\Omega_{D}-\gamma-2+\frac{\Omega_{k}}{3}-\left(\frac{\gamma-1}{\Omega_{D}}\right)\left(1+\frac{\Omega_{k}}{3}\right)-2b^{2}
    \left(\frac{1+\Omega_{k}}{\Omega_{D}}\right)}{\Omega_{D}-\gamma-1}
    -3\epsilon\right],\label{VphiQ}
\end{eqnarray}
\begin{eqnarray}
\dot\phi^2=\frac{\alpha^{2}\Omega_{D}}{3M^{2}_{P}(\Omega_{D}+\gamma-1)^{2}}
\left[\frac{\Omega_{D}-\gamma-\frac{\Omega_{k}}{3}+\left(\frac{\gamma-1}{\Omega_{D}}\right)\left(1+\frac{\Omega_{k}}{3}\right)+2b^{2}
    \left(\frac{1+\Omega_{k}}{\Omega_{D}}\right)}{\Omega_{D}-\gamma-1}
    +3\epsilon\right].\label{phidot2Q}
\end{eqnarray}
Integrating Eq. (\ref{phidot2Q}) with respect to $a$ and using
(\ref{OmegaDMGDE}) yields the modified ghost quintessence scalar
field as
\begin{eqnarray}
\phi(a)-\phi(a_i)=\sqrt{3}~M_{P}\int^{a}_{a_i}\sqrt{\Omega_{D}}
\left[\frac{\Omega_{D}-\gamma-\frac{\Omega_{k}}{3}+\left(\frac{\gamma-1}{\Omega_{D}}\right)\left(1+\frac{\Omega_{k}}{3}\right)+2b^{2}
    \left(\frac{1+\Omega_{k}}{\Omega_{D}}\right)}{\Omega_{D}-\gamma-1}
    +3\epsilon\right]^{\frac{1}{2}}
~\frac{{\rm d}a}{a}.\label{phiQ}
\end{eqnarray}
The above integral cannot be taken analytically. But with the help
of Eq. (\ref{omegak}) and numerical solution of the differential
equation (\ref{omegaprim}) one can obtain the evolutionary behavior
of the quintessence MGDE scalar field. In Figs. \ref{PhiQ-z-b2},
\ref{PhiQ-z-e} and \ref{PhiQ-z-b2e} we plot the variation of scalar
field, Eq. (\ref{phiQ}), versus redshift. Figures \ref{PhiQ-z-b2}
and \ref{PhiQ-z-e} present that: i) for a given $b^2$ or $\epsilon$,
$\phi$ increases during history of the universe. ii) For a given
$z$, $\phi$ decreases and increases with increasing $b^2$ and
$\epsilon$, respectively. Figure \ref{PhiQ-z-b2} shows that for
$b^2=0.0$, $0.02$ and $0.04$ at $z<-0.97$, $-0.57$ and $-0.41$,
respectively, $\phi$ becomes pure imaginary, i.e. $\dot{\phi}^2<0$,
and does not show itself in Fig. \ref{PhiQ-z-b2}. For
$\dot{\phi}^2<0$ the MGDE scalar field behaves as a phantom-type
scalar field \cite{Caldwell}. Figure \ref{PhiQ-z-b2e} clarifies that
interaction and viscosity have opposite effects on $\phi$. For
$b^2=\epsilon$, the effect of viscosity is more than the
interaction.

In Figs. \ref{VQ-z-b2}, \ref{VQ-z-e} and \ref{VQ-z-b2e}, the
variation of quintessence MGDE potential, Eq. (\ref{VphiQ}), versus
redshift is presented for different $b^2$ and $\epsilon$. Figures
\ref{VQ-z-b2} and \ref{VQ-z-e} illustrate that: i) for a given $b^2$
or $\epsilon$, the quintessence MGDE potential $V(\phi)$ decreases
during history of the universe. ii) For a given $z$, $V(\phi)$
increases and decreases with increasing $b^2$ and $\epsilon$,
respectively. Figure \ref{VQ-z-b2e} shows that in the presence of
both interaction and viscosity, although for $b^2=\epsilon$ at early
times the effect of viscosity is more than the interaction, at late
times they neutralize the effect of each other.

\subsection{Tachyon MGDE}

The tachyon field has been proposed as the source of DE and may be
described by effective field theory corresponding to some sort of
tachyon condensate with an effective Lagrangian density given by
\cite{Sen,Senb,Senc,Send,Sene,Senf}
\begin{equation}
{\mathcal{L}}=-V(\phi)\sqrt{1+\partial_{\mu}\phi
\partial^{\mu}\phi},
\end{equation}
where $\phi$ is a tachyon scalar field and $V(\phi)$ is a potential
of $\phi$. The energy density and pressure of the tachyon scalar
filed are as follows \cite{Sen,Senb,Senc,Send,Sene,Senf}
\begin{equation}\label{rhoT}
    \rho_{T}=\frac{V(\phi)}{\sqrt{1- \dot{\phi^{2}}}},
\end{equation}
\begin{equation}
p_{T}=-V(\phi)\sqrt{1- \dot{\phi^{2}}}.
\end{equation}
The EoS parameter of the tachyon filed reads
\begin{equation}\label{omegaT}
\omega_{T}=\frac{p_{T}}{\rho_{T}}=\dot{\phi^{2}}-1.
\end{equation}
From Eqs. (\ref{rho}) and (\ref{rhoT}), $\rho_{D}=\rho_{T}$ gives
the kinetic energy term
\begin{equation}\label{phidot2}
\dot{\phi^{2}}=\frac{\Omega_{D}-\gamma-\frac{\Omega_{k}}{3}+\frac{\gamma-1}{\Omega_{D}}\left(1+\frac{\Omega_{k}}{3}\right)+2b^{2}
    \left(\frac{1+\Omega_{k}}{\Omega_{D}}\right)}{\Omega_{D}-\gamma-1}
    +3\epsilon.
\end{equation}
Also using Eqs. (\ref{omega}) and (\ref{omegaT}),
$\omega_{D}=\omega_{T}$ yields the tachyon potential
\begin{equation}\label{VphiT}
   V(\phi)=\frac{\alpha^{2}\Omega_{D}}{3M^{2}_{P}(\Omega_{D}+\gamma-1)^{2}}
   \left[\frac{1-\frac{\Omega_{k}}{3}+\frac{\gamma-1}{\Omega_{D}}\left(1+\frac{\Omega_{k}}{3}\right)+2b^{2}
    \left(\frac{1+\Omega_{k}}{\Omega_{D}}\right)}{1+\gamma-\Omega_{D}}-3\epsilon\right]^{1/2}.
\end{equation}
Integrating Eq. (\ref{phidot2}) with respect to $a$ and using
(\ref{OmegaDMGDE}) gives the evolutionary form of the modified ghost
tachyon scalar field as
\begin{eqnarray}\label{e}
  \phi(a)-\phi(a_i)=\frac{3M_{P}^{2}}{\alpha}\int^{a}_{a_i}(\Omega_{D}+\gamma-1)
  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\nonumber\\\times\left[\frac{\Omega_{D}-\gamma-\frac{\Omega_{k}}{3}+\frac{\gamma-1}{\Omega_{D}}\left(1+\frac{\Omega_{k}}{3}\right)+2b^{2}
    \left(\frac{1+\Omega_{k}}{\Omega_{D}}\right)}{\Omega_{D}-\gamma-1}
    +3\epsilon\right]^{\frac{1}{2}}
~\frac{{\rm d}a}{a}.\label{phiT}
\end{eqnarray}
In Figs. \ref{PhiT-z-b2}, \ref{PhiT-z-e} and \ref{PhiT-z-b2e} we
plot the variation of $\phi$, Eq. (\ref{phiT}), versus redshift for
different $b^2$ and $\epsilon$. Also the evolution of tachyon MGDE
potential, Eq. (\ref{VphiT}), is shown in Figs. \ref{VT-z-b2},
\ref{VT-z-e} and \ref{VT-z-b2e}. Figures \ref{PhiT-z-b2} to
\ref{VT-z-b2e} clarify that the scalar field and potential of the
tachyon MGDE model behave like the quintessence one (see Figs.
\ref{PhiQ-z-b2} to \ref{VQ-z-b2e}).

\subsection{K-essence MGDE}

The scalar field model known as K-essence is also used to explain
the DE. It is well known that K-essence scenarios have
attractor-like dynamics, and therefore avoid the fine tuning of the
initial conditions for the scalar field
\cite{Chiba,Chibab,Chibac,Picon3,Picon3b}. This kind of models is
characterized by non-canonical kinetic energy terms, and are
described by a general scalar field action which is a function of
$\phi$ and $\chi=\dot{\phi}^2/2$ and is given by
\cite{Chiba,Chibab,Chibac,Picon3,Picon3b}
\begin{equation}
S=\int {\rm d} ^{4}x\sqrt{-{\rm g}}~p(\phi,\chi),\label{action}
\end{equation}
where the Lagrangian density $p(\phi,\chi)$ corresponds to a
pressure with non-canonical kinetic terms as
\begin{equation}
p_K(\phi,\chi)=f(\phi)(-\chi+\chi^{2}),\label{pk}
\end{equation}
and the energy density of the K-essence field $\phi$ is
\begin{equation}
\rho_K(\phi,\chi)=f(\phi)(-\chi+3\chi^{2}).\label{rhok}
\end{equation}
One of the motivations to consider this type of Lagrangian
originates from considering low energy effective string theory in
the presence of a high order derivative terms. The EoS parameter of
the K-essence scalar field is obtained as
\begin{equation}\label{omegaK}
    \omega_{K}=\frac{p_K}{\rho_K}=\frac{\chi-1}{3\chi-1}.
\end{equation}
Following Eqs. (\ref{omega}) and (\ref{omegaK}),
$\omega_{D}=\omega_{K}$ gives
\begin{equation}\label{chiK}
\chi=\frac{2+\gamma-\frac{\Omega_k}{3}+\frac{\gamma-1}{\Omega_{D}}\Big(1+\frac{\Omega_k}{3}\Big)
+2b^2\Big(\frac{1+\Omega_k}{\Omega_{D}}\Big)-\Omega_{D}+3\epsilon(\Omega_{D}-\gamma-1)}
{4+\gamma-\Omega_k+\frac{\gamma-1}{\Omega_{D}}(3+\Omega_k)+6b^2\Big(\frac{1+\Omega_k}{\Omega_{D}}\Big)
-\Omega_{D}+9\epsilon(\Omega_{D}-\gamma-1)}.
\end{equation}
Integrating this with respect to $a$ gives the scalar field of the
K-essence MGDE as
\begin{eqnarray}
  \phi(a)-\phi(a_i)= \frac{3M_{P}^{2}}{\alpha}\int^{a}_{a_i}(\Omega_{D}+\gamma-1)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\nonumber\\\times\
\left[\frac{4+2\gamma-\frac{2\Omega_k}{3}+2\Big(\frac{\gamma-1}{\Omega_{D}}\Big)\Big(1+\frac{\Omega_k}{3}\Big)
+4b^2\Big(\frac{1+\Omega_k}{\Omega_{D}}\Big)-2\Omega_{D}+6\epsilon(\Omega_{D}-\gamma-1)}
{4+\gamma-\Omega_k+\Big(\frac{\gamma-1}{\Omega_{D}}\Big)(3+\Omega_k)+6b^2\Big(\frac{1+\Omega_k}{\Omega_{D}}\Big)
-\Omega_{D}+9\epsilon(\Omega_{D}-\gamma-1)}\right]^{\frac{1}{2}}
~\frac{{\rm d}a}{a}.\label{phiK}
\end{eqnarray}
The evolutionary behavior of this field for different $b^2$ and
$\epsilon$ is plotted in Figs. \ref{PhiK-z-b2}, \ref{PhiK-z-e} and
\ref{PhiK-z-b2e}. Figures present that the K-essence MGDE scalar
field increases during history of the universe.

\subsection{Dilaton MGDE}

The dilaton scalar field model is also an interesting attempt to
explain the origin of DE using string theory. This model appears
from a four dimensional effective low energy string action and
includes higher order kinetic corrections to the tree level action
in low energy effective string theory \cite{Gasperini,Gasperinib}.

For the dilaton scalar field, the pressure and energy density take
the forms \cite{Gasperini,Gasperinib}
\begin{equation}
    p_{D}=-\chi+ce^{\lambda\phi}\chi^{2},\label{pD}
\end{equation}
\begin{equation}\label{rhoD}
    \rho_{D}=-\chi+3ce^{\lambda\phi}\chi^{2}.
\end{equation}
Here $c$ and $\lambda$ are two positive constants and
$\chi=\dot{\phi}^2/2$ is the dilaton kinetic energy. Therefore, the
dilaton EoS parameter reads
\begin{equation}\label{omegaD}
    \omega_{D}=\frac{p_{D}}{
    \rho_{D}}=\frac{-1+ce^{\lambda\phi}\chi}{-1+3ce^{\lambda\phi}\chi}.
\end{equation}
Identifying Eq. (\ref{omega}) with (\ref{omegaD}) reduces to
\begin{equation}
    ce^{\lambda\phi}\chi=\frac{2+\gamma-\frac{\Omega_k}{3}+\frac{\gamma-1}{\Omega_{D}}\Big(1+\frac{\Omega_k}{3}\Big)
+2b^2\Big(\frac{1+\Omega_k}{\Omega_{D}}\Big)-\Omega_{D}+3\epsilon(\Omega_{D}-\gamma-1)}
{4+\gamma-\Omega_k+\frac{\gamma-1}{\Omega_{D}}(3+\Omega_k)+6b^2\Big(\frac{1+\Omega_k}{\Omega_{D}}\Big)
-\Omega_{D}+9\epsilon(\Omega_{D}-\gamma-1)}.
\end{equation}
Using $\chi=\dot{\phi}^2/2$ one can take the integral of above
equation with respect to $a$. The result yields
\begin{eqnarray}
  \phi(a)=\frac{2}{\lambda}\ln\left\{e^\frac{\lambda\phi(a_i)}{2}+\delta\int^{a}_{a_i}(\Omega_{D}+\gamma-1)
\right.~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\nonumber\\\left.\times\
\left[\frac{4+2\gamma-\frac{2\Omega_k}{3}+2\Big(\frac{\gamma-1}{\Omega_{D}}\Big)\Big(1+\frac{\Omega_k}{3}\Big)
+4b^2\Big(\frac{1+\Omega_k}{\Omega_{D}}\Big)-2\Omega_{D}+6\epsilon(\Omega_{D}-\gamma-1)}
{4+\gamma-\Omega_k+\Big(\frac{\gamma-1}{\Omega_{D}}\Big)(3+\Omega_k)+6b^2\Big(\frac{1+\Omega_k}{\Omega_{D}}\Big)
-\Omega_{D}+9\epsilon(\Omega_{D}-\gamma-1)}\right]^{\frac{1}{2}}
~\frac{{\rm d}a}{a}\right\},\label{phiD}
\end{eqnarray}
where $\delta=\frac{3M_p^2\lambda}{2\alpha\sqrt{c}}$. In Figs.
\ref{PhiD-z-b2}, \ref{PhiD-z-e} and \ref{PhiD-z-b2e}, variation of
the dilaton MGDE scalar field, Eq. (\ref{phiD}), is illustrated for
different $b^2$ and $\epsilon$. Figures clear that the scalar field
of the dilation MGDE model behaves like the K-essence one (see Figs.
\ref{PhiK-z-b2}, \ref{PhiK-z-e} and \ref{PhiK-z-b2e}).
\section{Conclusions}
Here we studied the Veneziano MGDE model in the framework of
Einstein's gravity. We considered a spatially non-flat FRW universe
filled with interacting viscous MGDE and DM. We derived a
differential equation governing the evolution of the MGDE density
parameter and solved it numerically. We also obtained the EoS
parameter of the interacting viscous MGDE and the deceleration
parameter of the universe. Moreover, using the linear perturbation
theory we investigated the evolution of growth of structure in our
model. Furthermore, using a correspondence between the interacting
viscous MGDE and quintessence, tachyon, K-essence and dilaton scalar
field models of DE we reconstructed the dynamics $\phi$ and
potentials $V(\phi)$ of the aforementioned scalar field models
according to the evolution of the MGDE density. Our numerical
results show the following.

(i) The MGDE density parameter $\Omega_D$ for a given coupling
constant $b^2$, increases when the time increases. The evolutionary
behavior of $\Omega_D$ is independent of viscosity.

(ii) The variation of the deceleration parameter $q$ shows that the
universe transitions from an early matter dominant epoch to the de
Sitter era in the future, as expected. Also $q$ like $\Omega_D$ does
not depend on viscosity.

(iii) The EoS parameter $\omega_D$ of the MGDE model in the absence
of viscosity ($\epsilon=0$), varies from the quintessence phase to
the phantom regime. Whereas in the absence of interaction ($b^2=0$),
it behaves like the freezing quintessence model. The interaction and
viscosity have opposite effects on $\omega_D$.

(iv) The evolution of DM density perturbation $\delta_{\rm m}$ in
the presence of interaction shows a suppression of growth of
structure. Since the DM density is lower in the past relative to the
non-interacting case, leading to a suppression of growth of
structure. Here the viscosity constant $\epsilon$ does not affect
$\delta_{\rm m}$, because the DM has not the viscosity properties in
our model.

(v) The quintessence and tachyon MGDE scalar fields for a given
$b^2$ or $\epsilon$, increase during history of the universe. For a
given redshift, they decrease and increase with increasing $b^2$ and
$\epsilon$, respectively. The potentials of the quintessence and
tachyon MGDE models for a given $b^2$ or $\epsilon$, decrease with
increasing the time. The interaction and viscosity have opposite
effects on both the scalar fields and potentials of the
aforementioned models.

(vi) The K-essence and dilaton MGDE scalar fields for a given $b^2$
or $\epsilon$ increase during history of the universe.

\subsection*{Acknowledgements}

The authors thank the unknown referees for very valuable comments.
The work of K. Karami has been supported financially by Center for
Excellence in Astronomy \& Astrophysics of Iran (CEAAI-RIAAM), under
research project No. 1/3076.


\section{Introduction}
\label{S-intro}
Small-amplitude propagating intensity disturbances are often observed with imaging telescopes in open/extended loop structures close to active regions on the sun \citep{berghmans1999, 2000A&A...355L..23D, 2003A&A...404L...1K} and also in polar plumes \citep{1998ApJ...501L.217D, 2010A&A...510L...2M, 2011A&A...528L...4K}. These are identified as alternating slanted ridges in the space-time maps, which are thought to be signatures of slow magneto-acoustic waves owing to some of the observed properties such as periodicities and propagation speeds \citep[][and references therein]{2009SSRv..149...65D}. Assuming them to be slow-mode waves, their contribution to the coronal heating and the acceleration of the solar wind has been worked out by several authors \citep{2000A&A...355L..23D, 2006A&A...448..763M, 1998ApJ...501L.217D, 2007A&A...463..713O}. Propagating slow waves can also be used as a tool for coronal seismology. Some authors used these waves' properties to derive coronal parameters that are difficult to obtain from direct observations \citep{2009A&A...503L..25W}. But the recent detailed analyses using the spectral line profiles from these structures reveal periodically enhanced emission in blue wings, which is indicative of quasi-periodic upflows \citep{2010ApJ...722.1013D, 2011ApJ...727L..37T}. It is also argued that a quasi-static plasma component in the line of sight may induce the slow waves to produce similar spectral signatures \citep{2010ApJ...724L.194V}. This makes the general nature of these disturbances somewhat ambiguous. \citet{2011ApJ...737L..43N} reported a case of slow-mode waves propagating on top of a continuous outflow. They found that the wave signatures become dominant at larger heights. It is important to establish the nature of these disturbances before proceeding to evaluate their importance in heating or to other applications. To understand the nature of these disturbances, we study the distribution of power at different frequency bins for different open coronal structures.
\section{Observations}
\label{S-observations}
A three-hour-long imaging sequence taken by the Atmospheric Imaging Assembly \citep[AIA;][]{2012SoPh..275...17L} on-board SDO, on June 5, 2010 starting from 10:00:11 UT, was used in this analysis. We chose two different coronal channels centred at 171~\AA, and 193~\AA, with their peak temperature responses at 0.8 MK, and 1.25 MK, respectively. Images at level 1.0 were already subjected to basic corrections that involve removal of dark current, de-spiking, flat-fielding, and bad-pixel removal. These were then processed to co-align all images in different channels to a common centre, to correct the roll angles, and to rescale the images to a common plate scale using the \verb+aia_prep.pro+ routine (version 4.10), available in SolarSoftWare (SSW). The final pixel scale is $\approx$ 0.6\arcsec\ and the cadence of the data is 12~s.

Three different open-loop structures, one isolated on-disk plume-like structure, a set of fan loops extending off-limb, and plume/interplume structures above the north pole of the sun were chosen in this study. Square field sub regions of dimensions 350, 600, and 700 pixel units, respectively, were considered to cover these three target regions. In both channels, images corresponding to the on-disk region were tracked for solar rotation and then co-aligned using the intensity cross-correlation, taking the first image as reference. Images in the off-limb locations were co-aligned directly because it was not possible to correct for solar rotation.
\section{Analysis and results}
\label{S-results}
We performed a wavelet analysis for the time series at each pixel to determine the distribution of power at different frequencies. Keeping the total length of the time series (3~hrs) in mind, we decided to eliminate periods longer than 30~min. This was accomplished by detrending the time series with a 150 point ($\approx$ 30~min) smoothed original light curve before performing the wavelet analysis. Frequencies with power above the 99\% confidence level were chosen to be significant for a white noise process \citep{1998BAMS...79...61T}. We then considered three different periodicity ranges, which we refer to as short (2~min -- 5~min), intermediate (5~min -- 12~min), and long (12~min -- 25~min), and produced power maps corresponding to these by adding the total significant power within the periodicity limits at each pixel location.
\begin{figure}
\centering
\resizebox{\hsize}{!}{\includegraphics[angle=90]{171_off-limb.eps}}
\resizebox{\hsize}{!}{\includegraphics[angle=90]{193_off-limb.eps}}
\caption{Powermaps in different periodicity ranges for the off-limb loop structure in the 171~\AA\ (top) and 193~\AA\ (bottom) channels of AIA. The intensity image is shown in the top left panels. All powermaps show the total power above the 99\% confidence level in the periodicity limits.}
\label{171_193_off-limb}
\end{figure}
Figure~\ref{171_193_off-limb} displays the original intensity image (top left panels) along with powermaps for these periodicity ranges in the 171~\AA\ (top) and 193~\AA\ (bottom) channels of AIA for the off-limb loop structure.  The powermaps clearly show that the power in longer periods is significant up to larger distances along the loop. In other words, this could mean that the longer periods travel farther compared to the shorter ones, which become damped faster. This trend is followed in both channels. This is true even for the other two structures, the on-disk plume-like structure and the north polar plume/interplume structures. The powermaps for these structures are shown in Fig.~\ref{171_193_other}. Since this behaviour is the same, it is possible that the nature/source of these disturbances is the same.
\begin{figure*}
\centering
\includegraphics[angle=90,height=7.8cm]{171_on-disk.eps}
\includegraphics[angle=90,height=7.8cm]{193_on-disk.eps}
\includegraphics[angle=90,height=7.8cm]{171_north-pole.eps}
\includegraphics[angle=90,height=7.8cm]{193_north-pole.eps}
\caption{{\it Top row:} Powermaps constructed in different periodicity ranges for the on-disk plume structure in 171~\AA\ (left) and 193~\AA\ (right) channels of AIA. The intensity image is shown in the top left panels. All powermaps show the total significant power above the 99\% confidence level in the periodicity limits. {\it Bottom row:} Same, but for the north polar region}
\label{171_193_other}
\end{figure*}

We now focus on the on-disk plume structure for a more detailed analysis. We traced a portion of the loop, constructed space-time maps and processed them to enhance the visibility of the alternating ridges.
\begin{figure*}[ht]
\centering
\begin{minipage}{0.39\textwidth}
\includegraphics[height=0.95\textwidth,angle=90]{loc_plot.eps}
\end{minipage}
\begin{minipage}{0.60\textwidth}
\includegraphics[height=0.98\textwidth,angle=90]{171_speed.eps}\\
\includegraphics[height=0.98\textwidth,angle=90]{193_speed.eps}
\end{minipage}
\caption{{\it Left:} Snapshot displaying the on-disk plume-like structure. The white lines overplotted bound the chosen loop location. {\it Right:} Processed space-time maps constructed from this loop portion in the 171~\AA\ (top) and 193~\AA\ (bottom) channels. The overplotted slanted solid lines connect the final peak positions derived from the linear fit. The slopes of these lines give the propagation speeds. The horizontal dashed lines in each panel mark the location of the row used in the wavelet analysis for periodicity estimation.}
\label{on-disk_plume}
\end{figure*}
Fig.~\ref{on-disk_plume} shows the location of the loop on a snapshot displaying the plume like structure along with the processed space-time maps in the two AIA channels. The processing was performed by subtracting a background trend constructed from a 55 point ($\approx$ 11~min) smooth average of the original and then dividing by it. The trend subtraction removes or suppresses periods longer than 11~min, and division by the trend compensates for the intensity fall along the loop. We measured the propagation speed at three clean ridges following the method described in \citet{krishna2012}. This method involves identifying the peak positions in time from the space-time map at each spatial location along the ridge and then fitting them linearly. The peak positions were identified from the locations of maxima within a strip wide enough to cover the ridge alone. The apparent propagation speed and the error in its estimation could then be derived from the slope of the fit and the goodness of the fit. The boundaries of the chosen strips and the final fitted peak positions are shown as dotted and solid lines along the
\begin{table}
\begin{center}
\caption{Periodicity and propagation speeds in different AIA channels. Second-strongest periods are also listed in brackets.}
\label{values}
\begin{tabular}{ccccc}
\hline\hline
Channel & Periodicity & \multicolumn{3}{c}{Apparent speed (km~s$^{-1}$)} \\
\cline{3-5}
 (\AA)  &  (min)      & Ridge1 & Ridge2 & Ridge3 \\
\hline
 171    & 9.4 (15.7) & 88.7$\pm$9.7  & 80.5$\pm$5.7 & 70.6$\pm$5.8 \\
 193    & 9.4 (15.7) & 100.8$\pm$6.4 & 84.0$\pm$8.6 & 85.0$\pm$13.3 \\
\hline
\end{tabular}
\end{center}
\end{table}
ridge, respectively, in the space-time maps shown in Fig.~\ref{on-disk_plume}. Unlike the case of choosing two points manually on each ridge, this method takes several points along the ridge and uses the peak positions in the speed estimation and hence improves the error in it. Recently, \citet{2012SoPh..279..427K} used a similar method to determine the speeds of such propagating disturbances and discussed their temperature dependence. \citet{2012A&A...543A...9Y} also developed three different methods for measuring the apparent speed of propagating waves with considerably high accuracy. Apparent propagation speeds obtained from the two AIA channels are listed in Table~\ref{values}. The same ridges are chosen in both channels for comparison. The observed apparent propagation speed is consistently higher in all three ridges in the hotter 193~\AA\ channel.
\begin{figure*}
\centering
\includegraphics[angle=90,height=8.0cm]{5jun_AIA171_pix50_11min_smooth.eps}
\includegraphics[angle=90,height=8.0cm]{5jun_AIA193_pix50_11min_smooth.eps}
\caption{Wavelet plots for the time series from the row marked by the dashed line in Fig~\ref{on-disk_plume} in 171~\AA\ (left) and 193~\AA\ channels. In each of these figures top panels show the original light curve in solid line and the subtracted trend in dashed line. The middle panels display the trend subtracted and normalized light curves. The main wavelet plot is shown in the bottom left and the global wavelet plot in the bottom right panels. Overplotted dotted line in the global wavelet plot is the 99\% significance curve. Top two strongest significant periods are listed as P1 and P2 in the middle right portion.}
\label{171_193_wavelet}
\end{figure*}
To estimate periodicity, we performed a wavelet analysis at an arbitrarily chosen row marked by a horizontal dashed line in the space-time maps. The wavelet plots at this location in the 171~\AA\ and 193~\AA\ channels are shown in Fig.~\ref{171_193_wavelet}. The first two significant periods obtained from this analysis are labelled at the middle right portion of these plots. These periodicity values are also listed in Table~\ref{values}. As can be seen from these values the periodicity of these disturbances is 9.4~min. The values listed in brackets in the table are the second-strongest peaks. Note the presence of the longer periods despite the 11~min filtering. They are only suppressed but not completely eliminated.
\begin{figure}
\centering
\resizebox{\hsize}{!}{\includegraphics[angle=90, clip=true]{spatial_damping.eps}}
\caption{{\it Top:} Solid line in grey represents the intensity values along the loop, relative to the background, averaged over three successive time frames in 171~\AA\ channel. Plus symbols in blue are the 10 point ($\approx$ 4.3~Mm) smooth averages of these values. The damping sine curve fitted to these smoothed values is overplotted as a red solid line. Obtained fit parameters are labelled. {\it Bottom:} Similar curves and values generated from the 193~\AA\ channel.}
\label{damping}
\end{figure}

We also investigated the spatial damping along the loop. Fig.~\ref{damping} shows the intensity values along the loop (solid grey line) relative to the background, averaged over three successive time frames. Overplotted plus symbols in blue are the 10 point ($\approx$ 4.3~Mm) smooth average values. A damping sine wave function given by
\begin{equation}
	I(l)=A_{0}\mathrm{sin}(\omega l+\phi)\mathrm{exp}(- l/L_{d}) + B_{0}+B_{1}l
 \label{eq:function}
\end{equation}
was fitted to these values, which is also overplotted as a solid line in red. Here $A_{0}$ is the amplitude of the sine wave, $l$ is the length along the loop, $\omega$ and $\phi$ are frequency and phase of the wave, and $L_{d}$ represents the damping length. The last two terms were used to fit the background trend, if present. The first-order fit is found to be more than sufficient for this purpose. The best fits were obtained using a Levenberg-Marquardt least-squares minimisation \citep{1969drea.book.....B}. Top and bottom panels of Fig.~\ref{damping} display the plots from the 171~\AA\ and 193~\AA\ channels, respectively. Obtained fit parameters are labelled in the respective panels. Clearly, the relative amplitude of the oscillation and the damping length ($L_{d}$ value in equation~\ref{eq:function}) are lower and the wavelength is slightly longer in the hotter 193~\AA\ channel. The damping lengths are 38.9~Mm and 24.5~Mm in the 171~\AA\ and 193~\AA\ channels respectively. Lower damping lengths indicate faster damping in the 193~\AA~channel.
\section{Slow-wave model}
\label{S-model}
The farther propagation of the lower frequencies and the stronger dissipation in hotter channels may be understood as an effect of the damping of slow magneto-sonic waves by thermal conduction. We consider a simple 1D model \citep{demoortel2003,owen2009}, where radial structuring has been neglected. More complicated models have been studied by \citet{soler2008}. The dispersion relation of slow waves in such a simple configuration is given by
\begin{equation}
	\omega^3-\imath \gamma d c^2_\mathrm{s} k^2 \omega^2 - \omega k^2 c^2_\mathrm{s} + \imath d k^4 c^4_\mathrm{s}=0,
\end{equation}
where $\omega$ is the angular wave frequency, $k$ is the wave number, $\gamma$ is the ratio of specific heats and $c_\mathrm{s}$ is the sound speed. The parameter $d=\frac{(\gamma -1)\kappa_\parallel T_0}{\gamma c^2_\mathrm{s} p_0}$ expresses the importance of the thermal conduction. It has the dimension of time. The quantities $p_0$ and $T_0$ are the background pressure and temperature, respectively. $\kappa_\parallel=\kappa_0 T_0^{5/2}$ is the thermal conduction coefficient. Low $d$ means a weak thermal conduction, high $d$ corresponds to strong thermal conduction. Since the observed waves are presumably driven from near the foot points of the loop, we considered the frequency $\omega$ to be fixed by the driver, and solved the dispersion equation for the wave number $k$. For the vanishing thermal conduction ($d\omega = 0$), we find $k^2=\omega^2/ c^2_\mathrm{s}$, which corresponds to the sound wave solution. In the limit of weak thermal conduction ($d\omega \ll 1$), the solution may be found as
\begin{equation}
	k=\frac{\omega}{c_\mathrm{s}} \left( 1-\imath\frac{d\omega}{2}(\gamma-1)\right)=\frac{\omega}{c_\mathrm{s}}-\imath \frac{1}{L_\mathrm{cond}}, \label{eq:damping}
\end{equation}
where we defined the thermal conduction damping length $L_\mathrm{cond}$. This formula is only marginally applicable in the solar corona, since the thermal conduction in coronal loops is quite strong \citep{vd2011}. For instance, the results from \citet{vd2011} indicate $d\omega \approx 0.7$ which is less than 1 but not very low. However, a lot of the physics relevant for these observations may be understood from this asymptotic formula.

We first consider the effect of temperature. When the temperature increases (i.e. when we observe the structure in a higher temperature passband), the thermal conduction coefficient $\kappa_\parallel$ increases as well. This results in a higher value for $d$, and thus a lower value for the damping length $L_\mathrm{cond}$ in the limiting case of weak thermal conduction ($d\omega \ll 1$) that we currently consider. As found in the observations, the slow waves are more heavily damped in higher temperature passbands. Now, we consider the frequency dependence of the damping. In Eq.~\ref{eq:damping}, we have defined $L_\mathrm{cond}=2c_\mathrm{s}/d\omega^2(\gamma -1)$, where it is now obvious that the damping length is inversely proportional to the frequency squared, or directly proportional to the period squared. Indeed, in the observations we find that the longer periods propagate much higher (longer damping lengths). This simple model explains the qualitative features of the observations quite well. This model could be tested by measuring the damping lengths for different frequency bands, and showing that a power law with slope 2 is found. In principle, it also possible to use this simple model to measure the thermal conduction. By measuring the damping length in a loop for several frequency bands, the thermal conduction may be found as the intercept of the power law between the frequency and damping lengths.

\section{Conclusions}
\label{S-conclusions}
Open-loop or plume-like structures both off-limb and on-disk show that propagating disturbances with longer periods travel farther distances in the powermaps constructed in different periodicity ranges. The apparent speeds obtained for the propagating disturbance in the on-disk plume structure are higher in the hotter (193~\AA) channel. The amplitude of the disturbance is relatively low in this channel. The disturbance is also observed to be damped along the loop. The damping lengths were found to be 38.9~Mm and 24.5~Mm in 171~\AA\ and 193~\AA\ channels, respectively. The faster damping of the disturbance in hotter channels might indicate possible damping due to thermal conduction and matches the earlier theories on slow-mode wave damping. All these observed properties can be explained using a simple slow-wave model. We conclude that these disturbances are more likely to be signatures of slow mode waves than quasi-periodic high-speed upflows. Furthermore, these intensity waves are omnipresent at longer periods, higher up even in active regions.

\begin{acknowledgements}
We thank the referee for his/her valuable comments, which have enabled us to improve the quality of the paper. The AIA data used here are courtesy of SDO (NASA) and AIA consortium. DB wishes to thank the KULeuven research council for a senior research fellowship (SF/11/005).
\end{acknowledgements}

\bibliographystyle{aa}


\section{Introduction}
Because of progressive cooling with age, brown dwarfs do not obey a unique
mass-luminosity relation \citep{1997ApJ...491..856B, 2001RvMP...73..719B}.
Therefore, the determination of a brown dwarf mass requires either a good
knowledge of its age or a direct dynamical measurement. This, in turn, is
possible for substellar companions of stars or for those that are found in
multiple systems. An additional advantage is that the metallicity can be
inferred from the primary star. For solar-type stars the atmospheres are
much better understood than for very low-mass stars and brown dwarfs, given
the poor knowledge of opacities in cool atmospheres \citep{2006ApJ...652.1604B,
2005A&A...442..635B}. Coeval systems containing low-mass companions also
provide very useful constraints on evolutionary and atmospheric models
(e.g. \citealp{2006MNRAS.368.1281P, 2010ApJ...721.1725D}) as well as
offering a rather unique view on how the process of star formation works
at the very bottom of the main sequence \citep{2007prpl.conf..427B}. In
particular, brown dwarf companions with well determined metallicities are
benchmark objects allowing for a better understanding of the effects of
metallicity on the physical properties and evolution of substellar objects
\citep{2012MNRAS.422.1922P}. Unfortunately, substellar companions located
at wide separations ($>$50\,au) from their parent stars are relatively
rare, with an estimated frequency of less than a few per cent
\citep{2004AJ....127.2871M, 2007ApJ...670.1367L, 2007ApJ...662..413K}.

We are conducting a search for very low-mass common proper motion companions
of nearby ($\apprle 25$\,pc) stars, using the VISTA Hemisphere Survey (VHS)
(McMahon et al., in preparation) and Two Micron All Sky Survey (2MASS)
\citep{2006AJ....131.1163S}. Our sample of objects includes some of the known
multiple systems (\citealp{2010AJ....139..176F} and references therein). One
of the targets investigated so far was HD\,221356, already known to be a
triple system. The F8.0 V primary is a field star with slightly subsolar
metallicity $\left[\rm{Fe/H}\right]=-0.26$ \citep{2005ApJS..159..141V},
located at $26.12\pm0.37$\,pc \citep{2007AA...474..653V}. The main properties
of this star are given in Table~\ref{tab:tab1}. \cite{2000AJ....120.1085G}
reported that the secondary, initially described as a single object, has a
photometric distance consistent with that of the HD\,221356 star determined
by {\it Hipparcos}. The secondary was later resolved by \cite{2002ApJ...567L..53C}
into a binary separated by $0.57$\,arcsec ($\sim$14.9 au), using adaptive
optics on the Gemini North Telescope. Based on their photometric colours,
they estimated spectral types of M8.0\,V and L3.0\,V for each of the two
components, respectively. This binary (hereafter referred to as HD\,221356BC)
was also investigated by \cite{2007ApJ...667..520C}. Using data at epochs
separated by 48.3 years, he confirmed a common proper motion between
HD\,221356A and BC. He also measured a mean separation of
$\rho=451 \farcs 8 \pm 0 \farcs 4$ between both components, which corresponds
to a projected physical separation of nearly twelve thousand au, making it
one of the widest known systems with an L-type component (see also fig 11 of
\citealp{2010MNRAS.404.1817Z}).

In this article we present the identification and characterization of a fourth,
very low-mass companion of the HD\,221356 system. We outline the procedure and
results of our proper motion search together with the analysis of $I$ and
$YJHK_s$-band photometry and near-infrared spectroscopic data of the identified
companion. We derive the physical properties of the new object which turns out
to be very close to the hydrogen burning mass limit.

\begin{table}
\caption{Properties of HD\,221356A (a.k.a. HR\,8931, HIP\,116106).}
\centering
\begin{tabular}{l l c}
\hline
\hline
Parameter & Value & References\\
\hline
RA. (J2000) &  $23^{\rm h} 31^{\rm m} 31 \fs 62$ & -\\
Dec. (J2000) &  $-04\degr 05\arcmin 16 \farcs 78$ & -\\
$V$ (mag) & 6.50 & 1\\
Spectral type & F8.0 V & 3 \\
$\mu_{\alpha}\cos(\delta)$ & $178.7 \pm 0.9 $ mas/yr & 2\\
$\mu_{\delta}$ & $-192.8 \pm 0.9$ mas/yr & 2\\
Parallax & 38.29 $\pm$ 0.54 mas & 2\\
Distance & 26.12 $\pm$ 0.37 pc & 2\\
T$\rm{_{eff}}$ & 5976 $\pm$ 44 K & 1 \\
$\log(g)$ & 4.31 $\pm$ 0.06 cm/s$^2$ & 1 \\
Mass & 0.94 $\pm$ 0.13 M$_{\odot}$ & 1 \\
$\left[\rm{Fe/H}\right]$ & $-0.26 \pm 0.03$ & 1 \\
Age & $2.5-7.9$ Gyr & 1\\
\hline
1) \cite{2005ApJS..159..141V} & 2) \cite{2007AA...474..653V} \\
3) \cite{2007ApJ...667..520C}  \\
\hline
\end{tabular}
\label{tab:tab1}
\end{table}

\section{Identification and follow-up observations}
\subsection{VISTA Hemisphere Survey data}
The new low-mass companion of HD\,221356\,A was identified using the 2MASS and VHS
catalogues. The VHS is a near-infrared public survey
intended to cover the entire Southern hemisphere ($\sim20,000~\rm deg^2$) in
the $JK_s$ broad band filters with a sensitivity more than 3 mag deeper than
2MASS. It uses the 4.1-m telescope VISTA (Visible and Infrared Survey Telescope
for Astronomy) operating since 2009 at ESO's Cerro Paranal Observatory in Chile
\citep{2006Msngr.126...41E}. It is equipped with a wide-field near-infrared
camera (VIRCAM), comprising sixteen 2k$\times$2k pixel detectors with a mean
plate scale of $0.339$\,arcsec. The HD\,221356 system was observed with VISTA on
2010 November 25 and 26. Average seeing conditions were $1.4^{\prime\prime}$
and $0.9^{\prime\prime}$, respectively.

The VHS near-infrared images are processed and calibrated automatically by a
dedicated science pipeline implemented by the Cambridge Astronomical Survey Unit (CASU).
Standard reduction and processing steps include dark and flat-field corrections,
sky background subtraction, linearity correction, destripe and jitter stacking.
For more detailed description we refer to the CASU webpage http://casu.ast.cam.ac.uk/surveys-projects/vista,
as well as to \cite{2004SPIE.5493..411I} and \cite{2010ASPC..434...91L}.

\subsection{Proper motion}

The search for additional common proper motion companions of HD 221356 A
was done using the astrometry given in the VHS and 2MASS catalogues, which provide
a 12.18 yr baseline. The positions of 2MASS sources have an estimated accuracy
of 70-80 mas over the magnitude range of $9<K_s\leq16$ \citep{2006AJ....131.1163S}.
The astrometric solution for VHS observations is done automatically as part of the
CASU pipeline, using the 2MASS point source catalogue. The objects on the catalogues
extracted from each VISTA detector are matched to their counterparts in 2MASS using a
correlation radius of 1 arcsec. Because 2MASS has a high degree of internal
consistency it is possible to calibrate the world coordinate system of VISTA
images to better than 0.1 arcsecond.

\begin{figure}
\centerline{
\vbox{
\hbox{
\includegraphics[keepaspectratio=true,scale=0.75]{fig1.pdf}
}}}
\caption{Proper motion vector-point diagram for the HD\,221356 system. All correlated
objects within 15 arcmin from the primary are plotted as black dots, with HD\,221356
components labelled as A, BC and D. The primary is saturated in both surveys, its
proper motion value was taken from the literature. Time baseline between the 2MASS
and VHS epochs is 12.18 years.}
\label{fig:1}
\end{figure}

The search was performed using TOPCAT,\footnote[1]
{http://www.star.bris.ac.uk/$\sim$mbt/topcat/} a useful tool for analysis
and manipulation of source catalogues and other data tables, developed as part of the Virtual Observatory.
We retrieved the astrometric and photomeric data from both
2MASS and VHS catalogues, for all the objects within a radius of 15 arcmin corresponding
to $\sim$23,000 au around the examined star. To avoid some of the spurious
detections, we selected sources brighter than $J$=17 mag in 2MASS.
We have cross-matched 300 objects from both catalogues within
1 arcsec. The sources that remained unmatched were subsequently cross
correlated taking into account the proper motion of the primary star
provided by {\it Hipparcos} \citep{2007AA...474..653V}. We illustrate
the resulting proper motion vector-point diagram of HD 221356 on Figure \ref{fig:1}.

We have found that the proper motion of HD 221356 BC (Table~\ref{tab:tab2}) is
consistent with that of the primary HD 221356 A, thereby confirming the result of the
K\"{o}nigstuhl survey of \cite{2007ApJ...667..520C}. Individual
components of the BC pair, separated by only 0.57\,arcsec \citep{2002ApJ...567L..53C},
were not resolved in the VHS images. We also identified a new common
proper motion companion (2MASS J23313095-0405234, hereafter HD\,221356\,D), with $(\mu_{\alpha}\cos\delta,~\mu_{\delta})
=(153.48\pm21,-190.20\pm19)$ mas/yr, located $12.13\pm0.18$\,arcsec southwest from
the primary. We adopted a total astrometric uncertainty of $\sim 28$ mas/yr,
estimated using the standard deviation of proper motions for sources with $\mu\,<$100 mas/yr.
The proper motion of the new object is common with that of the HD 221356 system.
Measured separations, position angles and proper motions are listed in Table~\ref{tab:tab2}.

\subsection{Photometry}

The VHS catalogue provides aperture photometry for HD\,221356 in the $YJHK_s$ near-infrared
bands. The VISTA photometric system is calibrated using the magnitudes of colour selected 2MASS
stars converted onto the VISTA system using colour equations, including terms to account
for interstellar reddening\footnote[2]{http://casu.ast.cam.ac.uk/surveys-projects/vista/technical/photometric-properties}.
Photometric calibrations are determined to an accuracy of 1-2\%.

 \begin{figure}
 \centerline{
\vbox{
\hbox{
 \hbox{\includegraphics[scale=1]{fig2.png} }
}}}
\caption{False colour VISTA $J$-band image of HD 221356AD. Angular separation
is $12.13\pm0.18$ arcsec and the position angle of the identified companion
 is $221.8\pm1.7$ degrees. Saturation in the centre of the primary is visible. Field-of-view is
1$\times$1 arcmin, with North up and East to the left.}
\label{fig:2}
\end{figure}

The new faint companion HD 221356 D is well resolved in the VHS images (see Figure
~\ref{fig:2}), but within the glare of the primary. In order to minimize the
possible light contamination in the aperture photometry of the VHS catalogue, we
applied a method to suppress the PSF of the primary by subtracting the flipped
and rotated images from the original ones. Standard deviation of the background
at the same separation as the companion was a factor two lower in the PSF
subtracted images than in the original ones, which slightly improves the
detection of the object. We performed aperture photometry
on the resultant images with an aperture of one full-width half-maximum (FWHM).
The instrumental magnitudes were then calibrated to apparent magnitudes by adding
the aperture corrections determined using the VHS photometry of fourteen isolated
bright stars located within 10 arcmin from HD 221356A. The differences between
the new and the catalogue photometry in the $Y, J, H$ and $K_s$ bands are
$0.53\pm0.06$, $0.26\pm0.04$, $0.15\pm0.03$, $0.13\pm0.03$ mag, respectively.
Photometric values for the BC companion were taken directly from the VHS source catalogue.

\begin{table} \scriptsize
\caption{Proper motion, separations and position angles of low mass components of the HD 221356 system.}
\centering
\begin{tabular}{c c c c c c c}
\hline
\hline
\multirow{2}{*}{Comp.} &  $\mu_{\alpha}\cos(\delta)$  & $\mu_{\delta}$ & $\rho$*   & $\theta$* & $r$  \\
                      &  (mas/yr)                    & (mas/yr)       & (arcsec) & (deg)    & (au) \\
\hline
 BC & $176\pm21$&$-167\pm19$ & $451.10\pm0.18$ & $261.77\pm0.04$ & $11900\pm50$ \\
 D & $154\pm21$&$-190\pm19$ & \hspace*{0.7mm}$12.13\pm0.18$ & $221.8\pm1.7$ & \hspace*{1.2mm}$317\pm9$  \\
\hline
\end{tabular}

\raggedright \hspace*{2mm} * epoch (MJD) = 55525.12460836; $\rho$, $\theta$ and $r$ are measured with respect to the primary.
\label{tab:tab2}
\end{table}

On November 15th, 2011 we performed follow-up $I$-band imaging of HD 221356 AD.
Observations were carried out using the IAC80 telescope equipped with a E2V
2048$\times$2048 back illuminated CCD detector with a plate scale of 0.304
arcsec/pix, which provides a 10.4$\times$10.4 arcmin field of view.
We selected the 12 images with best seeing (FWHM$<$1.2 arcsec) to minimize interference from the primary star.
We reduced the images applying standard techniques, including bias and
flat-field correction, using \textsc{iraf} routines. Individual exposure times were 5s. For each image, we
performed a similar method for the PSF subtraction of the primary as used
for VISTA images and subsequently aligned and combined all of them.
We obtained the PSF-fit photometry using the \textsc{daophot} package in \textsc{iraf} and calibrated the instrumental
magnitude of our object using 11 bright stars in the field with DENIS \citep{1999A&A...349..236E} $I$-band data available.
We note that the photometric system used (Cousin) is not the same as the one of DENIS,
and that some differences may appear for very cool objects, however in our previous photometric
calibrations we found that the zeropoint between IAC80 and Denis $I$-band has small colour dependence \citep{2005A&A...443.1021C}.

Additionally, we acquired $I$-band observations of HD\,221356\,A using FastCam,
mounted on the 1.5-m Carlos S\'{a}nchez Telescope at the Teide Observatory on January 31, 2012.
FastCam is a lucky imaging instrument, designed to perform high spatial
and time resolution observations \citep{2008SPIE.7014E.137O}. Optics provide
a plate scale of 43.5 mas/pix, and a field of view of $\sim22\times22$
arcsec$^2$. We obtained 17 blocks of 1000 images of 50\,ms individual exposure
times. Images were bias corrected, aligned and stacked into the final image
using the software provided by the FastCam team. We derived
the $I$-band aperture photometry of the primary, since we found that the
literature values based on photographic plates are not reliable. Instrumental
magnitudes were calibrated using photometric standard stars from \cite{1992AJ....104..340L}
observed at different airmasses along the night under photometric conditions.
We also explored the inner region to search for the presence of
additional companions, but none were identified. We may exclude the presence
of an equal mass companion to the primary at separations greater than
0.2$^{\prime\prime}$ ($\sim$5\,au) and companions with $\Delta I < 5$\,mag at separations
greater than 1$^{\prime\prime}$ ($\sim$26\,au).

The photometric data are listed in Table~\ref{tab:tab3}. The $I$-band
magnitudes of the distant BC pair were taken from \cite{2003AJ....125.3302G}
who imaged individual components with Hubble Space Telescope WFPC2. The
integrated $JHK_s$ photometry of the BC component from VHS were
also decomposed into individual magnitudes using the flux ratios derived by
\cite{2002ApJ...567L..53C}. The $JHK_s$ band photometry of the primary,
also given in Table~\ref{tab:tab3}, is from the 2MASS catalogue.

\begin{table*}
\caption{Photometric data for the components of the multiple system HD 221356.}
\centering
\begin{tabular}{c c c c c c}
\hline
\hline
Component & $I$ (mag) & $Y$ (mag) & $J$ (mag) & $H$ (mag) & $K\rm{_s}$ (mag) \\
\hline
 A &  5.953$\pm$0.018 & ... & 5.488$\pm$0.019 & 5.264$\pm$0.038 & 5.150$\pm$0.017  \\
BC & 15.536$\pm$0.027 & 13.695$\pm$0.025 & 12.852$\pm$0.010 &12.261$\pm$0.026 &11.946$\pm$0.026  \\
 B & 15.568$\pm$0.018 & ... & 12.933$\pm$0.011 & 12.353$\pm$0.026 & 12.055$\pm$0.028  \\
 C & 19.363$\pm$0.039 & ... & 15.713$\pm$0.042 & 14.993$\pm$0.061 & 14.495$\pm$0.114  \\
 D & 16.70$\pm$0.10 & 14.933$\pm$0.059 & 13.763$\pm$0.038 & 13.209$\pm$0.026 & 12.755$\pm$0.025 \\
\hline
\multicolumn{6}{p{.665\textwidth}}{\footnotesize Notes: $I$-band mag of A is from the FastCam, $JHK_s$ mag are from 2MASS.
$I$ mag of the BC is from \protect\cite{2003AJ....125.3302G}, $JHK_s$ mags are from VHS and were decomposed using the flux
ratios derived by \protect\cite{2002ApJ...567L..53C}. $I$ mag of D is from IAC80 measurement, $JHK_s$ mags are from our photometry on VHS images.}
\end{tabular}
\label{tab:tab3}
\end{table*}

\begin{figure*}
\centering
\includegraphics[scale=0.77]{fig3.pdf}
\caption{Near-infrared $H$- and $K_s$- band spectra of the new companion compared
with the M8 standard LP 412-31 \citep{1995AJ....109..797K}, observed with the same instrumentation, and the L1 (2MASS\,J14392836+1929149)
and L3 templates (2MASS\,J15065441+1321060) taken from the IRTF library \citep{2005ApJ...623.1115C}.
Spectra were normalised at 1.7 microns and offsets have been added for clarity. The grey area
indicates the region of high telluric absorption. The most prominent molecular and
atomic features are indicated.}
\label{fig:3}
\end{figure*}

\subsection{Near-IR spectroscopy}
We also obtained near-infrared spectroscopy of the HD\,221356 AD system and
the M8V spectroscopic standard LP 412-31 \citep{1995AJ....109..797K} using LIRIS, with the
$HK$ grism and the 1K$\times$1K Hawaii detector at the 4.2 m William Herschel
Telescope (WHT) on December 30 2011. This instrumental configuration provides
a nominal dispersion of 9.7\,\AA\,pix$^{-1}$ and a wavelength coverage of
1.4--2.4~$\mu$m. A slit width of 0.75" was used rotated to the direction
along the AD system and the final resolution of the spectrum was 26\,\AA (R\,$\sim$\,600).
Total integration time was 2240\,s, divided into individual
exposures of 160s. A nodding pattern of two positions (AB) was used to subtract
the sky background. Weather conditions were photometric and the average seeing
was 0.9". Data were dark corrected, sky subtracted, aligned and combined at
each nodding position.
Flat-field correction using a tungsten lamp was not applied due to a spurious feature in the $K$-band tungsten spectrum.
After subtracting the contribution of the primary wings, spectra of HD\,221356\,D
were optimally extracted using the \textsc{apall} routine and wavelength calibrated using
ArXe arc lines. We finally combined the spectra at both AB positions and
corrected for telluric absorption lines by dividing them by the A3V star
HR8840, observed at a similar airmass, and multiplying by a blackbody of the
corresponding effective temperature of 8500\,K. Spectroscopic data of
LP 412-31 were reduced and analized in a similar way to HD\,221356\,D,
but telluric correction was done using the A3V star HR\,1036, also observed at
a similar airmass. Spectra of HD\,221356\,D and LP 412-31 in comparison
with other standard objects and the main spectroscopic features are shown in
Figure~\ref{fig:3}.

HD\,221356\,D displays stronger water vapor absorption bands than the
standard M8 dwarf observed with the same instrumental setup and sky
conditions and reduced in the same manner as our target. This implies that
HD\,221356\,D has a cooler spectral type, quite likely within the L domain.
Aimed at deriving the spectral type of HD 221356D, we have
compared its LIRIS spectrum with data extracted from the IRTF libraries
\citep{2005ApJ...623.1115C}, see Fig. 3.
We note that these data correspond mostly to dwarfs with solar composition.
The overall $HK$ spectral energy distribution (SED) of HD\,221356D is better
reproduced by a spectral type of L0-L1. However we determine
a spectral type of L1--L3 if we consider the different water indexes at
$\sim$ 1.7 or 2.0 ~$\mu$m \citep{2001ApJ...552L.147T, 2007ApJ...657..511A, 2004ApJ...610.1045S}.
The differences can be explained if the object is slightly metal poor (consistent with the primary),
because the $K$-band flux should be reduced by the $H_2$ collisional induced absorption.
The strength of the sodium feature at 2.2 microns in object D compares better with solar metallicity
M9 dwarf than with the early L-type templates, however we prefer to rely on the general SED,
which is better fit by L1 templates, than on a single feature that can be uncertain bearing in mind the
poorly understood effects of low metallicity.

We also note that the water vapour absorption band at $\sim$ 1.5 ~$\mu$\,m, in the blue part of the
$H$-band is more intense in this object than in early L dwarfs and more resembles mid/late L dwarfs.
This is not an instrumental effect since it does
not appear in the M8.5 object, but it may be due to a larger contamination by the primary at blue
wavelengths than at red wavelengths. We therefore cannot, with confidence,
assign this feature as unusual for this object. In summary, we
estimate that HD\,221356\,D  is a slightly blue early L dwarf with a spectral type of L0--L2.

\begin{figure*}
\centering
\includegraphics[scale=0.8]{fig4a.pdf}
\includegraphics[scale=0.825, trim=0pt 10pt 0pt 0pt]{fig4b.pdf}
\caption{Left panel: $M_I$, $I-J$ colour-magnitude diagram of the HD\,221356 system.
The positions of the four components are marked with points and labelled with
the corresponding letters. Pleiades low-mass stars and brown dwarfs (crosses) from \protect\cite{2010AA...519A..93B},
M and L dwarfs from \protect\cite{2006PASP..118..659L} with available parallaxes
and M dwarfs (dots) from \protect\cite{2000ApJ...535..965L} are also plotted.
The 5 Gyr isochrones of the NextGen models for solar and
low metallicity stars \citep{1998A&A...337..403B}, represented by a dotted and dashed line,
respectively, and of the DUSTY model for solar metallicity \citep{2000ApJ...542..464C},
shown as a solid line, are also included. The masses in solar mass units from the DUSTY
model are also indicated and marked with open circles in the corresponding isochrone.
Right panel: $M_J$, $J-K$ of BCD components of the HD\,221356 system. We have also added
L and T dwarfs from \protect\cite{2004AJ....127.2948V} and \protect\cite{2011ApJS..197...19K} and
12 subdwarfs with measured parallaxes from \protect\cite{2012ApJ...752...56F}, depicted by squares.
The mean L and T near-infrared photometric sequence from \protect\cite{2004AJ....127.2948V}
is represented by a dashed line.}
\label{fig:4}
\end{figure*}
\section{Physical properties of HD 221356D}
Assuming a distance of 26.12$\pm$0.37 pc to the system, we calculated the absolute
magnitudes of individual components and constructed the $\left(M_I,~I-J\right)$
and $\left(M_J,~J-K_s\right)$ colour-magnitude diagrams, shown in Fig. \ref{fig:4}.
All four objects clearly follow a well defined photometric sequence,
with the new companion located between the known B and C components. All the
colours and magnitudes of the new object are in good agreement with its physical
membership to the HD 221356 system. To better illustrate the position of the
three low mass companions, we added in both panels of Fig. \ref{fig:4} the Pleiades
low-mass stars and brown dwarfs (at $\rm{d=120~pc}$) from \cite{2010AA...519A..93B},
M and L dwarfs from \cite{2006PASP..118..659L} with available parallaxes and
the field M dwarfs from \cite{2000ApJ...535..965L}. The Pleiades cluster offers a
homogeneous collection of objects with similar age and metallicity.
The right panel also includes L and T dwarfs from \cite{2004AJ....127.2948V},
\cite{2011ApJS..197...19K}, and 12 subdwarfs with spectral types
from M7 to L7 with measured parallaxes \citep{2012ApJ...752...56F}.
In the $\left(M_I,~I-J\right)$ diagram we have also plotted the 5 Gyr isochrones
of the NextGen models for solar and low metallicity $(\rm{[M/H]}=-0.5)$
\citep{1998A&A...337..403B} and of the DUSTY model for solar metallicity
\citep{2000ApJ...542..464C}.

We have depicted a veriety of objects with similar spectral types, that span from
young ages ($\sim$120 Myr) and solar composition, to old stars with solar metallicities
and field subdwarfs with metal-poor abundances ([Fe/H]$\sim-0.5$, \citealp{2007ApJ...669.1235L}).
The low-mass components of the HD\,221356 system have slightly bluer colours than typical
field stars, placing them on the blue edge of the photometric sequence defined
by late M and L dwarfs.
However their colours are redder than known field subdwarfs, implying an intermediate
metallicity, being thus in good agreement with the metallicity determination of the primary.
 The $J-K_s$ and $I-J$ colours of the new companion,
which are $1.01\pm0.06$ and $2.94\pm0.14$ mag respectively, correspond to the
typical values of M7--M9 field dwarfs\citep{2002ApJ...564..452L, 1994AJ....107..333K}.
These photometric colours suggest an earlier spectral type than that determined
in the spectroscopic analysis. We attribute this difference to the sub-solar
metallicity of the system, which affects the spectral energy distribution.
In particular, the flux suppression in the $K_s$-band was already
recognized in a group of L dwarfs, as a low metallicity feature
\citep{2010ApJS..190..100K, 2011ASPC..448..531W}. Because of that, the spectral
classification of the wide binary components of the HD\,221356 system, which is based on
photometry, may be uncertain.
\begin{table}
\caption{Physical properties of low mass companions in HD 221356 system.}
\centering
\begin{tabular}{c c c c c c c}
\hline
\hline
Comp. & SpType & $\log~L/L_{\odot}$ &  $\rm{T_{eff}}$ (K) & $\rm{M/M_{\odot}}$ \\
\hline
 B & $\rm M8.0\pm1.5$ & $-3.28\pm0.09$ & 2351$-$2717 & $0.090\pm0.008$ \\
 C & $\rm L3.0\pm1.5$ & $-4.35\pm0.09$ & 1756$-$2141 & $0.065^{+0.007}_{-0.016}$ \\
 D & $\rm L1.0\pm1.0$ & $-3.61\pm0.10$ & 2079$-$2304 & $0.079\pm0.006$ \\
\hline
\end{tabular}
\label{tab:tab4}
\end{table}

The luminosities of the B, C and D components were derived from their
$JHK_s$-band magnitudes, using the trigonometric distance of the primary,
bolometric corrections from \cite{2004AJ....127.3516G} and spectral
type-colour relations from \cite{2004AJ....127.2948V}. The effective
temperature ranges were calculated adopting the temperature scale for
high-gravity field dwarfs given by \cite{2004AJ....127.3516G}, assuming
a spectral type range of L0--L2 for the new companion and M6.5--M9.5,
L1.5--L4.5 for B and C, respectively. The resulting values are given in
Table \ref{tab:tab4}. The masses of the new companion and B and C
components were estimated using the DUSTY model from the Lyon group
\citep{2000ApJ...542..464C}, which is only available for solar metallicity.
We adopted a wide range of ages, using the 1, 5 and 10 Gyr isochrones.
Masses were derived from their luminosities by interpolating the
mass-luminosity relations given in the models.
The differences in mass calculated for different ages are lower than the errors in mass
determination resulting from the uncertainties of luminosities. For the
new companion we finally adopted a mass of $0.079\pm0.006 M_{\odot}$, which
is the average value obtained using $JHK_s$-band magnitudes,
a spectral type range of L0$-$L2, and an age range of 1$-$10 Gyr. To take account
of the differences in models for low metallicity stars, we have checked the mass
of object D using the 5\,Gyr NextGen model for $[M/H]=-0.5$ \citep{1998A&A...337..403B}. We obtained
a mass of $0.083\pm0.002$, which is slightly larger, but still
within the uncertainties of that determined for solar metallicity models.
\section{A benchmark sub-solar metallicity multiple system}
The age of the primary star in the HD\,221356 multiple system was estimated by
\cite{2005ApJS..159..141V} to be 2.5$-$7.9
Gyrs based on isochrone analysis. The lithium abundance of the primary star
($\log n(Li)=2.5$ in the usual scale of $\log n(H)=12$) is typical of late F-type
stars in clusters with ages in the 2$-$8 Gyr range \citep{2005A&A...442..615S}.
The chromospheric activity of the star is also typical
of a moderately old main sequence star \citep{2005ApJS..159..141V}.

HD 221356 is a slightly metal-poor stellar system whose components have masses just
above and below the hydrogen burning limit. For sub-solar metallicities the
stellar--brown dwarf borderline is expected to be shifted to higher masses, e.g.
to $\sim0.079~M_{\odot}$ at $[M/H]=-0.5$ ($0.072~M_{\odot}$ at $[M/H]=0$)
\citep{1998A&A...337..403B}.
In such a coeval, old system it becomes particularly interesting to investigate
the lithium abundances of the very low mass components. While the M8 star (component B)
and the L0$-$L2 (component D) should have fully burnt their original lithium, component
C with 0.065 solar mass may have preserved some amount of the initial lithium content.
Theoretical models for solar metallicity predict full lithium depletion for such a mass,
however this may not be the case for sub-solar metallicity, since models also predict
a less efficient depletion at low metallicities \citep{1997A&A...327.1039C}.
Observations of the Li abundance in the three low-mass components of this system will constrain
both the evolutionary models and the age of the system.

Multi epoch measurements of the system will allow to detect the orbital motion
of companion D. Although the estimated orbital period (assuming
circular orbit) is of the order of 5500 years, the relative change of the
position would be up to $\sim 14~\rm mas/yr$, which is measurable using modern
high spatial resolution imaging (e.g., adaptive optics, lucky imaging). The
expected semi-amplitude of radial velocity variation of the primary induced
by the presence of companion D will be of the order of 130 m/s (for 90$^{\rm o}$
inclination), however the orbital period is too long to allow a full determination
of the three-dimensional orbit. Maximum annual variations of roughly 0.15 m/s
are expected, which may be explored with the new generation of ultra high
precision spectrographs \citep{2012Natur.485..611W}.
\section{Conclusions}
Using the VHS and 2MASS surveys we have identified a new very low-mass
companion (HD 221356 D) in the slightly metal-poor HD 221356 system, which thus becomes a
quadruple. The new object is located
at a projected distance of $\sim$ 312 AU from the F8 primary.
The four components of the system follow a well defined photometric sequence.
From near-infrared spectroscopy
we determined L0--L2 spectral type for the D companion. Based on
evolutionary models its mass is estimated at $0.079\pm0.006~\rm M_{\odot}$,
and its effective temperature is in the range from 2100 to 2300 K.
The $J-K_s$ and $I-J$ colours of the low-mass components are slightly bluer
than field counterparts of the same spectral type.
We interpret this as a result of the low metallicity of the system, which may
become a reference for the spectral classification of metal poor M and L-type field objects.
Since the distance and metallicity of the HD 221356 system are well known,
the detailed study of its ultracool
companions, which are located above and below the frontier
between stars and brown dwarfs, can provide valuable constrains on evolutionary models and,
in particular, shed light on the properties of objects on the transition from stellar to
substellar regime.

{\bf Acknowledgements.}
Based on observations obtained as part of the VISTA Hemisphere Survey, ESO Program; 179.A-2010 (PI: McMahon).
The VISTA Data Flow System pipeline processing and science archive are described in \cite{2004SPIE.5493..411I}
and \cite{2009MNRAS.399.1730C}.
This article is based on observations made with the TCS and IAC80 telescope operated on
the island of Tenerife by the IAC in the Spanish Observatorio del Teide,
and with the William Herschel Telescope operated on the island of
La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos
of the Instituto de Astrof\'{i}sica de Canarias. This publication makes use of data
products from the Two Micron All Sky Survey, which
is a joint project of the University of Massachusetts and the Infrared Processing
and Analysis Center/California Institute of Technology, funded by the National
Aeronautics and Space Administration and the National Science Foundation. This
research has benefitted from the SpeX Prism Spectral Libraries, maintained by
Adam Burgasser at http://pono.ucsd.edu/~adam/browndwarfs/spexprism.
This research has been supported by the Spanish Ministry
of Economics and Competitiveness under the projects
AYA2010-21308-C3-02, AYA2010-21308-C03-03 and AYA2010-20535.
V.J.S.B. and N. L. are partially supported by the Spanish Ram\'{o}n y Cajal program.
\bibliographystyle{mn2e}