Variable protostellar mass accretion rates in cloud cores

TL;DR

This study revisits spherical hydrodynamic models with a polytropic EoS to address the luminosity conundrum in protostar formation, showing that decreasing mass accretion rates over time can explain observed luminosities.

Contribution

It demonstrates that polytropic models with gamma > 1 naturally produce decreasing accretion rates, offering a solution to the luminosity problem in star formation models.

Findings

Mass accretion rate decreases over time in polytropic models with gamma > 1.

Protostellar luminosity declines as accretion rate decreases.

Radial density profiles support the use of polytropic spheres in modeling star formation.

Abstract

Spherical hydrodynamic models with a polytropic equation of state (EoS) for forming protostars are revisited in order to investigate the so-called luminosity conundrum highlighted by observations. For a molecular cloud (MC) core with such an EoS with polytropic index $\gamma$ >1, the central mass accretion rate (MAR) decreases with increasing time as a protostar emerges, offering a sensible solution to this luminosity problem. As the MAR decreases, the protostellar luminosity also decreases, meaning that it is invalid to infer the star formation time from the currently observed luminosity using an isothermal model. Furthermore, observations of radial density profiles and the radio continua of numerous MC cores evolving towards protostars also suggest that polytropic dynamic spheres of $\gamma$ > 1 should be used in physical models.

Equations

$$\frac{\partial \rho }{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}(r^2\rho u)=0,\frac{\partial M}{\partial t}+u\frac{\partial M}{\partial r}=0,$$

$$\frac{\partial M}{\partial r}=4\pi r^2\rho ,\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial r}=-\frac{1}{\rho }\frac{\partial p}{\partial r}-\frac{GM}{r^2},$$

$$\left(\frac{\partial }{\partial t}+u\frac{\partial }{\partial r}\right)\ln \left(\frac{p}{{\rho }^{\gamma }}\right)=0,$$

$$p=K(r,t){\rho }^{\gamma },$$

$$r=k^{1/2}{t}^{n}x\equiv \bar{r}(t)x,$$

$$u=k^{1/2}{t}^{n-1}v(x)\equiv \bar{u}(t)v(x),$$

$$N=\frac{\alpha (x)}{4\pi \mu m_{\mathrm{H}}G{t}^2}\equiv \bar{N}(t)\alpha (x),$$

$$p=\frac{kt^{2n-4}}{4\pi G}\alpha(x)^{\gamma}\equiv\bar{p}(t)\alpha(x)^{\gamma}\ ,\$$

$$M=\frac{k^{3/2}{t}^{3n-2}}{(3n-2)G}m(x)\equiv \bar{M}(t)m(x),$$

$$m(x)=\alpha (x)x^2[nx-v(x)].$$

$$T=\frac{\mu m_{\mathrm{H}}p}{k_{\mathrm{B}}\rho }=\frac{\mu m_{\mathrm{H}}}{k_B}k{t}^{2n-2}\alpha (x{)}^{\gamma -1}\equiv \bar{T}(t)\alpha (x{)}^{\gamma -1},$$

$$\displaystyle v=0\ ,\quad\qquad\alpha=\bigg{[}\frac{n^{2}}{2(2-n)(3n-2)}\bigg{]}^{-1/n}x^{-{2}/{n}}\ ,$$

$$\displaystyle m=n\bigg{[}\frac{n^{2}}{2(2-n)(3n-2)}\bigg{]}^{-1/n}x^{(3n-2)/n}\ ,$$

$$\alpha =A{x}^{-2/n}-3(1-1/n)AB{x}^{-3/n},$$

$$v=B{x}^{(n-1)/n}$$

$$\displaystyle+\bigg{[}\frac{2(2-n)A^{1-n}}{n}-\frac{nA}{(3n-2)}+\frac{(1-n)}{n}B^{2}\bigg{]}x^{{(n-2)}/{n}}\ ,$$

$$v=-\bigg{[}\frac{2m_{0}}{(3n-2)x}\bigg{]}^{1/2},\qquad\ \alpha=\bigg{[}\frac{(3n-2)m_{0}}{2x^{3}}\bigg{]}^{1/2},$$

$$\dot{M_{0}}=k^{3/2}t^{3(n-1)}m_{0}/G\equiv\bar{\dot{M_{0}}}(t)m_{0}\$$

$$\bar{r}=k^{1/2}{t}^n=4\times 10^3\mathrm{AU},$$

$$\bar{N}=(4\pi G\mu m_{\mathrm{H}}{t}^2{)}^{-1}=9\times 10^4\mathrm{c}{\mathrm{m}}^{-3},$$

$$\bar{u}=k^{1/2}{t}^{n-1}=0.30\mathrm{km}{\mathrm{s}}^{-1},$$

$$\bar{M}=\frac{k^{3/2}{t}^{3n-2}}{(3n-2)G}=\frac{0.41}{(3n-2)}{M}_{\odot },$$

$$\bar{T}=\mu m_{\mathrm{H}}k{t}^{2n-2}/k_B=21\mathrm{K},$$

$$\overline{{\dot{M}}_0}=k^{3/2}{t}^{3(n-1)}/G=6.4\times 10^{-6}{M}_{\odot }{\mathrm{yr}}^{-1},$$

$$\displaystyle\dot{M_{\star}}(t)=\epsilon\bigg{[}\frac{t}{0.6\times 10^{5}~{}{\rm yr}}\bigg{]}^{3(1-\gamma)}6.4\times 10^{-6}m_{0}~{}M_{\odot}~{}{\rm yr}^{-1}.$$

$$\displaystyle M_{\star}(t)=\frac{\epsilon t}{4-3\gamma}\bigg{[}\frac{t}{0.6\times 10^{5}~{}{\rm yr}}\bigg{]}^{3(1-\gamma)}6.4\times 10^{-6}m_{0}~{}M_{\odot}.$$

$$L_{\mathrm{acc}}(t)=\frac{G{M}_{\star }(t)\dot{M_{\star }}(t)}{R_{\star }}$$

$$\displaystyle=\frac{\epsilon^{2}t}{(4-3\gamma)}\bigg{[}\frac{t}{0.6\times 10^{5}~{}{\rm yr}}\bigg{]}^{6(1-\gamma)}4.2\times 10^{-4}m_{0}^{2}~{}L_{\odot}\ ,$$

Full text

1 Introduction

The expansion wave collapse solution (EWCS) of Shu (1977) is well known for low-mass star formation111Formation of brown dwarfs (André et al., 2012) and super Jupiter planets (Sumi et al. 2011) can be also described by gaseous spherical collapse models at certain stages (Lou & Shi 2014)., while by extensive observations of actual star-forming molecular cloud (MC) cores, the protostellar mass accretion rate (MAR) disagrees with Shu’s isothermal model of a constant MAR (Kenyon et al., 1994). By observing dust extinction and (sub)millimetre continua, it is realized that the radial profiles of MC core temperature and density do not fit an isothermal model (e.g., Shirley et al. 2002; Harvey et al., 2003; Nielbock et al., 2012). For non-isothermal MC cores, hydrodynamic models with polytropic equations of state (EoSs) have been studied (Fatuzzo et al., 2004; Lou & Gao, 2006; Wang & Lou, 2008) and the special $\gamma=4/3$ case has also been explored (Goldreich & Weber, 1980; Lou & Cao, 2008, Lou & Shi 2016, Li & Lou 2016, Lou & Li 2016). A combination of the singular isothermal sphere (SIS) and polytropic sphere has also been attempted (Curry & McKee, 2000). In massive star formation and the description of giant MCs, static polytropic spheres with $\gamma<1$ (e.g., McKee & Tan, 2002), or even logotrope spheres (e.g., McLaughlin & Pudritz, 1997) were proposed to mimic turbulence effects. Polytropic sphere models with a range of $\gamma$ can broaden applications of the isothermal EWCS, especially in accommodating current observations of MARs and MC core radial density profiles.

For forming protostars in molecular filaments, we recently analyzed dynamic collapses of polytropic cylinders under self-gravity with or without magnetic fields (Lou 2015; Lou & Xing 2016; Lou & Hu 2016a,b). Such cylindrical collapses with axial uniformity and axisymmetry may further break up into segments and clumps along the axis by Jeans instability, leading to chains or binaries of collapsed objects.

One crucial observational test for star formation models is the so-called luminosity problem. Kenyon et al. (1994) found that the protostellar MAR derived from the observed luminosity is too low to accumulate $\sim 1M_{\odot}$ mass in the typical embedded phase of $\sim 10^{5}$ yr in the Taurus-Auriga region, if a constant MAR of the isothermal EWCS is adopted. This luminosity dilemma has been confirmed and aggravated by the Spitzer legacy project “From Molecular Cores to Planet-forming Disks” or “Cores to Disks,” (i.e. c2d) and Gould Belt data (e.g., Evans et al., 2009; Dunham et al., 2013). Radiative transfer analyses involving the EWCS model, a star-disk system, accretion from clumps, or episodic MARs have been performed (Terebey et al., 1984; Myers et al., 1998; Young & Evans, 2005; Dunham et al., 2010; Myers, 2011; Offner & McKee, 2011; Dunham & Vorobyov, 2012; Padoan et al., 2014; Vorobyov & Basu, 2015) towards solving the luminosity problem. Here the episodic MARs as suggested by Kenyon & Hartmann (1995) are supported by theoretical studies via mechanisms of gravitational instabilities (Boss, 2002; Kratter et al., 2010), magnetically driven bursts (Tassis & Mouschovias, 2005), or magneto-rotational instabilities (Zhu et al., 2009; Stamatellos et al., 2012); and by observations of protostellar accretion bursts (Jørgensen et al., 2015) and D/H ratio of water (Owen & Jacquet, 2015). Besides, evolving luminosity has also been studied by slightly changing the SIS model (Offner & McKee, 2011; Dunham et al., 2014) or including turbulence effects in the mass infall of MC cores (Padoan et al., 2014). Alternatively, Lou & Gao (2006) have emphasized that the polytropic gas sphere model gives evolving protostellar MARs, which would offer a physical solution to this luminosity conundrum.

Star formation models are also constrained by the radial profiles of density and temperature from observations of dust extinction (Alves et al., 2001) and (sub)millimeter radio continuum emissions (Motte & André, 2001).222A combination of dust extinction, submillimeter continuum and molecular line features further constrains dynamic polytropic collapse models Lou & Gao (2011); Fu, Gao & Lou (2011). Different radial mass density profiles have been inferred by fitting dust radio continuum observations in the (sub)millimeter bands with assumed temperature profiles (van der Tak et al., 2000; Shirley et al., 2002) or by fitting dust extinction data without prior assumptions of temperature profiles (Kandori et al., 2005; Hung et al., 2010). Various globally static models have been adopted in fitting these observations, e.g., Bonnor-Ebert spheres (Alves et al., 2001; Kandori et al., 2005), power-law models (van der Tak et al., 2000; Shirley et al., 2002; Mueller et al., 2002; Hung et al., 2010; Miettinen & Offner, 2013), double power laws (Beuther et al., 2002), or SISs (Hogerheijde & Sandell, 2000; Shirley et al., 2002; Harvey et al., 2003; Kurono et al., 2013). For the feasibility of dynamic self-similar polytropic sphere models with various $\gamma$, we show their mass density and temperature profiles and make comparison with previous analyses in this Letter.

2 Dynamic Polytropic Gas Spheres

To present physical properties of star-forming MCs and to compare them with observations, the general polytropic (GP) self-similar protostar formation model (Wang & Lou, 2008; Lou & Cao, 2008; Cao & Lou, 2009; Lou & Hu, 2010, Lou & Shi 2014) is adopted here without magnetic field. Hydrodynamic partial differential equations (PDEs) of spherical symmetry are:

[TABLE]

[TABLE]

[TABLE]

with mass density $\rho$, radial bulk flow velocity $u$, thermal gas pressure $p$, and enclosed mass $M$ within radius $r$ at time $t$; $G$ is the gravitational constant. PDEs (1) are equivalent mass conservations, the second of PDE (2) is the radial momentum conservation and PDE (3) is the conservation of specific entropy along streamlines for the GP EoS in the form of

[TABLE]

with $K(r,\ t)$ being a dynamic coefficient depending on $t$ and $r$. We now consider the conventional polytropic (CP) EoS with $n+\gamma=2$ where $n$ is a self-similar scaling index in equation (5) (Lou & Gao, 2006). For the CP EoS, $K(r,\ t)$ above is a global constant.

For a self-similar transformation, all physical variables are expressed as products of dimensional scaling functions and dimensionless self-similar variables, i.e., the radius $r$, radial velocity $u$, particle number density $N$, thermal gas pressure $p$, and enclosed mass $M$ bear the following forms:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with $\mu$ and $m_{\rm H}$ being the mean molecular weight and the hydrogen mass, $n>2/3$ and $\gamma\neq 4/3$. For the CP EoS, we have $n+\gamma=2$ and $n=\gamma=1$ is an isothermal gas. The dimensionless independent self-similar variable $x$ is a combination of $r$ and $t$, and $\alpha(x)$, $m(x)$ and $v(x)$ are dimensionless reduced mass density, enclosed mass and radial velocity, respectively. By mass conservation (1), we derive

[TABLE]

The gas temperature $T$ from the ideal gas law is

[TABLE]

where $k_{\rm B}$ is the Boltzmann constant. All barred factors of $t$ are of the respective physical dimensions.

There are three analytic and asymptotic solutions for CP spheres that are important when comparing our hydrodynamic model with observations.333There exists also a Larson-Penston type asymptotic solutions in the regime of small $x$ for GP gaseous sphere without a central point mass (Lou & Shi, 2014). The first is the singular polytropic sphere (SPS) with

[TABLE]

which is a globally static equilibrium solution for a CP MC core yet singular as $x\rightarrow 0^{+}$. The second one is the asymptotic dynamic solution for $x\rightarrow+\infty$ with

[TABLE]

where $A$ and $B$ are two integration constants (Lou & Shi, 2014). The last one is the dynamic asymptotic free-fall solution towards the MC core centre as $x\rightarrow 0^{+}$:

[TABLE]

where integration constant $m_{0}$ is the reduced enclosed mass as $x\rightarrow 0^{+}$, representing the dimensionless protostellar mass or MAR [see eq. (15)]. Invoking central free-fall solution (14) and recalling eqns (9) and (10), we derive the central protostellar MAR as $x\rightarrow 0^{+}$:

[TABLE]

(see Lou & Gao, 2006; Wang & Lou, 2008, for details).

To estimate the time-dependent dimensional scaling factors of variables in eqns. (5)$-$(9), (11) and (15), we need empirical information for MC cores. According to Myers (2005) and Evans et al. (2009), a typical radius of a star-forming MC core is $\sim 0.01-0.1$ pc ($\sim 10^{3}-10^{4}$ AU); and a typical mean number density in MC cores is estimated by $\sim 10^{4}-10^{5}~{}{\rm cm}^{-3}$. Noting that the dimensionless parts of both radius and number density are around unity in relations (5) and (7), we may choose the length scale and number density scale as

[TABLE]

[TABLE]

respectively. Scales of other physical variables can be expressed in the following manner accordingly:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the mean molecular weight is $\mu\cong 2$ for H2 MCs and $M_{\odot}$ is the solar mass. The relevant dynamic time-scale $t_{d}=0.6\times 10^{5}~{}$yr is estimated by density scale (17).

3 Decreasing Mass Accretion Rate

Referring to MAR scaling (21) at $t_{\rm d}=0.6\times 10^{5}$ yr, the protostellar ($r\rightarrow 0^{+}$ for any time $t>0$, i.e., $x\rightarrow 0^{+}$) MAR (15) can be estimated by $\dot{M_{\star}}(t)=\epsilon\dot{M_{0}}(t)$, viz.

[TABLE]

Here the efficiency coefficient $\epsilon\sim 0.3$, as suggested by Alves et al. (2007), describes the fraction of materials that ends up onto a protostar (see also Evans et al., 2009). By $t$ integration of eq. (22), the protostellar mass $M_{\star}(t)$ becomes

[TABLE]

Bypassing detailed energy transfer processes, we estimate the protostellar accretion luminosity by

[TABLE]

where $L_{\odot}=3.9\times 10^{33}\ {\rm erg\ s}^{-1}$ is the solar luminosity, and the protostellar radius $R_{\star}$ is assumed at $\sim 3R_{\odot}$ with $R_{\odot}=6.9\times 10^{10}\ {\rm cm}$ being the solar radius (Evans et al., 2009). Realistically, $R_{\star}$ may vary during a mass accretion process and $3R_{\odot}$ here is just a first cut. For simplicity, we consider four different dynamic CP spheres with near-static outer envelopes ($B=0$ in eq. (13)) and central free-fall collapsing cores (14), with dimensionless central mass $m_{0}=0.975,1.25,1.15$ and 0.26 for $\gamma=1.0,1.1,1.2$ and 1.3, respectively (Gao, Lou & Wu, 2009). Time evolution of protostellar MAR $\dot{M}_{\star}$, protostellar mass $M_{\star}$ and accretion luminosity $L_{\rm acc}$ are shown in Fig. 1.

A solid conclusion of the top panel in Fig. 1 and MAR (22) is that CP protostellar MARs decrease with increasing time $t$ for $\gamma>1$, and in contrast to constant MARs of isothermal EWCSs ($\gamma=1$ and $n=1$), offer a sensible resolution to the luminosity problem. The evolution of protostellar masses (middle panel in Fig. 1) shows that most protostar masses are accumulated in their early phase. For $\gamma=1.3$, nearly $\sim$ 80% of the stellar mass is accumulated in the first $\sim 10^{4}$ yrs of formation. By protostellar accretion luminosity (24), only for $\gamma>7/6$ dynamic CP spheres do accretion luminosity decrease with increasing time (bottom panel in Fig. 1). Observations show that bolometric luminosity of star-forming MC cores do roughly decrease in their early phases (Class 0, I) of star formation (fig. 13 in Evans et al., 2009), suggesting $\gamma>7/6$ dynamic CP collapses.

With dynamic CP accretion model in hand, we could evolve either from core mass function (CMF) for MC cores to protostellar mass function (PMF) or from initial mass function (IMF) for stars back to PMF; we follow the latter path. By setting up the relation between protostellar luminosity function (PLF) and PMF, we derive model PLF to compare with the observed PLF. By sampling the dynamic CP model results of Fig. 1, one may fit PLF of protostars with envelopes (Evans et al., 2009; Kryukova et al., 2012). By overall energy conservation, the accretion power may be grossly converted to the observed bolometric luminosity with $L_{\rm bol}\sim L_{\rm acc}$. As one cannot follow the luminosity evolution of a single protostar, a statistical analysis of a large enough sample would hint at evolution track of a single protostar formation. We assume a continuous and steady start rate of star formation. Then the IMF of stars (Chabrier, 2005) is convolved (Offner & McKee, 2011; Myers, 2011). Separate comparisons of our model fits of $\gamma=1.3,1.25$ and 1.2 CP spheres with the c2d results (fig. 14 of Evans et al., 2009) are shown in Fig. 2. The intermediate part of the luminosity spectrum is better fitted by our CP dynamic sphere of $\gamma=1.25$, while higher and lower luminosity parts may involve small contributions from $\gamma=1.2$ and $\gamma=1.3$ dynamic spheres, respectively. The shift of peak sites for different $\gamma$ models is due to respective model mean accretion luminosity, related to different accretion histories for pertinent $\gamma$ as shown in Fig. 1. This fitting suggests that a $\gamma=1.25$ sphere is prominent for the early stage of PMF (\;\lower 4.0pt\hbox{{\buildrel\displaystyle<\over{\sim}}}\;0.2~{}{\rm Myr}), with properly weighted contributions from other cores of $1.2<\gamma<1.3$. For illustration, a quick trial (light solid in Fig. 2) weighted for the three $\gamma$s leads to a much better overall PLF fit. Note that the PLF for $L<1L_{\odot}$ has a width of $0.1L_{\odot}$ while for $L>1L_{\odot}$ the width is $1L_{\odot}$. The double peaks in the $\gamma=1.25$ profile and the weighted fit are due to this bin width difference.

In addition to protostellar accretion processes in Fig. 1, accretion termination by bipolar outflows (Shu et al., 1988) or global expansions (Gao & Lou, 2010, Fu et al. 2011) in the late evolution phase for a dynamic CP MC core would also affect the PLF. The observational completeness limit of $L_{\rm comp}=0.05-0.10~{}L_{\odot}$ plays a role in data fit. The lower luminosity limit $L_{\rm ter}=0.08~{}L_{\odot}$ is adopted. More realistic fits to the observed PLF require better determination of the start rate of star formation, the ending mechanism of mass accretion and the CP EoS $\gamma$ distributions.

4 Radial Mass Density Profiles

A successful theoretical model for protostar formation should predict radial profiles of density, temperature and velocity consistent with observations. The radial density profile and the temperature profile vary for CP spheres with different $\gamma$. Power-law radial profile ($\rho\propto r^{-p}$) of cloud density is an empirical model used to fit observed dust continuum distribution. A power-law index of density distribution achieved by assuming certain temperature profiles when fitted to data ranges from $p=1$ to $p=2$ (van der Tak et al., 2000; Shirley et al., 2002), the boundaries of which are logotropic spheres (McLaughlin & Pudritz, 1997) and isothermal spheres (Shu, 1977), respectively. However, without prior temperature assumptions, radial density profiles inferred from dust extinctions would be more realistic (Kandori et al., 2005, Hung et al. 2010), and $p\sim 2.5$ for outer envelopes of cloud cores is inferred by Hung et al. (2010).

With our CP dynamic models, radial profiles of density, temperature and velocity are plotted for $t=t_{d}$ in Fig. 3. Clearly shown in the top panel of Fig. 3, radial density profiles of CP spheres with near-static envelopes and free-fall collapsing cores are broken into two parts: the inner core described by $\rho\propto r^{-1.5}$ and the outer envelope by $\rho\propto r^{-2/(2-\gamma)}$ (see eqns (13) and (12)). The $\gamma=1.2$ case gives a density profile of $\rho\propto r^{-2.5}$ in the outer envelope consistent with extinction data in Harvey et al. (2003), Hung et al. (2010) and Miettinen & Offner (2013). Double power-law fittings to a large sample of high-mass star forming clouds in Beuther et al. (2002) also show that a large fraction of sources have $p>2$ density indices for MC outer envelopes, also hinting at $\gamma>1$ CP spheres. Fig. 3 shows that a CP MC core model of $\gamma=1.2$ naturally provides an increasing temperature towards the centre. Only two typical profiles of $\gamma=1.0$ and $1.2$ are shown in Fig. 3; one can show CP models of $1.0<\gamma<4/3$ bearing similar profiles but with varying slopes as referred to the $\gamma=1.2$ case. And as a generalization to the velocity field of Fig. 3, infalling or expanding envelopes (Lou & Gao, 2006; Wang & Lou, 2008) would accommodate more possible density and temperature radial profiles than near-static envelopes.

5 Conclusions and Discussion

The CP dynamic MC core model of forming protostars with $\gamma>1$ is advanced as a sensible scenario for low mass star-forming MC cores in their early ages, as (1) it offers a reasonable resolution to the luminosity conundrum by allowing a decreasing luminosity with increasing time as a protostar dynamically accretes from surrounding envelope self-similarly; and (2) it features density and temperature radial profiles more realistic as compared with observations of low mass star-forming MC cores. The empirical fact that the accretion luminosity decreases with increasing time corresponds to dynamic CP cores with $\gamma>7/6$, also consistent with the constraint of averaged radial density profiles. For such dynamic CP cores of $\gamma>7/6$, most protostellar masses are actually accumulated in their early phases. This would naturally avoid the luminosity problem.

Our CP dynamic MC core model does not necessarily exclude other concurrent physical processes. As suggested recently (e.g., Kryukova et al., 2012; Vorobyov & Basu, 2015), the actual MAR history of an emerging protostar may vary for low and high mass stars and is most likely a combination of the decreasing MAR due to a CP EoS and bursty accretions arising from, e.g., disc instabilities and shocks (e.g., Shen & Lou 2004, 2006). The latter may be more important in late stages of star formation after the appearance of a disc, which might also be relevant to the low luminosity fit (Evans et al., 2009). Also, either global or bipolar outflows in late stages of star formation provides a possible means of accretion termination. As the CP EoS is more adaptable for low mass star formation, the corresponding decreasing trend of MAR should be more important in the formation history of low mass stars, consistent with the results of Padoan et al. (2014). Further studies connecting CMF, PMF, PLF, IMF, radiation transfer and consequences of discs are desirable.

Acknowledgments

This research was supported in parts by MOST grant 2012CB821800, Tsinghua CCE and THCA, by NSFC grants 10373009, 10533020, 11073014, 11473018 and 51206088, by MOST grant 2014YQ030975, by THISRP (20111081008), and by the Yangtze Endowment, the SRFDP 20050003088, 200800030071 and 20110002110008 as well as 985 grants from MoE and AMD scholarships at Tsinghua Univ. YQL acknowledges support of the China-Chile Scholarly Exchange Program administered by CASSACA.