Adding a Suite of Chemical Abundances to the MACER Code for the Evolution of Massive Elliptical Galaxies
Zhaoming Gan Ena Choi Jeremiah P. Ostriker Luca Ciotti Silvia Pellegrini



        

^1Shanghai Astronomical Observatory, Chinese Academy of Sciences,
        80 Nandan Road, Shanghai 200030, China
        
^2Department of Astronomy, Columbia University,
        550 W, 120th Street, New York, NY 10027, USA
^3Department of Physics and Astronomy, University of Bologna,
        via Piero Gobetti 93/2, 40129 Bologna, Italy




We add a suite of chemical abundances to the  
(Massive AGN Controlled Ellipticals Resolved) 2D code, 
by solving 12 additional continuity equations for
H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe and Ni respectively 
with sources from AGB stars and supernovae of type Ia and II
with metal yields based on standard stellar physics. 
New stars, formed in Toomre unstable circumnuclear disks 
(of a size ≲150 parsec), are 
assumed to have a top-heavy initial mass function with  a power index of 1.65. 
The metal dilution effects due to cosmic accretion are also included. 
With the high resolution of few parsecs in central regions, 
resolved black hole accretion and 
AGN feedback,  we can track the metal enrichment, 
transportation and dilution throughout the modeled massive elliptical galaxy of velocity dispersion 
∼280 km/s. We retrieve the chemical composition of 
the BAL winds launched by the central AGN, 
synthesize the X-ray features of the hot ISM,
and find that 
(1) the simulated metallicity  in the BAL winds could be up to ∼ 8 Z_⊙, 
while that of the hot ISM in the host galaxy is ∼ 2.3 Z_⊙, matching well with 
SDSS observations of BLR gas;
(2) the X-ray emitting hot gas is metal enriched with a typical value ∼2.5 Z_⊙;
(3) the circumunuclear cold gas disk,
where the metals are condensed, further enriched and recycled,
plays a critical role in the metal enrichment;
(4) the black hole accretion rate Ṁ_ BH linearly correlates with
the star formation rate Ṁ_⋆^+ in the circumnuclear disk, 
i.e, Ṁ_⋆^+ ∼ 7.7 Ṁ_ BH, but lagged in time by roughly 10^6 years.









§ INTRODUCTION
 
We have developed the  (Massive AGN Controlled Ellipticals Resolved) code 
over some time as an instrument for exploring the evolution of massive elliptical galaxies 
at high spatial resolution including a relatively complete set of physical processes. 
A recent paper <cit.> outlines the details of how rotation, 
massive dark matter halos, infall of cosmological gas and other features were 
recently added and tested. With this additional physical infrastructure 
we discovered features expected in the normal evolution of elliptical galaxies 
that surprised us but for which there is ample observational evidence: 
specifically the formation of circumnuclear cold gaseous discs that are Toomre unstable 
and which form massive stars in episodic bursts that fuel the central black hole and 
explain the enigmatic “E+A” phenomenon. The code, with resolution of parsecs 
in the central region, resolves the fiducial Bondi radius and hence can treat black hole accretion 
and AGN (Active Galactic Nucleus) feedback 
in some detail including both radiative (UV & X-ray) and 
BAL (Broad Absorption Line) winds in a fashion that imitates the observed output of black holes 
in both their low and high output states.

In this introductory paper we describe how we can add 12 chemical species produced 
by standard stellar physics and injected into the interstellar medium (ISM) by 
asymptotic branch stars (AGBs), type I and type II supernovae (SNe I, SNe II). 
This will enable us to address a whole new range of questions. For example, 
“What is the chemical composition expected of the BAL winds?” or 
“What iron abundance should be seen in the hot X-ray emitting gas?”. 
These are questions that most cosmologically based galaxy evolution codes would have 
difficulty addressing because they typically do not have the spatial resolution 
needed to resolve the central regions where the winds are generated 
and a significant part of the new metals are produced.

Half of the observed elliptical galaxies are known to contain 
cold gaseous discs in their centers
(a.k.a circumnuclear disks, e.g., 
). 
As we demonstrated in <cit.>, 
the formation of circumnuclear disk is inevitable
due to the angular momentum barrier, as rotation is the general case. 
Since the length scale of an AGN is very small when compared to the galactic scale,


most of the infalling gas (due to radiative cooling) 
will settle in a flat, rotationally supported disk-like structure 
before it reaches the galaxy center (see also ).
Moreover, the gas will cool down further in the circumnuclear disk,
which makes it ideal for star formation.


As cold gas condenses and accumulates on the circumnuclear disk,
such a cold gaseous disk would become highly over-dense, and we found that 
some of the over-dense disk rings would be gravitationally unstable 
(a.k.a. the Toomre Instability; ).
Consequently, angular momentum transfer and star formation are allowed 
at the same time due to the Toomre instability
<cit.>. 
Therefore, one may expect a near coincidence of star formation and AGN activity 
— the former is an important metal enriching source 
(by SNe II; ), 
and via accretion onto the supermassive black hole from the circumnuclear disk 
and then the AGN wind feedback, 
some of the metal enriched gas will be recycled
in form of BAL winds — 
in this paper, we will also demonstrate the key role of the circumnuclear disks
in the chemical evolution of massive elliptical galaxies, i.e.,
where the metals are condensed, further enriched and recycled.
	 
In the next section 
of this paper we will provide the details of 
how we implement the chemical abundance inputs to the code. 
We outline briefly in <ref> the model improvements
we have made to facilitate the introduction of diverse metallicity components.
In <ref> we will present some results of computations, 
including the chemical distribution in the modeled galaxy, 
the radiative features of the metal rich gas, composition expected in the BAL winds. 
Conclusions are reserved for final sections.
		
		


§ CHEMICAL ABUNDANCES
 
In the  simulations, we solve the continuity equation of the ISM, 
together with the conservation laws of momentum and energy, 
including source terms due to the stellar evolution of
AGBs  (ρ̇_⋆), 
SNe Ia (ρ̇_ I) and
SNe II (ρ̇_ II) 
(), 
i.e., 

    ∂ρ/∂ t + ∇·(ρ v) + ∇·ṁ_ Q  
       = ρ̇_⋆ + ρ̇_ I + ρ̇_ II - ρ̇_⋆^+,

where  ρ and v are the mass density and fluid velocity, respectively.
∇·ṁ_ Q and ρ̇_⋆^+ are the mass advection 
due to the Toomre instability <cit.> 
and the mass sink term due to star formation, respectively.
We refer to our code paper ( and references therein) 
for the full description of the physics and equations we solve in the simulations. 

In this paper, we intend to track the chemical evolution of the metals 

by using 12 tracers 
X_i (i=1, 2, ..., 12; mass of each element per unit volume) 
for H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe and Ni, respectively 
(see also ,  
for large scale simulations using SPH).
We solve 12 additional continuity equations of the tracers, assuming the chemical species 
co-move once after they are injected into the ISM, i.e.,

    ∂ X_ i/∂ t + ∇·(X_ i v) + ∇·ṁ_ Q,i 
            = Ẋ_⋆,i + Ẋ_ I,i + Ẋ_ II,i - Ẋ_⋆,i^+,

where

    ṁ_ Q,i = (X_ i/ρ) ·ṁ_ Q,   Ẋ_⋆,i^+ = (X_ i/ρ) ·ρ̇^+_⋆,

and v is obtained by solving the hydrodynamical equations as ususal.

In Equation <ref>, the passive stellar evolution, i.e., AGBs (Ẋ_⋆,i), 
SNe Ia (Ẋ_ I,i) and SNe II (Ẋ_ II,i),  serves as metal-enriching sources, 
while with different metal compositions (see Table <ref>).  
The advection terms ∇·(X_ i v) and ∇·ṁ_ Q,i 
describe the transport and the mixing of ISM with different metal abundances.
Star formation (Ẋ_⋆,i^+) is treated as a sink of local metals, 
but not to change the local abundance. 










We start the simulations with an initial stellar population of an age of 2 Gyr, 
its mass distribution is determined by a galaxy dynamics model 
().
As cooling flows develop, a circumunuclear disk forms naturally in the galaxy
center, in which the gas is cold and over-dense. Star formation is then inevitable. 
Therefore, we have two stellar populations in our modeled galaxy:
(1) the initial stellar population which is 2 Gyr old at the beginning of the simulation, 
and in which all massive stars have died, 
while its secular evolution (i.e., AGBs and SNe Ia) is considered in our model.
(2) a new stellar population actively forming
in the central cold gas circumnuclear disk during the simulations
and producing SNe II that plays an important role in the metal enrichment.


We calculate time-dependent nucleosynthesis output returned
to the ISM by evolving stars of a simple stellar population (SSP)
with “CELib”, an open-source software library for chemical evolution 
<cit.>.

For the initial stellar population, 
we assume a SSP that consists of stars of identical age and chemical
composition, and follows the <cit.> initial
mass function (IMF) with a stellar mass range of 0.1 ­- 120. 
The SSP is assumed of metallicity Z=1.5Z_⊙ 
(where Z_⊙=0.0134 is the solar metallicity as quoted in Table <ref>).
This abundance matches that expected for an elliptical as massive 
as the modeled one, which should be supersolar  (e.g., ).
The metallicity-dependent stellar lifetime 
is taken from <cit.>.  The stellar “yields”, 
the amount of newly synthesized and ejected elements via 
stellar evolution, are adopted from <cit.>, 
<cit.> and <cit.> 
for high mass stars (Type II SNe), massive AGB stars with 
M>6, and low mass AGB stars respectively. Our 
time-dependent chemical enrichment calculation shows that the 
nucleosynthesis output enriches the ISM mainly via two phases,
an early (t ≲ 10 Myr) phase driven by Type II SNe, and 
a subsequent late phase (t ≳ 40 Myr) driven by AGB stars 
(cf. Equation <ref>). 
Since our simulation starts at t_ age=2 Gyr, there will be only AGBs,
and its time-dependent abundance is used in our model as shown in Figure <ref>, 
while the averaged metal abundance of  AGB winds

is also summarized in Table <ref> for the readers' reference.









The initial stellar population also contributes metals via Type Ia SNe 
(cf. Equation <ref>). 









We adopted the ejecta mass and the nucleosynthesis yields of SN Ia from the N100 model of <cit.>.  According to three-dimensional high resolution hydrodynamical simulations of SN Ia explosion model of <cit.>, each SNIa explosion distributes 1.35 Msolar ejecta in total, and the chemical composition of the ejecta is Fe, Silicon and Oxygen dominated as summarized in Table <ref>.

The new stellar population contribute metals, unlike the initial stellar population,
mostly in its early phase evolution via massive stars. 
For simplicity, we only consider the metals ejected by Type II SNe.
For such new stellar population formed by gravitational instability 
in Toomre unstable circumnuclear disks, we adopt
the  IMF, 
but with the top-heavy index of 1.65 advocated by 
<cit.> and <cit.>, 
which is based on the observed stellar disks in the 
MW and M31. We adopt a stellar mass range of 1 ­- 50 
(cf. Equation <ref>).
The time-averaged metal abundance of its SN II phase

is summarized in Table <ref>.



As in <cit.>, we also consider the cosmic accretion (i.e., CGM infall) 
of low-metallicity gas (cf. Table <ref>) through the galaxy outskirts, 
which will dilute the metallicity as it falls into the galaxy and mixes with the ISM there.  
In addition, we also recycle the metals via the BAL winds which are injected back 
to the galaxy by the central AGN.

Those two processes do not appear explicitly in the equations above, but are 
implemented via the designed outer and inner boundary conditions, respectively.

In this way, we can track the chemical evolution, i.e., the metal enrichment, 
transportation and dilution, throughout massive elliptical galaxies.

After AGN bursts there are significant mass outflows from the galaxy.


To evaluate the metal abundance of the BAL winds, we keep tracking the chemical composition
in our AGN “sub-grid” model <cit.>. 
As  the mass source of the AGN fueling is accreted via 
the inner boundary, we record all the metals that pass through to the galaxy center, 
so that we are able to calculate the averaged metal abundance of the black hole accretion disk 
at any given time. The BAL winds are injected (at the inner edge) back to the galaxy 
with the instantaneous metal abundance of the black hole accretion disk. 

Finally, we assume the initial ISM is of the solar abundance  <cit.>.

We must emphasize that we have not yet corrected for the depletion of metals
onto dust, and thus our characterization of element abundance is the total
(dust plus gas phase) mass in that element per unit volume in units
of the local mass density. Especially in the cold gas component, 
the fractional depletion of refractory elements onto dust grains
may be a large correction <cit.>.  We reserve it to our future work. 




§ MODEL IMPROVEMENTS
 
Compared to <cit.>, the model setups have been modified as follows 
(besides adding the metal tracers) to facilitate the metallicity changes, 
to improve the model self-consistency, and to better fit the observed galaxy properties:



I. the initial black hole mass is increased by 
a factor of 2  (i.e., to M_ BH=6.7×10^8 M_⊙) 
so as to better match the observed M_ BH-σ relation 
<cit.>, 
while the rest of the galaxy model parameters are the same as in <cit.>, 
i.e., total stellar mass at the present time M_⋆=3.35×10^11M_⊙, 
the effective radius R_e=7 kpc, the galaxy ellipticity η=0.2,
and the dimensionless parameters for the total gravity mass =20 
and for the length scale ξ=20 of the modeled galaxy, 
which results in a central projected stellar velocity dispersion of ∼280 km/s.


II. in order to describe the ordered rotational velocity 
of the stellar component we adopted the usual Satoh decomposition
in our previous papers, with a constant rotation parameter k
<cit.>. 
Here we add some flexibility to the modeling by considering 
a central value (k_0) for the k parameter, 
and we reduce it in the outer parts of the galaxy as

    k = k_0 · e^-r/R_e,

where we adopt k_0=0.25.
In this way we can have some control over the amount 
of angular momentum stored in the external regions of the galaxy.


III. we adopt the stellar (AGB)
mass loss formula in <cit.> 
 and adopt the coefficient A=3.3 to match the assumption of the  IMF, i.e.,
 
    Ṁ_⋆ = 10^-12 A ·M_⋆/M_⊙(t_ age/12  Gyr)^-1.3   M_⊙/ year,

where t_ age is the age of the initial stellar population, 
and  the fitting formula is valid when  t_ age≥ 2 Gyr.
As usual, the SN Ia rate is evaluated as 
 
    R_ SN = 0.16×10^-12 h^2 L_ B/L_⊙(t_ age/12  Gyr)^-1.1 year^-1,

where h is the Hubble constant in units of 70  km  s^-1 Mpc^-1, 
and the B-band luminosity L_ B is derived by adopting a mass-to-light ratio 
Γ≡ M_⋆/L_ B=5.8 in units of the solar value,
as appropriate for an old stellar population with a Kroupa IMF 
<cit.>.






IV. the initial mass function of the newly formed stars in the Toomre unstable disks 
is assumed to have a top-heavy profile based on both observations and
theoretical expectations
<cit.>, we assume 

    dN/dM = N_0/M_⊙(M/M_⊙)^-1.65

with a mass range of 1 ­- 50.
Such an IMF gives that ∼60% 
of the total new star mass is in massive stars (M>8M_⊙), 
which will turn into SNe II in a timescale of ∼2×10^7 year.
We assume each SN II leaves a neutron star of 1.4M_⊙ and 
ejects 10^51 erg of energy into the ISM (cf. ).





V. we allow for the fact that the atomic heating/cooling function (S_ line)
is approximately proportional to the local metallicity (Z ≡ 1-(X_1+X_2)/ρ)
in units of the solar value Z_⊙ = 0.0134, 
i.e., in the net heating rate H - C = n^2 (S_ comp + S_ brem + S_ line),  


    S_ line =  10^-23a  + b  (ξ/ξ_0)^c 1 + (ξ/ξ_0)^cZ  Z_⊙,

S_ comp and S_ brem 





are as usual the Compton heating/cooling and the bremsstrahlung loss, respectively
(detailed description of the heating/cooling functions above can be found
in ).


VI. we use the same CGM infall profile as in <cit.>, 
however, since we start the simulation when the galaxy is 2 Gyr old, 
we now also start the CGM infall at that time. 
The accumulated mass due to the CGM infall before 
t_ age = 2 Gyr is now taken as the initial remnant ISM,
which is ∼ 10^10 M_⊙.


VII. the AGN feedback wind efficiency ϵ_w is reduced from 0.005 to 0.004.




Finally, we solve the ISM hydrodynamical equations, together with the metal tracers, using the grid-based
 code <cit.> in spherical coordinates. 
We assume axi-symmetry but allow rotation (a.k.a. 2.5-dimensional simulation).  
The outer boundary is chosen as 250 kilo-parsec to enclose the whole massive elliptical galaxy, 
the inner boundary R_ in is set to be 2.5 parsec to resolve the Bondi radius.
We use a logarithmic grid (Δ r_ i+1/ Δ r_ i = 1.1) to 
divide the radial axis into 120 discrete cells. The azimuthal angle θ 
is divided into 30 uniform cells and covers an azimuthal range from 0.05π to 0.95π. 
The numerical solver for the gas dynamics is composed 
by the combination of the HLLE Riemann Solver, the PLM reconstruction 
and the second-order van Leer integrator.





§ RESULTS

Before we address the simulation details, it is helpful to discuss the 
“closed-box problem”, 
i.e., what metal composition could be expected if the metal enrichment is only contributed 
by the secular evolution of the initial old stellar population, i.e., via AGBs and SNe Ia?
We can see from Figure <ref> that the metal composition of AGB winds 
provides only weak evolution after t_ age=2 Gyr (when the simulation starts),
and we use a fixed metal composition for SN Ia ejecta.
Moreover, Equations <ref> & <ref> give that
the ratio between the mass return rates of SNe Ia and AGBs evolves weakly 
from 0.008 (at t_ age=2 Gyr) 
to 0.012 (at t_ age=13.7 Gyr). 
Therefore, the secular stellar evolution (AGBs + SNe Ia) 
alone will result in a characteristic metallicity in a narrow range of 2.1 Z_⊙ < Z < 2.4Z_⊙ 
throughout the host galaxy.





In addition, as mentioned in <ref>, SNe II contribute to the metal enrichment 
when there is star formation,
while star formation mainly occurs in the circumnuclear disk, 
where most of the infalling gas is circularized because of angular momentum barrier.
Ideally, local metallicity could be up to 12.7 Z_⊙ if SNe II dominate the metal
enrichment (cf Table <ref>). 

We will demonstrate in the rest of this section that the metal enrichment 
in the circumnuclear disk is the key to
understanding the chemical evolution of the modeled galaxy.
Therefore, it is worthwhile to analyze the circumnuclear disk in advance 
(see <ref>). 

Following the physical sequence, i.e., 
formation of the circumnuclear disk, star formation, mass inflow, BAL winds, and
then metal transportation, we perform in <ref> detailed analysis 
on the star formation and black hole feeding processes.

In <ref>, we present the overall gaseous mass and metal budget.

In <ref>, the spatial distribution of metals is presented, 
and we will see the role of the BAL winds in transporting metals throughout the modeled galaxy.
Finally, the radiative features of the metal-enriched hot ISM are calculated in <ref>.











 §.§ Circumnuclear disk
 
In Figure <ref>, we present the hydrodynamical properties 
of the ISM during the burst event at t_ age=12.1 Gyr. 
The gaseous counterpart of the circumnuclear disk is found
of a size ∼150 parsec, and of mass ∼2.5×10^8M_⊙

(see Figure <ref> for the stellar counterpart). 
It is shown in the polar regions that the high-speed BAL winds 
launched by the central AGN
can be heated up to 10^8 K, which should be capable in emitting X rays. 
We have designated as BAL wind gas flowing out
of the galaxy center with radial velocity ≥ 1000 km/s.
Also, metal-rich gas,  recycled from circumnuclear disk, will be transported 
through the galaxy by virtue of the high-speed BAL winds.



In Figure <ref>, we show the cumulative star formation 
in the circumnuclear disk, which  is ∼ 6.5×10^9 M_⊙ at t_ age=12.1 Gyr. 
The stellar disk of the newly formed stars is of a thin geometry, 
and its size is ≲150 parsec, which coincides with the cold gaseous disk 
(cf. Figure <ref>). 
Recall that we have two star formation algorithms
in our model setup: one is based on the Toomre instability 
and the second on the local Jeans timescale,
with the latter allowed only in the densest gaseous zones 
with hydrogen number density n_H>10^5  cm^-3 <cit.>.  
So, most of the star formation is confined to the innermost region of the circumnuclear disk.
In the simulation, we are able to track the age and the in-situ metallicity 
of the newly formed stars, as shown in Figure <ref>, 
the averaged age and metallicity of the new stars are 7.52 Gyr 
and Z/Z_⊙=5.82, respectively.












 §.§ Star formation vs. mass accretion
 
As we discussed previously, 
it is the Toomre instability that primarily triggers star formation 
and allows mass inflow simultaneously.
So it is unavoidable to expect a coincidence between AGN bursts and star formation
<cit.>.

In Figure <ref>, we plot the history of star formation rate Ṁ_⋆^+ (blue line) 
and the mass inflow rates (that through the inner boundary Ṁ_ in in yellow line, 
and that onto the black hole horizon Ṁ_ BH in green line, respectively). 
As expected, the AGN activities are usually coincident roughly with star formation, 
and the instantaneous star formation rate Ṁ_⋆^+ is typically much larger than 
the mass inflow rate through the inner boundary Ṁ_ in
(or the black hole accretion rate Ṁ_ BH).
That is, star formation consumes most of the cold gas within the circumnuclear disk 
before it can fall onto the central supermassive black hole.
This also means that SNe II in the new stellar population can play a very important
role in enriching metals (see <ref> for detailed analysis). 


In Figure <ref>, we present the scatter plot 
between the black hole accretion rate Ṁ_ BH
and star formation rate Ṁ_⋆^+,  
where we bin the timing data in Figure <ref> with adjacent time intervals 
of Δ t = 10^6 and 10^7 year, respectively, 
as shown with the blue (open) and yellow (filled) cycles. 
We find a tight linear correlation of Ṁ_⋆^+ ∼ 7.7 Ṁ_ BH. 
Note that no correlation can be found before binning the data,
due to the time lags between star formation and black hole accretion.
To evaluate the time lags, we perform the calculation below,

    F(Δ t_ lag) = ∫Ṁ_⋆^+(t)·Ṁ_ BH(t+Δ t_ lag) dt/T/∫Ṁ_⋆^+(t) dt/T ·∫Ṁ_ BH(t) dt/T

where T=12 Gyr is the duration of the simulation. 
The results are shown in Figure <ref>. 
It indicates that the time lag between star formation and black hole accretion
is about few 10^6 year. 







 §.§ Overall chemical budget
 
The cumulative chemical yields from the initial (old) stellar population 
and the new (young) stellar population 
are shown in Figure <ref>
with the solid and dashed lines, respectively.  During the simulation, 
the old stellar population ejects in total ∼ 3.4×10^10 M_⊙ 
of total mass into the ISM (∼ 10%M_⋆; cf Equations 
<ref> & <ref>). 
The cumulative star formation is ∼ 6.7×10^9M_⊙ in total, 
of which 58% is injected into the ISM via SNe II (cf Equation <ref>), 
i.e., the new stellar population contributes ∼1/10 of the total stellar mass loss
(see the black lines in the top panel), 
while it dominates the chemical yields of Neon and 
contributes significantly to Oxygen as shown in the middle panel. 
The new stellar population also contributes a significant fraction of heavy metals 
when compared to that of the initial stellar population as shown in the bottom panel.
Note that star formation mainly occurs in the circumnuclear disk.
















With the metal sources and the flow pattern as described above, we are ready to analyze the 
metal enrichment and its transportation throughout the modeled galaxy.
As shown in Table <ref>, 
SN II ejecta is enriched in alpha elements (especially Ne and O);
SN Ia ejecta is of the highest mass fraction of Fe and Si; 
AGBs produce a relatively high fraction of Nitrogen.
So it is possible for us to track the stellar evolution by using the elements mentioned above.













 §.§ Spatial distribution of metals
 




In Figure <ref>, we present the spatial distribution of the metallicity Z 
and individual chemical species (Fe, Ne, C, O, S) during the burst event at t_ age=12.1 Gyr.
We find that the characteristic metal abundances are different from place to place, 
accordingly, the modeled galaxy can be divided into three parts: 
(1) the circumnuclear disk, i.e., (r=0.10 kpc, θ=π/2) and the surrounding region; 
(2) the BAL region, i.e., (r=0.05 kpc, θ≤π/6); 
(3) the main body of the galaxy on the length scale of 10 kpc, i.e., (r=10 kpc, θ=π/2).
We sample the metal abundances from the representative locations above 
and show the results in Table <ref>. 
The high metallicities found in the BAL region and in the cold disk 
are consistent with the SDSS observations 
(e.g., ).















As shown in top-middle panel of Figure <ref>, 
the main body of the galaxy shows the highest Fe abundance 
(e.g. at r∼10 kpc, θ=π/2).
As we discussed previously,  the secular stellar evolution 
of the old stellar population dominates the metal enrichment
since there is no star formation in the main body of the galaxy. 
The typical value of the metallicity Z is ∼2.3Z_⊙, 
and it follows the abundance pattern of AGBs + SNe Ia (cf. Table <ref>)
as discussed at the beginning of this section.

As shown in top-left panel of Figure <ref>, 
the circumnuclear disk has the highest metallicity,  
with a typical value up to ∼8 Z_⊙ in the innermost region
and decreasing to ∼2.3 Z_⊙ at its outer edge. 
Obviously, SNe II play a crucial role in the metal enrichment within the circumnuclear disk.
Note that the disk here seems much more puffy than its geometry seen 
in the hydrodynamical counterpart (cf Figure <ref>),
which is because the massive SN II explosions contaminate the chemical composition 
above and below the disk.
We use Neon to track the contributions of SNe II to the metal enrichment 
in the upper-right panel, we can see that the spatial distribution of Neon 
follows well with the biconical structure of the BAL winds 
(cf Figure <ref>; 
alpha elements shown in the bottom panels of Figure <ref>
also share the similar patterns). 
As we mentioned, the AGN feeding process is mainly via the circumnuclear disk. 
Therefore, the BAL winds would naturally inherit the high metallicity of the circumnuclear disk, 
and spread the metal-rich gas in a fashion of biconical winds.
The metallicity of the BAL winds is up to ∼8 Z_⊙ and 
on average (mass-weighted) it is ∼3.96 Z_⊙, because it is diluted by the ISM while the wind is propagating through its host galaxy. Detailed metal composition of the BAL winds (during three representative burst events, at t_ age=4.0, 8.1, 12.1 Gyr, respectively) can be found in Table <ref>,
as previously, we define the BAL winds as the components 
moving with radial velocity ≥1000 km/s.

We also see at the early epochs that the galaxy outskirts has extremely low metallicity
because of the dilution due to the low-metallicity CGM infall (not shown here). 
When such gas cools down and falls onto the galaxy, 
it would also dilute the metal abundance in the circumnuclear disk. 
























 §.§ Metals seen in X-rays
 
As a result of metal enrichment, the radiative 
cooling and photoionization heating of the ISM are significantly enhanced.
We pay special attention to finding  observable features which can be used to test our model. 
In the rest of this section,
we will perform detailed analysis on the radiative features of the hot ISM, 
and illustrate the advantages of computing detailed metallicity. 
The post processing is based on the hydrodynamical data. 
The frequency-dependent radiation is calculated 
using the atomic database  (version 3.0.9; 
assuming collisional equilibrium),
for example, at any given time and location, 
we retrieve the ISM properties, i.e., its density, temperature,
and the detailed abundance of the 12 chemical species.
Then, we take such output from the  simulation 
as the input of ,
ultimately, the latter would give us the emissivity in the energy band of interest. 










In Figure <ref>,  we plot the representative X-ray spectra 
for the hot ISM and BAL winds as sampled in Table <ref>.
The BAL sample is of a typical temperature ∼2×10^7 K,
and that of the hot ISM is ∼5×10^6 K. 
With the metallicity-dependent emissivity, 
we are able to synthesize the radiative features of the hot gas.


In the top two panels of Figure <ref>, 
we plot the X-ray luminosities and the mass of the hot ISM. 
We can see that  the ISM X-ray luminosity decreases
systematically from redshift z=3.2 to 0, which is mainly because of 
the decrease of gas content in the modeled galaxy. 



In the lower panels of Figure <ref>, we analyze  
the X-ray-emission weighted properties of the hot ISM. 
We average a quantity of interest P in the manner below,

    <P> ≡∫  j_ν P  d V/∫  j_ν d V.

where j_ν is the X-ray emissivity in the energy band of 0.3-8 keV. 
In the volume integrals above, we exclude the disk region 
(r<200 parsec and π/3 < θ < 2/3π) 
and the densest zones with n_ H>10^2  cm^-3,
i.e., which corresponds approximately to a column density of 10^23  atom/ cm^2 
and may be optically thick to X-rays <cit.>; 
The X-ray-emission weighted ISM temperature is presented in the second panel from bottom, 
the typical temperature is ∼ 0.7 keV which is comparable to the virial temperature 
on the length scale of 10 kpc. 
The averaged metallicity is shown in the bottom panel, 
of which the typical value is ∼2.5Z_⊙,
while with Fe abundance approximately ∼6× the solar value.









§ CONCLUSIONS AND DISCUSSIONS
 
We have developed the  code as an instrument for exploring 
the evolution of massive elliptical galaxies at high spatial resolution 
down to and within the fiducial Bondi radius, which enables us to evaluate the 
black hole feeding and feedback processes in a self-consistent way, 
and to track the coevolution between the supermassive black holes
and their host galaxies. 
In the  simulations, we paid special attention to
the internal secular stellar evolution 
which plays a crucial role in driving the galaxy evolution
as important sources of mass and energy for the ISM. 
As an increment to the  code, in this paper, 
we track the recycled mass from the stellar physics 
with its detailed chemical abundances. 
Such an increment  provides us unprecedented details of the ISM properties, 
which are testable by direct comparison with observations. 

To trace the metal enrichment, transportation and dilution processes, we solved 12 additional continuity equations dedicated to H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe and Ni respectively. 
The metal yields, from AGBs and supernovae of type Ia and II,  
are calculated based on standard stellar physics. 
The chemical species are assumed to co-move once after they are injected into the ISM, 
and they are naturally mixed following the fluid motion. 




As expected, a Toomre unstable circumnuclear disk forms in the galaxy center 
with a size of ∼150 parsec, 
which plays a crucial role in the chemical evolution of its host galaxy
— it is where the metals are condensed, further enriched and recycled.
The half mass radius of the new stars is roughly 20 parsec.
More specifically, the massive stars formed in the cold gas disk 
supplement a significant fraction of metals.
From Figure <ref> we see that of order 10^6 massive stars
will be formed over the lifetime of a massive galaxy. 
In their death throes, they will produce core collapse SNe II spewing out 
alpha rich chemical products. 
As a result, the metallicity of the disk in the innermost region reaches up to ∼8 Z_⊙
(Figure <ref>). 
Such metal-rich gas will be captured by the supermassive black hole, 
and then much of it will be recycled back to its host galaxy by virtue of high-speed BAL winds. 
The latter can be readily observed. 
We found that the simulated metallicity  in the BAL winds could be up to ∼ 8 Z_⊙, 
while that of its host galaxy is ∼ 2.3 Z_⊙. 
The X-ray emitting hot gas is very metal enriched with a fluctuating value 
typically near 3 Z_⊙.


Our results are well consistent with the clear correlation found observationally 
between quasar luminosities and their nuclear gas metallicity 
<cit.>.
<cit.> analyzed the metallicity of quasar broad line region (BLR)
with a large SDSS sample. They found that the metallicity in the most massive galaxies
is ∼ 2 Z_⊙, while the BLR metallicity is 0.3 ∼ 1.0 dex larger 
than their host galaxies and it does not evolve with cosmic time.
It indicates that the metal enrichment is due to recent star formation
rather than secular stellar evolution, as we found in this paper.
In the  simulation, the spatial resolution is close to the BLR length scale,
we also found clear correlation between star formation and AGN bursts
(but lagged in time by roughly 10^6 years),
and thus significant metal enrichment in the galaxy center.
The ratio of star formation in the circumnuclear disk to 
accretion towards to the central supermassive black hole is 2.3, 
and of that amount the bulk (62%) is blown out in the BAL winds
with a fraction of 38% finally accreted onto the black hole.


Our simulated metallicity, both for the innermost central region 
and for the main body of the modeled galaxy, is larger than 
derived from X-ray observations 
(e.g., see ,  for the abundance of the hot ISM of early-type galaxies). However, there are a number of uncertainties when deriving abundances from the X-ray spectra (see  for a review). Also, 
dust depletion of metals could be a significant factor 
(e.g., see  for a recent analysis). Additionally we assumed
a metallicity of 1.5Z_⊙ which seemed reasonable given that
the velocity dispersion of our test galaxy is roughly double that in the MW spheroid.
However, as noted earlier, <cit.> would give our test galaxy 
a Fe abundance (in the old stars) slightly less than the solar value.


One important caveat should be noted in our discussion of the expected
chemical abundances in the X-ray emitting gas. We allowed for CGM infall
of metal poor gas, but the amount of infall in our fiducial simulation was
only 1/12 of the stellar mass of the initial galaxy, and even somewhat less than
that (∼1/20) when we corrected for starting at 2 Gyr. This could be too low
an estimate and, in any case, will vary from case to case depending on the
detailed cosmic environment. To check on the sensitivity to this component
we have run a case with three times higher CGM inflow (1/4 of the initial 
stellar mass corrected to ∼1/7 of the stellar mass at t = 2 Gyr). In this test run
the CGM inflow dilutes the metallicity enrichment to such an extent that the
final X-ray-emission weighted metallicity 
drops to ∼1.5 times solar by roughly a factor
of two from our fiducial run. Thus the expected metallicity of the hot Xray
emitting gas has an almost factor of two uncertainty depending on the 
cosmic environment — an uncertainty that can only be addressed by simulations
that are both cosmologically correct and have very high internal spatial
resolution.
We leave it to our future work to explore the consequences of higher 
CGM inflow (at 1/7 of the initial stellar mass) as indicated by cosmological simulations.


The massive stars formed in the Toomre unstable central cold gas disk 
should have dramatic observable consequences, besides the metal enrichment.
For example, the SN II going off in the central disk will also produce
a copious X-ray output. 
However, those SN remnants occur in such dense regions 
that even the X-ray may not be observable
except from the “runaway” star explosions and these are not 
allowed for in the current simulations.

In addition, the (million) neutron star and black hole condensed remnants of the exploded
early type stars must be considered. If they survive in the discs, then 
accretion onto them should have major consequences and they would
be dramatically visible in the “E + A” phases when the gas discs have
largely dissipated. But there is another phenomenon which may intervene.
The embedded condensed remnants will interact with the dense gas
disk  and be dragged into the central supermassive black hole 
in a manor labeled Type I migration first discussed by 
 
(; which is important in planet formation). 
Using the Equation 70 of <cit.>, 
we estimate that the migration time is comparable to the Hubble time especially for
10M_⊙ black holes. Using Type I migration underestimates the migration rate
since it does not allow for the fact that the discs are so dense that in places
they are marginally Toomre unstable. This predictable set of processes 
would produce high mass ratio black hole captures 
(a.k.a. “Extreme Mass Ratio Inspiral", or in short EMRI) which, 
while undetectable by LIGO, 
would be candidates for LISA detection (e.g., ). 
A zeroth order estimate of the rate gives 10^2.5  events/year within 1000 mpc.



In this paper, we did not consider the effect of runaway stars 
in the star-forming circumnuclear disk, which could  be 
important in spreading metals in the galaxy center. 
We also did not  yet consider the dust effects in depleting metals and obscuring radiation, 
which could alter the radiative features of the circumnuclear disk significantly.
In the near future, we will introduce the above effects into 
the  simulations, aiming to interpret/predict 
the observed/observable radiative features of the “E+A” phenomena. 





§ ACKNOWLEDGEMENT

We thank Kengo Tomida for helping us to add tracers into the  code.

We thank Gregory S. Novak for sharing the first 2D version of 
the  code in 2011, which was using .

We thank Jeremy Goodman, James Stone, Pieter van Dokkum,
Nadia Zakamska, Takayuki Saitoh,Tuguldur Sukhboldfor and Charlie Conroy 
for useful discussions. 

ZG is supported in part by the Natural Science Foundation of Shanghai 
(grant 18ZR1447200) and by the Chinese Academy of Sciences 
via the visiting scholar program. 

This work was done during ZG's visit to the department of astronomy 
in Columbia University. 

We acknowledge computing resources from Columbia University's 
Shared Research Computing Facility project,
  which is supported by NIH Research Facility Improvement 
  Grant 1G20RR030893-01, and associated funds from 
  the New York State Empire State Development, Division of Science Technology 
  and Innovation (NYSTAR) Contract C090171, both awarded April 15, 2010. 
Some of the simulations presented
were performed with the computing resources made available via the 
Princeton Institute for Computational Science and Engineering.



  
apj
