Rapid Secular Radial Mass Accretion in N-Body Simulated Galaxy Disks
$^{1}$Department of Physics and Astronomy
Fairfax
VA 22030
USA
Recent studies have indicated that large mass radial inflow rates are possible in observed galaxies with strong density wave patterns. Yet, numerical simulations have generally failed to account for such high rates. Here it is shown that the reason for the discrepancy is the treatment of ``softening'', the artificial parameter inserted by numerical simulators into the formula for gravitational potential, to control the magnitude of relaxation in simulations with small numbers of particles compared to a real galaxy. Excess softening reduces the collective effects underlying the significant secular evolution inferred for physical galaxies. Less softening, coupled with an increase in number of particles, allow the N-body simulations to reveal significant morphological transformation of galaxies over a Hubble time.
ms.tex
\documentclass[useAMS,usenatbib]{mn2e}
\def\mic{$\mu$m\ }
\def\br{\hfill\break}
\def\etal{{\it et al.}}
\def\arcs{\rlap{.}$^{\prime\prime}$}
\def\arc{$^{\prime\prime}$}

\title[Rapid Secular Evolution in N-Body Galaxy Disks] 
{Rapid Secular Radial Mass Accretion in N-Body Simulated Galaxy Disks}
\author[Xiaolei Zhang]{X. Zhang$^{1}$\thanks{E-mail: xzhang5@gmu.edu} 
\\
$^{1}$Department of Physics and Astronomy, George Mason University, 
4400 University Drive, Fairfax, VA 22030, USA}

\begin{document}

\date{Accepted . Received ; in original form }

\maketitle

\label{firstpage}

\begin{abstract}
Recent studies have indicated that large mass radial inflow rates are possible 
in observed galaxies with strong density wave patterns. Yet, numerical 
simulations have generally failed to account for such high rates. Here 
it is shown that the reason for the discrepancy is the treatment of 
``softening'', the artificial parameter inserted by numerical simulators
into the formula for gravitational potential, to control the magnitude 
of relaxation in simulations with small numbers of particles 
compared to a real galaxy. Excess softening reduces the collective effects 
underlying the significant secular evolution inferred for physical galaxies. 
Less softening, coupled with an increase in number of particles, 
allow the N-body simulations to reveal significant morphological
transformation of galaxies over a Hubble time.
\end{abstract}

\begin{keywords}
galaxies: evolution; galaxies: structure \end{keywords}

\section{INTRODUCTION}

The prospect of significant transformation of the 
Hubble type of a galaxy during its lifetime
has important implications on our understanding of the
origins of the Hubble sequence and on broader issues
in cosmology.  Despite the coming of age of this field
of research, there is as of now no consensus on how
important this process is, and through what dynamical
mechanism this transformation is mainly achieved (see
the review of Kormendy \& Kennicutt 2004).

Zhang (1996, 1998, 1999) proposed and demonstrated a collective
dissipation process in galaxies that could serve as the underlying
dynamical mechanism for secular evolution of galaxies along
the entire Hubble sequence, by incorporating the stellar component
in the mass inflow process rather than just the interstellar medium as
proposed in the pseudo-bulge formation scenario (Kormendy 1979).  
Despite the confirmation of the theoretically
derived mass flow rates in the accompanying N-body simulations, which
established the viability of the analytical approach, the simulated
mass flow rates were small, far from being adequate to lead to 
significant morphological transformation of a galaxy within
a Hubble time.  This small mass flow rate was attributed
to the small amplitudes of the spiral patterns formed in these
simulations, with the state-of-the-art in simulating spontaneously
formed modes far from being adequate to produce the extremely nonlinear
density wave modes observed in physical galaxies. 

In order to circumvent this difficulty, Zhang \& Buta (2007; 2012)
used near-and-mid-infrared images of galaxies, coupled
with the analytical rate equations derived in Zhang 
(1996, 1998, 1999) to calculate the torques and
mass flow rates in these physical galaxies, and found that
mass flow rates from a few $M_{\sun}$ per year to over
one hundred $M_{\sun}$ per year were obtained, depending
on the Hubble types of galaxies and their interacting environment,
which together determined the density wave amplitudes and
pitch angles in these galaxies.
These mass flow rates for real galaxies far exceed those obtained from the
past N-body simulations, partly as a result of the fact
that the mass flow rate is proportional to the effective
wave amplitude squared (Zhang 1998), and the amplitudes in physical
galaxies are oftentimes a factor of 3-8 times higher
than obtained in N-body simulations, implying a difference
in evolution rate of a factor of 10 - 100 between physical
and previously-simulated galaxies.

Despite the confirmation of the importance of secular
evolution process in physical galaxies in works such as
Zhang \& Buta (2007,2012), one could not help but
wonder what exactly were the factors that have prohibited
the simulated disk galaxies from obtaining the level of mass
flow rates and the magnitude of wave amplitudes in physical
galaxies.  It is to this question that the result of the
current Letter partly addresses.  Incidentally, the
main discovery of this work was made in a totally
serendipitous fashion -- i.e. the result was stumbled upon
while the author was searching in a different direction
for the cures of the low mass flow rates in N-body simulations\footnote{The
author used to think that the low mass flow rates were mainly
due to the two-dimensional (2D) nature of her simulations, and once the bulge
and halo were allowed to be active, the mass flow rates
will become higher.  This turned out to be true,
but only moderately true: i.e., having about 15\% active
bulge and halo in 2D simulations increases the mass inflow rate by about 10\%,
but having still larger active bulge and halo fraction leads rather
to outflows of mass over the main parts of the disk (likely
due to the lowering of Q and the emergence of transient
local instabilities), contrary
to what is needed for bulge building. 
Note that the above conclusion
about the fraction of allowable active bulge/halo mass applies only in
the 2D cases where these tests had been made.  In three
dimensional (3D) studies the entire bulge and halo mass of course
can be made active with proper pressure and rotational support.}).
The lesson learned is that one has to be cautious
in taking the results of N-body simulations
literally without questioning sacred ``rules of thumb'' and
accepted wisdom, especially in circumstances where new physical
processes are encountered -- in our situation, that is when the
collective effects become important.

\section{ROLES OF SOFTENING: THEORETICAL AND NUMERICAL}

To come to the core result of this Letter immediately: the main culprit 
to the low mass accretion rate in N-body simulations of disk galaxies 
turned out to be the improper choice of softening parameter used to 
control relaxation effects in the simulation of these nominally 
collisionless systems.

The softening parameter is generally defined as
follows:  The gravitational potential in an N-body system is
calculated as
\begin{equation}
\Phi= {{-Gm} \over {\sqrt{(r^2 + {a_{soft}}^2)}}}
\end{equation}
where $a_{soft}$  is the so-called softening parameter.  This parameter is
introduced into the N-body simulations due to the fact
that the numerical setups usually have significant less
number of particles than the physical systems they try to emulate.
A finite softening parameter reduces the amount of unrealistic close encounters
which would be more often for a small N system, and thus reduces
the level of artificial relaxation in these small N systems.
The hope is that a collisionless configuration would then result
from a well-matched N and $a_{soft}$.
In particle-mesh based N-body approaches, the grid size serves as an additional
softening parameter that acts in consort with $a_{soft}$ in controlling 
relaxation (see Sellwood 1987 and the references therein).

However, with the realization that galaxies which possess large
scale density wave modes are governed by collective dissipation
processes (Zhang 1996), one soon faces the fact that it is impossible
to completely avoid collision-like behavior,
no matter how large the particle numbers are used: After all, the
collective effects in these systems {\em depend} on the near-collision
or scattering of particles in the global instabilities to set up
long-range correlations to achieve self-organization
and to induce secular evolution.  In some sense, these galactic systems
are in a {\em forced relaxation configuration}, with the forcing
accomplished by the density wave collisionless shocks, which force
the effective Q in the spiral arms to be less than one and lead
to inter-particle interactions and correlations.  Thus, in hindsight,
it is also obvious that excess softening reduces the very inter-particle
interaction that is the backbone support of collective effects.
In essence, overly-softened gravity is {\em modified gravity}, and thus
is modified Newtonian interaction (not to be confused with the
Modified Newtonian Dynamics proposed to account for the flat rotation curves). 
The system one simulates with softened gravity is no longer exactly
the same system governed by Newton's gravitational law, but rather 
a more sluggishly-interacting system.  This might not be a serious concern 
if one's interest is in obtaining a modal morphology that mimics the observed
galaxy morphology (as shown in the work of Donner \& Thomasson 1994, 
hereafter DT94, which is the forerunner of long-lasting N-body spiral modes), 
but as we will show here, it is detrimental
to determining the realistic secular mass flow rates that
are relevant to physical galaxies.

In Fig.~\ref{fg:Fig1}, we show the result of a set of test runs
of N-body simulations using the basic state first explored in
DT94, and subsequently used for the study of collective effects 
and secular evolution in Zhang (1996,1998,1999).
The simulations are done in a 2D configuration, using
a particle-mesh approach on a polar grid, with an
active disk which has a modified exponential
surface density profile, as well as a rigid bulge and a rigid halo.
The grid size is 220 by 256 in the radial and azimuthal direction, respectively.
The mass in the active disk, the inert halo and 
inert bulge are 0.4, 0.5, 0.1, respectively, in the normalized unit,
similar to that used in Zhang (1998). The scaling of the time unit
is such that 1256 time steps represent one rotation period at r=20.
Both the particle numbers and the time resolution used are higher
than that in Zhang (1998).  These grid and time step resolutions
are used for the rest of the simulations presented in this
Letter as well, only particle numbers and the softening
paramters are changed and will be noted in each case.

The curves in Fig.1 are for the evolution of enclosed mass within the
central $r=3.5$ region (in the scale of this disk which has
corotation radius around 20-25 in normalized unit, this corresponds to
the central bulge region).  The three curves, from bottom to
top, correspond to softening parameter value of 1.5 (as used
in DT94 and Zhang 1998), 0.75, and 0.25 in the unit of the
smallest grid size of 1.  As we can see, the choice of softening
of 1.5 produced a barely noticeable mass inflow (to really see
the small amount of mass inflow for this choice of softening, 
one is referred to Figure 10 of Zhang 1998 which zoomed
in on the vertical scale, even though that was for the enclosed mass within
central $r=13.5$).  For that particular choice of softening,
which is of standard usage in the N-body simulation community
as being adequate in suppressing unwanted relaxation 
effects as well as in matching the grid softening, the central
mass growth of only 6.5\% was observed over about 25 pattern
rotation periods (Zhang 1998), which is obviously quite inadequate in
transforming the Hubble type of a galaxy.  
Gradually reducing softening is seen to
systematically lift up the mass inflow rate.  For the range
of softening tried by the current author (1.5 - 0.1), the
trend of increasing mass inflow with reduced softening
continues.  It is possible that realistic galaxies, which
have zero softening in the form of gravitational law but
do have much larger numbers of particles, as well as the
finite thickness of the disk to keep relaxation in check,
the mass flow rates much exceed those we have seen in these
recent simulations are possible -- as we indeed observed in
the mass flow rates calculated for the sample in Zhang \& Buta
(2007, 2012).

\begin{figure}
\vspace{150pt}
\centerline{
\special{psfile=Mt_asoft.ps
         hoffset=-100 voffset=0 hscale=25 vscale=25}
}
\caption{Evolution of enclosed disk mass within central r=3.5 radius
of a set of 2D N-body simulations
with different softening parameters. Solid: $a_{soft}=1.5$,
Dashed: $a_{soft}=0.75$, Dotted:i $a_{soft}=0.25$. 1 million
particles are used to represent the active disk.}
\label{fg:Fig1}
\end{figure}

\section{A NEW SET OF N-BODY SIMULATIONS WITH SMALL SOFTENING:
ACCELERATED BULGE BUILDING}

The price one pays for reduced softening in small-N simulations
is the increase in relaxation rate.  To compensate, one needs
to use a correspondingly increased number of particles, for
otherwise the collision-induced heating will swamp the
density wave activity and reduce the wave amplitude.
How much the particle number needs to be increased for a given
choice of softening can be found empirically (recalling
the fact that the traditional rules of thumb in cases
with collective effects may not always apply, since we do not
want to kill the right kind of collisional relaxation due to
collective effects).

In what follows, we show a set of simulations with changing particle numbers
while holding the softening parameter $a_{soft}=0.1$.
Fig.~\ref{fg:Fig2} shows the morphological evolution of a spontaneously
formed spiral-bar mode for one of these runs using 20 million particles.
Note that even though the duration of run for this mode
is similar to that used in Zhang (1998) (note due to the different
temporal resolution the 400 step here corresponds to 100 step in
the older simulation), the modal morphology
here appears to have changed much more that that seen in Zhang (1998),
i.e. the pattern in the older simulation retained the spiral morphology
throughout, whereas the pattern in the current simulation evolved
from a spiral-like morphology to a more bar-like morphology.
This is natural to expect, since the much-higher radial mass
accretion rate in the current simulation
resulted in a significant increase in central
mass concentration at the end of the run, 
this further means that the basic state is change so much
that the mode it supports needs also to change
correspondingly from that of spiral to that of bar with possibly
nested inner spiral (note that a heavy bar is
the natural modal shape for a centrally-heavy disk - the observed
early-types disks certainly contain many massive bars).
The pattern speed, however, appears to be little affected
by the change in morphology, as can be seen later in the
clean concentration of the power spectrum around a single
pattern speed throughout the duration of the run.

\begin{figure}
\vspace{250pt}
\special{psfile=n0.ps
          hoffset=0 voffset=130 hscale=25 vscale=25}
\special{psfile=n8000.ps
          hoffset=120 voffset=130 hscale=25 vscale=25}
\special{psfile=n16000.ps
          hoffset=0 voffset=0 hscale=25 vscale=25}
\special{psfile=n24000.ps
          hoffset=120 voffset=0 hscale=25 vscale=25}
\caption{N-body morphology of a spiral/bar mode at different
time steps. The rotation period at r=20 is about 1256 steps.
20 million particles are used to represent
the active disk, the softening parameter is $a_{soft}=0.1$.  
About 10,000 randomly-selected particles are plotted.}
\label{fg:Fig2}
\end{figure}

In Fig.~\ref{fg:Fig3}, we present the enclosed mass within the central r=3.5
radius, with the different runs having a constant $a_{soft}=0.1$ but varying the
number of particle used.  For the smaller particle number runs,
the rapid mass inflow is seen to saturate at an earlier time step,
which corresponds to the time when spiral activity is damped by
heating due to the insufficient number of particles for this
extremely small softening,
especially in the outer disk region where the surface density
is low and the particle
numbers are small.  For large particle runs, the rapid mass
inflow is seen to remain at a constant rate all the way until
the end of the run, which corresponds to 25 rotations.

Furthermore, we note that the mass inflow rate observed here
corresponds to about 60\% increase in enclosed mass within r=3.5
(the bulge region) over 10
rotation periods, or roughly 1/5 of a Hubble type.  This
level of mass accretion is more than sufficient to transform
the Hubble type of a galaxy by several stages in a Hubble time.
Plus, the mass inflow rate conceivably could be even higher
for even smaller choices of $a_{soft}$.  The natural limit
apparently is set only when one has the actually physical galaxy
in its 3D representation with the full inter-particle interactions
taken into account (i.e. we know what the actual density wave
amplitudes in galaxies are because we observed them).

\begin{figure}
\vspace{150pt}
\centerline{
\special{psfile=Mt4.ps
         hoffset=-100 voffset=-10 hscale=30 vscale=30}
}
\caption{Enclosed mass within central r=3.5 of the N-body disk
for runs with different numbers of particles and $a_soft={0.1}$.
Solid: 1 million particles. Dashed: 10 million particles.
Dotted: 20 million particles. Dash-Dotted: 40 million particles.}
\label{fg:Fig3}
\end{figure}

In Fig.~\ref{fg:Fig4},  the disk surface density at the beginning of the run,
and at the end of the four runs with different numbers of
particles (as used in Fig~\ref{fg:Fig3}, but with different line
style due to the need to represent the ``before'' surface density) 
are presented.  It can be seen that
the large inflow of the large-particle-number runs indeed builds
up a more substantial bulge.  However, the 20 million
particle run in fact produced more substantial bulge-building
than the 40 million particle run.  This is consistent
with the trend observed in Fig.~\ref{fg:Fig3}, where it is seen that the
40 million particle run has somewhat reduced mass inflow rate
near the end of the run.  This may be related to the emergence
of different kinds of local instabilities in the 
different particle runs (a small knot of globular-cluster-like 
morphology can be seen in the outer disk of the
second and third frame of Fig.~\ref{fg:Fig2}).  This kind of
cluster-forming activity in association with density wave modal
activity may be a reflection of the actual process happening in
physical galaxies.  By step 32000 (not shown in Fig~\ref{fg:Fig2}),
the 20 million particle run has developed two prominent dumb-bell
like structures at the end of the massive bar which might
account for the increased mass inflow for this run.

\begin{figure}
\vspace{160pt}
\centerline{
\special{psfile=BS5.ps
         hoffset=-100 voffset=-0 hscale=30 vscale=30}
}
\caption{Basic state surface densities at the beginning (solid), and the
end of four N-Body runs with different number of Particles.
Dashed: end surface density for 1 million particles.
Dotted: end surface density for 10 million particles. 
Dash-and-Single-Dotted: end surface density for 20 million particles.
Dash-and-Triple-Dotted: end surface density for 40 million particles.
$a_{soft}=0.1$}
\label{fg:Fig4}
\end{figure}

In Fig.~\ref{fg:Fig5}, the m=2 power spectrum together with the $\Omega$,
$\Omega \pm \kappa/2$ curves are presented for two runs with
two different numbers of particles.  These contour
plots clearly show a single dominant mode throughout the
time interval of the runs.  These are the cleanest power
spectra (albeit obtained using potential m=2 components
rather than the more conventional surface density) 
yet observed for this particular
mode, compared to those runs with larger softening and
smaller number of particles (note that due to higher
temporal resolution by a factor of 4 compared to that
used in Zhang 1998, the power spectrum should be multiplied by
by the same factor of 4 to compare with older results.
After taking this into account, the pattern speed still
appears to be slightly higher than before -- even the
40 million run appears to have a slightly higher pattern
speed than the 20 million particle run).

\begin{figure}
\vspace{180pt}
\centerline{
\special{psfile=power_pot_20m_0_32768.ps
         hoffset=-100 voffset=100 hscale=30 vscale=20}
\special{psfile=power_40m_0_32768.ps
         hoffset=-100 voffset=0 hscale=30 vscale=20}
}
\caption{Power spectra from m=2 Potential at different radii during
time interval step=0 and step=32768. Top: 20 million particle run.
Bottom: 40 million particle run. The solid line indicates
galactic rotation speed $\Omega$, the two dashed lines
are $\Omega \pm \kappa/2$.}
\label{fg:Fig5}
\end{figure}

In Fig.~\ref{fg:Fig6}, we plot the mass inflow (positive part of the curve)
and outflow (negative part of the curve) versus radius at
time step 8000.  This mass flow rate curve was calculated based
on the volume torque prescription given in Zhang (1996, 1998).
The major transition from inflow to outflow is
close to r=20 for the inner disk, which is close to
(but not exactly the same, due probably to the fast
evolving nature of the surface density) the
corotation radius predicted from pattern speed (Fig.~\ref{fg:Fig5}).
Note that the larger number-of-particle case the power
spectrum predicts a slightly smaller corotation radius
than the smaller number-of-particle run.
The level of peak accretion, $2-4 \times 10^{-7}$
per time step for the inner disk region, matches quite well
the average accretion rate shown
from the slope of the Fig.~\ref{fg:Fig3} 20 million particle curve,
which corresponds to $2 \times 10^{-7}$ per time step,
within the central  r=3.5.  It should be noted
that the small peaks and valleys of the mass flow curve in the inner disk
do change with time and in general the peaks move inward with time.

\begin{figure}
\vspace{150pt}
\centerline{
\special{psfile=accretion_8000.ps
         hoffset=-100 voffset=-0 hscale=30 vscale=30}
}
\caption{Mass inflow (positive) and outflow (negative) versus radius
at time step 8000}
\label{fg:Fig6}
\end{figure}

Thus this new set of simulations show that the broad modal
nature of the pattern as well as the predictions of the
collective dissipation formalism
persisted even with the drastic reduction in softening
length.  The dominant mode morphology does evolve a lot faster in
appearance, which is to be expected in a fast changing 
basic state.

Finally we note that with the drastically-reduced $a_{soft}$
in the equation for the N-body gravitational potential in
our new set of simulations, the
corresponding grid softening has not been much reduced
(the newer grid is only about a factor of 4 finer in linear
resolution than used in Zhang 1998, whereas the softening is
reduced by a factor of 15). In fact Zhang 1996 already
explored a grid twice as fine as used in DT94 (or in Zhang 1998) and
the mass inflow rate was not found to have changed much.
Thus it appears that particle softening is the main
inhibitor to obtaining higher mass inflow rates in N-body
disks containing collective modes, the grid softening
has only marginal effect in this regard.  This is reasonable as the main
contributor to collective effects is inter-particle interactions,
and the finite amount of grid softening perhaps also helped
to maintain the desirable finite thickness effect in these
2D simulations.

\section{DISCUSSION}

\subsection{Are the Increased Mass Flow Rates in
Less-Softened N-Body Runs Due to the Collective
Instabilities or Due to Few-Body Encounters?}

The answer is that it is due to collective instabilities associated
with the density wave modes.  The hint for 
this conclusion can be found again in Fig.~\ref{fg:Fig3}.
It is seen there that when not enough
number of particles are used for a given $a_{soft}$
(which generally {\em increases} the few-particle relaxation effect),
the mass inflow tapers off at an earlier time step.
Whereas for very large particle numbers (which
generally {\em decreases} the few-particle relaxation effect), both
the spiral/bar activity and the mass inflow continues
until the end of the run.  This correlation between the number
of simulation particles and the longevity of effective mass
inflow shows that the survivability of the density wave
modes is key to the continued mass inflow.  Few-particle
effect causes negligible mass inflow once the spiral/bar
activity ceases. 

\subsection{The Remaining Unrealistic Elements
in the Secular Bulge Building Simulations}

The simulations described in this Letter serve to highlight
an important, and so far neglected, aspect of N-body simulations
that has a huge impact on the ability of these simulations
to reproduce the observed level of mass inflow for bulge building
as inferred from the MIR observations (Zhang \& Buta 2007, 2012).
However, since the simulation is conducted in 2D, obvious
physical effects that are missed include the continued growth
of the disk through halo mass accretion,
the role of active disk and halo (though as discussed in the
footnote, we have a general idea of what this will do to the
enhancement of the mass inflow rate).
Furthermore, since this is a pure particle simulation, the
role of gas is ignored, though the star-gas two component
simulation conducted in Zhang (1998) shows that the general
features of the collective effects are preserved when gas is added, and
Zhang \& Buta (2012) showed that the overall effect of stellar
accretion is more important than gas accretion due to the larger stellar
surface density, for all galaxies except the very late types.

\section{CONCLUSIONS}

The low mass inflow rates found in the past N-body
simulations of disk galaxies had cast doubt on
the effectiveness of secular evolution as an
important dynamical process for the morphological
transformation of galaxies along the Hubble sequence.
In this work it is shown that such low numerical mass inflow
rates are chiefly the result of the artificial
``softening'' of gravity which is a common
practice in galaxy simulations to avoid rapid
relaxation due to the small number of particles
employed.  By decreasing the amount of softening
and simultaneously increasing the number of simulation
particles, realistic level of mass inflow comparable
to those inferred from observed galaxy images
through torque calculations are achieved in N-body
simulations of disk galaxies, which reinforce
the importance of secular evolution as an extremely relevant
process in transforming the morphologies of galaxies
during the past Hubble time.


\section*{ACKNOWLEDGMENTS}

The strong objections of J. Sellwood when serving as the
referee to a separate manuscript
coauthored by Zhang \& Buta, which used mid-infrared galaxy
images to calculate torques and mass flow rates, but did
not contain new N-body simulations to corroborate these
observationally-derived results, became the direct motivation for
the author to carry out a new set of N-body simulations
to explore possible means to improve the mass flow behavior of her
old codes.  She thanks both Sellwood and her collaborator on the
observational study R. Buta for serving as the catalyst
which made this breakthrough possible.

\section*{REFERENCES}


\noindent
Donner, K.J., \& Thomasson, K., 1994, A\&A., 290, 785 

\noindent
Kormendy, J., 1979, ApJ, 227, 714

\noindent
Kormendy, J., \& Kennicutt, R.C.Jr. 2004, ARA\&A, 42, 603

\noindent
Sellwood, J.A. 1987, ARA\&A, 25, 151

\noindent
Zhang, X. 1996, ApJ, 457, 125 

\noindent
Zhang, X. 1998, ApJ, 499, 93

\noindent
Zhang, X. 1999, 518, 613 

\noindent
Zhang, X. \& Buta, R. 2007, AJ, 133, 2584 

\noindent
Zhang, X. \& Buta, R. 2012, MNRAS submitted, arXiv:1203:5334 

\end{document}
