format: text

text:
# Enhancement of 𝛼-Elements in Dynamical Models of Elliptical Galaxies

## 1 Introduction

Twenty years after the seminal papers on the formation of elliptical galaxies (EGs) by Larson (1975) and Toomre (1977) there is still much disagreement among astronomical community on both the process of formation and evolution of early-type galaxies.

Are EGs formed in a single collapse during a short epoch of activity (monolithic scheme) or continuously by hierarchical mergers at similar levels by several mergers (hierarchical scheme)? Once created, is the stellar population of EGs just passively evolving or do minor merging/accretion events drastically change their characteristics frequently? What is the influence of density environment?

In the monolithic scenario EGs are supposed to suffer dominant star forming activity and consequent chemical enrichment at very early epochs followed by quiescent evolution, marginal episodes of star formation of both internal and external origin may occur under suitable circumstances (e.g. interactions).

In the hierarchical scenario EGs are supposed to be formed by mergers (one to several) of smaller units, each episode inducing star formation and chemical enrichment.

Both scenarios are able to reproduce only part of the observed properties of EGs. For instance the “monolithic scheme”, on which current semi-analytical models of chemical and spectro-photometric evolution of EGs (see below) are based, cannot account for the wide morphological and kinematical variety of EGs, for instance the presence of counter-rotating cores, small disks, recent stellar activity (Schweizer 1990; Schweizer & Seitzer 1992; Longhetti et al. 1998a, b, 2000).

In contrast, despite its success in modeling the dynamical structure of EGs, the “hierarchical scenario” suffers a point of difficulty as far as the chemical properties are concerned. Indeed, EGs often have different ratios between the abundances of some chemical species. They are in general more metal-rich and enhanced in $`\alpha `$-elements with respect to spiral galaxies. In the merging scenario this could be justified by an additional chemical enrichment with strong enhancement in $`\alpha `$-elements must take place. See the reviews by Chiosi (2000, 2001).

Recently Chiosi & Carraro (2001) have investigated by means of N-Body, Tree-SPH models whether the monolithic scheme of galaxy formation and evolution may interpret the observational properties of EGs. They have studied in particular how the history of star formation depends on the initial conditions (density) and/or mass of the proto-galaxy showing that the monolithic scheme can indeed reproduce many fundamental properties of EGs, such as colors, mass to light ratios, the color magnitude relation (CMR), and the so-called Fundamental Plane (see Chiosi & Carraro 2001, for all other details).

In this study we intend to examine whether the same models would also predict chemical properties compatible with the observational data. In particular we estimate their degree of enhancement in $`\alpha `$-elements and check whether they can account for the \[$`\alpha `$/Fe\]-mass relation of EGs (see below).

The plan of the paper is as follows. In Section 1 we discuss the internal contradiction between the kind of constraint imposed by the CMR and the \[$`\alpha `$/Fe\]-mass relationship. In Section 3 we present the main results from the dynamical models, outline the companion chemical models, and report on the results for the chemical abundances and abundance ratios. In Section 4 we present an ad hoc experiment aimed at pinning down the main parameter driving the level of $`\alpha `$-enhanced in model galaxies. In Section 5 we summarize a few basic relations defining the enhancement factor of a single stellar population (SSP) of assigned chemical composition \[X,Y,Z\], present a simple algorithm to calculate the spectral indices $`\mathrm{H}_\beta `$, $`\mathrm{Mg}_2`$, $`\mathrm{Fe}`$ and $`\mathrm{Mg}_\mathrm{b}`$, apply it to our model galaxies, and derive their mean indices. Finally, these latter are compared with the observational data paying particular attention to the $`\mathrm{Mg}_2`$ vs $`\sigma `$ (velocity dispersion) relation and the $`\mathrm{H}_\beta `$ vs $`\mathrm{Mg}_2`$ and $`\mathrm{H}_\beta `$ vs $`[MgFe]`$ planes. In Section 6 we draw some concluding remarks.

## 2 A long lasting contradiction

EGs show important properties that are used as key constraints to theoretical models aimed at understanding the formation mechanism and the past evolutionary history. Chief among these are the CMR and the enhancement in $`\alpha `$-elements suspected to exist in the brightest (most massive) ellipticals.

### 2.1 The Color-Magnitude Relation

The colors of EGs get redder at increasing luminosity and/or mass (Matthews & Baker 1971; Larson 1974; Bower et al. 1992a, b). The relation is tight for cluster EGs (Bower et al. 1992a, b), whereas it is more disperse for nearby field EGs and those in small groups (Schweizer & Seitzer 1992).

The CMR is conventionally interpreted as a mass-metallicity sequence (Faber et al. 1977; Dressler 1984; Vader 1986): massive galaxies reach higher mean metallicities than less massive ones. Starting from the original suggestion by Larson (1974), the mass-metallicity sequence is attributed to the so-called supernova driven galactic wind (SDGW) mechanism halting star formation at an age that depends on the galaxy mass. In short, the overall duration of the star forming activity should be longer at increasing mass of the EGs ($`\mathrm{\Delta }\mathrm{t}_{\mathrm{SF}}\mathrm{M}_\mathrm{G}`$).

The tightness of the CMR is thought of to reflect the age of a galaxy stellar population. The CMR of the Virgo and Coma clusters, where the dispersion is very small, is compatible with most of the stars in EGs being formed at redshift $`z>2`$ which roughly corresponds to a lookback time of about 10 Gyr for a Universe with $`q_0=0.5`$ and $`H_0=70`$ km s-1 Mpc-1. The absolute value of the age is less of a problem here. Nevertheless, secondary bursts of star formation at $`z=0.5`$, forming about 10% of the stellar mass could be also accommodated (Bower et al. 1992a, b). The possibility of secondary activity has also been invoked by Kuntschner & Davies (1998) and Kuntschner (2000) on a different ground. In contrast the more dispersed CMR of the field and small group EGs is compatible with ages of the component stars spread over several billion years. This has also been suggested using the line index diagnostic by González (1993) and Trager et al. (2000a, b) with a tendency for lower velocity dispersion (mass) systems to be younger (see also Bressan et al. 1996).

### 2.2 The \[$`\alpha `$/Fe\]-Mass Relation

Gradients in line strength indices $`\mathrm{Mg}_2`$ and $`\mathrm{Fe}`$ (and others) have been observed in EGs (Worthey 1992; González 1993; Davies et al. 1993; Carollo et al. 1993; Carollo & Danziger 1994a, b; Balcells & Peletier 1994; Fisher et al. 1995, 1996). The hint arises that the indices are stronger toward the center. There is also the controversial question whether the gradients in $`\mathrm{Mg}_2`$ and $`\mathrm{Fe}`$ are the same. According to (Carollo et al. 1993; Carollo & Danziger 1994a, b) the gradient in $`\mathrm{Mg}_2`$ is steeper than the gradient in $`\mathrm{Fe}`$, which is interpreted as indicating that the abundances of $`\alpha `$-elements (Mg, O, etc.) with respect to iron are enhanced – \[$`\alpha `$/Fe\]$`>`$0 – in the central regions. However, Kuntschner (1998) reaches the opposite conclusion and Davies et al. (2001) show for NGC4365 that the ratio \[Mg/Fe\] stays constant within the data range.

Furthermore, limited to the nuclear regions, indices vary passing from one galaxy to another (González 1993; Trager et al. 2000a, b). Looking at the correlation between $`\mathrm{Mg}_2`$ and $`\mathrm{Fe}`$ (or similar indices) for the galaxies in the above quoted samples, $`\mathrm{Mg}_2`$ increases faster that $`\mathrm{Fe}`$ which is once more interpreted as due to enhancement of $`\alpha `$-elements in some galaxies. In addition to this, since the classical paper by Burstein et al. (1988), the index $`\mathrm{Mg}_2`$ is known to increase with the velocity dispersion (and hence mass and luminosity) of the galaxy. Standing on this body of data the conviction arose that the degree of enhancement in $`\alpha `$-elements ought to increase passing from dwarf to massive EGs (Faber et al. 1992; Worthey et al. 1994; Matteucci 1994, 1997; Matteucci et al. 1998).

Understanding these properties rests on stellar nucleosynthesis and star formation. In brief, all heavy elements observed in the interstellar medium (ISM) are produced and ejected by stars, in particular via supernova explosions. Furthermore, iron is mainly produced by the long-lived Type Ia supernovae (accreting white dwarfs in binary systems in the most popular scheme) whereas $`\alpha `$-elements are mainly generated by the short-lived Type II supernovae. Therefore, to enhance the $`\alpha `$-elements there are several avenues:

* (i) A different ratio between Type Ia and Type II supernovae.
* (ii) A different time scale of star formation.
* (iii) A pre-enriched interstellar medium.

The simplest and most widely accepted interpretation is the one based on the different duration of the star forming period and the different contribution to chemical enrichment by Type Ia and Type II supernovae. In this context, it is worth recalling that the minimum mean time scale for a mass accreting white dwarf in a binary system to get the supernova stage is $``$0.5 Gyr. Therefore, the iron contamination by Type Ia supernovae occurs later as compared to the ones of Type II (Greggio & Renzini 1983). It follows from this that with the standard SDGW model and classical initial mass function (IMF), the time scale of star formation and galactic wind must be shorter than about 1 Gyr not to decrease the initial \[$`\alpha `$/Fe\]$`>`$0 (when $`\alpha `$-elements are mostly produced) to \[$`\alpha `$/Fe\]$``$0 (when iron is predominantly ejected). In other words, to reproduce the observed trend of the \[$`\alpha `$/Fe\]-mass relationship, the total duration of the star forming activity ought to scale with the galaxy mass according to $`\mathrm{\Delta }\mathrm{t}_{\mathrm{SF}}\mathrm{M}_\mathrm{G}^1`$, which opposes to the expectation from the CMR.

This point of contradiction cannot be avoided by the standard chemo-spectro-photometric models – both the closed-box ones by Bressan et al. (1994) and those with infall of primordial gas at early epochs by Tantalo et al. (1996, 1998) – because they all stand on the SDWG scheme and are tailored to reproduce the CMR. Therefore, they fail to fit the companion mass-$`\alpha `$-elements relationship. This issue has been addressed by Tantalo & Chiosi (2001) who have carefully examined the effect of key parameters of the classical models, namely the type and efficiency of star formation and the role played by the rate of infall on the degree of enhancement in $`\alpha `$-elements. The main result of this analysis is that the duration of star formation activity and age at which galactic winds set in do not bear very much on the degree of $`\alpha `$-enhancement. This in contrast is determined by the detailed time dependence and intrinsic efficiency of star formation. In other words, a galaxy whose stellar activity predominately occurs at very early epochs (say within the first Gyr), even if followed by a long tail of star formation at smaller efficiency (it does not matter whether continuously or in recurrent episodes), would nowadays appear as enhanced in $`\alpha `$-elements. Obviously, the opposite would hold for a galaxy whose star formation is diluted at low, yet comparable levels of intensity over long periods of time.

Finally, another short comment is worth here. Following Ferrini (2000), the above point of contradiction could stem from the hypotheses that (i) SFR and consequent rate of supernova explosions (mostly type II) are proportional to the gas content (and total galaxy mass $`M_G`$ in turn); (ii) the binding energy of the gas is proportional to $`M_G`$; and finally (iii) the star forming activity stops after the galactic wind has occurred. All of the above hypotheses are oversimplified descriptions of reality, for instance SF is known to occur in molecular clouds, which are a subcomponent of the gas content (in other words the SFR can be a non linear function of the gas content), SF not necessarily halts after the first galactic wind episode and finally the galactic winds can spread over significant intervals of time. All this could destroy the simple relation between the duration of SF and galaxy mass and perhaps alleviate the discrepancy between the CMR and the $`\alpha `$-element problem (Ferrini 2000).

The present approach based upon realistic dynamical models of early type galaxies, in which neither star formation comes to an abrupt halt nor galactic winds suddenly set up, see Chiosi & Carraro (2001) for all details, is perhaps a first step toward a self-consistent formulation of the problem.

## 3 Enhancement of $`\alpha `$-elements in dynamical models

### 3.1 The dynamical models

As we have already anticipated, Chiosi & Carraro (2001) adopting the monolithic scheme have calculated fully hydrodynamical models of early-type galaxies with different total mass $`\mathrm{M}_\mathrm{T}`$. The models included dark and baryonic mass, $`\mathrm{M}_\mathrm{D}`$ and $`\mathrm{M}_\mathrm{B}`$ respectively, in cosmological proportions ($`\mathrm{M}_\mathrm{D}=0.9\times \mathrm{M}_\mathrm{T}`$ and $`\mathrm{M}_\mathrm{B}=0.1\times \mathrm{M}_\mathrm{T}`$), star formation, chemical evolution, energy feed-back, and cooling by radiative processes.

In order to analyze the effect of the initial density of the proto-galaxy three sets of models were considered: (i) a first group (labeled A) whose initial density is given by $`\rho _0200\times \rho _u(z)`$ (where $`\rho _u(z)`$ is the mean density of the Universe as a function of red-shift); (ii) a second group (labeled B), in which the initial density has been arbitrarily set equal to $`\rho _0\rho _u(z)`$; (iii) finally a third group limited to a unique value of the total mass ($`\mathrm{M}_\mathrm{T}=10^9\mathrm{M}_{}`$) in which the initial density has slightly varied above and below the case B value. They are hereinafter indicated by LD (lower density) and HD (higher density). For all other details the reader to refer to Chiosi & Carraro (2001). The main results of this study can be summarized as follows:

(i) Independently of the total mass, galaxies of high initial density (case A), undergo a prominent initial episode of star formation ever since followed by quiescence.

(ii) The same applies to high mass galaxies of low initial density (case B), whereas the low mass ones undergo a series of burst-like episodes that may stretch over a considerable fraction of their lifetime. The details of their star formation history are very sensitive to the value of the initial density as shown by the models of the third group.

(iii) The mean and maximum metallicity increase with the galaxy mass. Therefore these models can account for the CMR of EGs.

(iv) The model galaxies suffer galactic winds ejecting conspicuous amounts of mildly enriched gas into the intra-cluster medium. However, they do not follow the simple SDGW mechanism envisaged by Larson (1974): when the total thermal energy of the gas heated up by supernova explosions (radiative cooling included) exceed the gravitational potential energy of it, the remaining gas is suddenly ejected and star formation is halted (see for instance Bressan et al. 1994; Tantalo et al. 1996, 1998, for more details on this type of models). In dynamical models the situation is more complicated. Looking at radial velocity $`\mathrm{V}_\mathrm{r}`$ of the gas and the escape velocity $`\mathrm{V}_{\mathrm{esc}}`$, at any time there is a layer a which some gas particles may meet the condition $`\mathrm{V}_\mathrm{r}>\mathrm{V}_{\mathrm{esc}}`$. All these gas particle can freely leave the galaxy in form of galactic wind. Therefore galactic winds are a kind of process taking place continuously, at least as long as there are gas particles reaching and exceeding the above condition. Furthermore, the total amount of gas lost in the wind $`\mathrm{\Delta }\mathrm{M}_{\mathrm{gw}}`$ depends on $`\mathrm{M}_\mathrm{T}`$. The following useful relations are found $`\mathrm{\Delta }\mathrm{M}_{\mathrm{gw}}/\mathrm{M}_\mathrm{B}=0.09\mathrm{logM}_\mathrm{T}+1.256`$ for models B and $`\mathrm{\Delta }\mathrm{M}_{\mathrm{gw}}/\mathrm{M}_\mathrm{B}=0.06\mathrm{logM}_\mathrm{T}+0.99`$ for models A. Low mass galaxies are somewhat more efficient emitters of partially processed material.

(v) How does it happen that low mass galaxies may eject large quantities of gas and yet form stars for long periods of time? In a normal or giant elliptical whose gravitational potential well is deep, the gas tends to remain inside the potential well despite being heated up by energy injection by the supernova explosions and stellar winds accompanying star formation. A sort of equilibrium is reached in which gas heating and cooling by radiative processes balance each other, forcing gas to consumption by star formation. In contrast, in a low mass galaxy, owing to the shallower gravitational potential well, the same processes heat up gas, part of it spills over the potential well and part flies to very low density regions so that star formation in the inner regions declines. Eventually, the gas pushed away to large distances but still bound to the galaxy cools down and falls back into the central furnace thus re-fueling star formation. Therefore, this latter keeps going for long periods of time in a series of burst-like episodes till either gas is exhausted or what remains is heated up so much (its gravitational binding energy in meantime gets smaller and smaller compared to the thermal energy) that it eventually leave the galaxy.

(vi) Finally, thanks to their star formation history these models have the potential of improving upon one of the most embarrassing difficulties encountered by the standard supernova driven galactic wind (SDGW) model as far as the chemical abundances and abundance ratios are concerned, i.e. the high $`[\alpha /\mathrm{Fe}]`$ ratio inferred in massive EGs and the long duration of the star forming period required by the SDGW model to explain the CMR.

In Table 1 we summarize a few properties of the models by Chiosi & Carraro (2001) that are relevant to this study (we follow the same identification coding used by Chiosi & Carraro (2001) for the sake of an easy comparison). Column (1) identifies the model. Column (2) gives the total mass $`\mathrm{M}_\mathrm{T}`$ in solar units. Column (3) gives the mean initial density $`\rho _0`$ in units of $`\mathrm{M}_{}\mathrm{kpc}^3`$. Columns (4) and (5) list the total duration of the star formation activity $`\mathrm{\Delta }\mathrm{T}_{\mathrm{SF}}`$ in Gyr, and the age $`\mathrm{T}_{\mathrm{PSF}}`$ in Gyr at which the maximum stellar activity occurs, respectively. Column (6) lists the mean star formation rate $`\mathrm{\Psi }`$ in units of $`\mathrm{M}_{}/\mathrm{yr}`$. Column (7) gives the present-day mass in stars in solar units. Column (8) contains the central velocity dispersion, $`\sigma `$ of baryonic matter in $`\mathrm{km}\mathrm{s}^1`$.

### 3.2 The detailed chemistry

Although the description of the chemical enrichment adopted by Chiosi & Carraro (2001) was fully adequate to the purposes of N-body Tree-SPH simulations, the details of chemical abundances and abundance ratios were out of reach owing to the numerical complexity of those calculations.

To copy with this obvious limitation of the dynamical simulations, we adopt here an indirect approach which albeit simple is sufficiently accurate, i.e. we insert the SFR of dynamical models into the chemical-evolution code of Portinari & Chiosi (1999).

Since chemical evolution affects only the baryonic mass of the galaxy it is sufficient to follow the evolution of an ideal galaxy having the same baryonic mass $`\mathrm{M}_\mathrm{B}`$, initially in form of gas. The chemical models are in the so-called closed-box formulation (see Tinsley 1980) as the dynamical ones. Finally, care has been paid to secure that the yield of heavy elements and gas per stellar generation are the same (or nearly the same) as those used in the dynamical models so that the present-day metallicity and mass in stars of the dynamical and chemical models with the same baryonic mass and star formation history coincide.

Sketch of the chemical models. The equation governing the abundances of elemental species in chemical models of galactic evolution are well known so that no details are given here. The reader is referred to Portinari et al. (1998) for the detailed description of the equations used also in the present study. Suffice here to recall a few important ingredients:

(i) The equations are cast as function of the dimension-less variables

$$\mathrm{G}_\mathrm{g}(\mathrm{t})=\frac{\mathrm{M}_\mathrm{g}(\mathrm{t})}{\mathrm{M}_\mathrm{B}},\mathrm{G}_\mathrm{s}(\mathrm{t})=\frac{\mathrm{M}_\mathrm{s}(\mathrm{t})}{\mathrm{M}_\mathrm{B}},\mathrm{G}_{\mathrm{g},\mathrm{i}}(\mathrm{t})=\mathrm{G}_\mathrm{g}(\mathrm{t})\times \mathrm{X}_\mathrm{i}(\mathrm{t})$$

where $`\mathrm{X}_\mathrm{i}(\mathrm{t})`$ are the abundances by mass of the elemental species $`i`$.111 Fifteen chemical elements are followed in detail, namely $`\mathrm{H}`$, $`{}_{}{}^{4}\mathrm{He}`$, $`{}_{}{}^{12}\mathrm{C}`$, $`{}_{}{}^{13}\mathrm{C}`$, $`{}_{}{}^{14}\mathrm{N}`$, $`{}_{}{}^{16}\mathrm{O}`$, $`{}_{}{}^{20}\mathrm{Ne}`$, $`{}_{}{}^{24}\mathrm{Mg}`$, $`{}_{}{}^{28}\mathrm{Si}`$, $`{}_{}{}^{32}\mathrm{S}`$, $`{}_{}{}^{40}Ca`$, $`{}_{}{}^{56}Fe`$, and the isotopic neutron-rich elements $`nr`$ obtained by $`\alpha `$-capture on $`{}_{}{}^{14}\mathrm{N}`$, i.e. $`{}_{}{}^{18}\mathrm{O}`$, $`{}_{}{}^{22}\mathrm{Ne}`$, $`{}_{}{}^{25}\mathrm{Mg}`$. By definition $`_\mathrm{i}\mathrm{X}_\mathrm{i}(\mathrm{t})`$=1.

(ii) The restitution fractions of the elements $`i`$ from a star of mass $`\mathrm{M}`$ at the end of its life are obtained by means the matrix formalism of (see Talbot & Arnett 1971, 1973, 1975).

(iii) No use of the instantaneous recycling is made (i.e. we take the life time $`\tau _\mathrm{M}`$ of stars of mass $`M`$ into account.

(iv) We consider the separated contributions to nucleosynthetic ejecta by single and binary stars. Single stars can terminate their life quietly becoming white dwarfs and ejecting part of their mass in form of stellar winds or exploding as Type II supernovae. Accreting white dwarfs in binary stars are thought of to be the site of Type Ia supernovae as long ago suggested by Greggio & Renzini (1983) and Matteucci & Greggio (1986). To evaluate the number of Type Ia supernovae we need do fix the minimum ($`\mathrm{M}_\mathrm{B}`$m) and the maximum ($`\mathrm{M}_\mathrm{B}`$M) value for the total mass ($`\mathrm{M}_\mathrm{B}`$M) of the binary systems, the distribution function ($`\mathrm{f}(\mu )`$) of their mass ratios, where $`\mu _{\mathrm{min}}`$ is the minimum value of it, and finally the fraction $`\mathrm{\Lambda }`$ of binaries with respect to the total. It is assumed here that binary stars as a whole obey the same IMF of single stars. We adopt $`\mathrm{B}_\mathrm{m}=3\mathrm{M}_{}`$, $`\mathrm{B}_\mathrm{M}=12\mathrm{M}_{}`$, and $`\mathrm{\Lambda }=0.02`$.

(v) The stellar ejecta are from Marigo et al. (1996, 1998) and Portinari et al. (1998), and incorporate the dependence on the initial chemical composition as described by Portinari et al. (1998).

(vi) Finally, the stellar lifetimes $`\tau _\mathrm{M}`$ are from Bertelli et al. (1994). We take also into account the dependence of $`\tau _\mathrm{M}`$ on the initial chemical composition.

The rate of star formation. The rate of star formation (SFR) of the dynamical models is according to the Schmidt (1959) law applied to each gas particle. Therefore, the total star formation rate is also the Schmidt law applied to the whole system:

$$\frac{\mathrm{dM}_\mathrm{g}(\mathrm{t})}{\mathrm{dt}}=\nu \mathrm{M}_\mathrm{g}(\mathrm{t})^\kappa $$
(1)
where $`\nu `$=1 and $`\kappa `$=1 as in Chiosi & Carraro (2001). The SFR of the dynamical models must be rescaled according to the normalization adopted in the equations of the chemical models. Upon normalization, the star formation rate becomes:

$$\mathrm{\Psi }(\mathrm{t})=\frac{\mathrm{d}}{\mathrm{dt}}(\frac{\mathrm{M}_\mathrm{g}}{\mathrm{M}_\mathrm{B}})$$
(2)
The initial mass function. We have adopted here the same IMF used by Chiosi & Carraro (2001), i.e. the Miller & Scalo (1979) over the mass interval $`\mathrm{M}_\mathrm{l}`$=0.1 to $`\mathrm{M}_\mathrm{u}`$=120 $`M_{}`$. It is worth of interest to look at the fractionary mass stored by the IMF in the mass interval whose stars effectively contribute to nucleosynthesis in a time interval shorter than the age of the galaxy. A suitable choice is the mass range $`\mathrm{M}_\mathrm{n}`$=1$`M_{}`$ to $`\mathrm{M}_\mathrm{u}`$. The above fraction $`\zeta `$ is

$$\zeta =\frac{_{\mathrm{M}_\mathrm{n}}^{\mathrm{M}_\mathrm{u}}\varphi (\mathrm{M})\times \mathrm{M}\times \mathrm{dM}}{_{\mathrm{M}_\mathrm{l}}^{\mathrm{M}_\mathrm{u}}\varphi (\mathrm{M})\times \mathrm{M}\times \mathrm{dM}}$$
(3)
With the slope of the Miller & Scalo (1979) IMF and the above mass limits, $`\zeta 0.25`$.

### 3.3 Chemical abundances and abundance ratios

The chemical evolution of all the models presented in Table 1 has been followed and the results of interest here are best represented by the series of Figs. 1 through 3.

In Fig. 1 we show the time variation of the SFR (top panel), the metallicity Z (mid panel), and the abundance ratio \[O/Fe\] (bottom panel) for the models 3B with $`\mathrm{M}_\mathrm{T}=10^9\mathrm{M}_{}`$ and 6B-LD with $`\mathrm{M}_\mathrm{T}=10^9\mathrm{M}_{}`$ as indicated. As expected, during the period of stellar activity, the metallicity increases, whereas when star formation is over the metallicity remains either constant or slightly decreases by dilution with gas ejected by old stars of low metal content. Looking at the time variation of the abundance ratio \[O/Fe\] – a measure of the $`\alpha `$-enhancement level – we note that

(i) Within the first 0.1 Gyr the abundance ratios are lowered to values comparable to the mean ones under the contamination by the first Type Ia supernovae.

(ii) Even a small amount of stellar activity changes the initial abundance ratios to much lower values.

(iii) Subsequent star formation does not alter the abundance ratios significantly. See for instance the case of the $`10^{12}\mathrm{M}_{}`$ galaxy in which the prominent burst occurring from 2 to 3.5 Gyr changes (increases) the \[O/Fe\] ratio by less than $`\mathrm{\Delta }[\mathrm{O}/\mathrm{Fe}]`$=0.2 (less than a factor of two).

(iv) Once star formation has ceased, almost immediately the effect of Type II supernovae stops and that by Type Ia takes over first rapidly and later slowly decreasing \[O/Fe\] towards the solar value or below this.

In order to answer the question whether or not the present-day stellar content of these galaxies would appear to be enhanced in $`\alpha `$-elements, we look at the cumulative distribution of the number of living stars in different bins of \[O/Fe\]. From an operational point of view a galaxy will be considered dominated by an $`\alpha `$-enhanced stellar content if $`\mathrm{N}([\mathrm{O}/\mathrm{Fe}]>0)/\mathrm{N}_\mathrm{T}>0.5`$, where $`\mathrm{N}([\mathrm{O}/\mathrm{Fe}])`$ is the number of stars in bins with \[O/Fe\]$`>`$0 and $`\mathrm{N}_\mathrm{T}`$ is their total number. The opposite if $`\mathrm{N}([\mathrm{O}/\mathrm{Fe}]<0)/\mathrm{N}_\mathrm{T}<0.5`$. The situation for the two galaxies under examination is shown in Fig. 2 ($`10^{12}\mathrm{M}_{}`$) and Fig. 3 ($`10^9\mathrm{M}_{}`$).

The massive galaxy (Fig. 2), characterized by a rather short ($`2.5`$ Gyr) and prominent burst of star formation activity has about 98% of its stellar content in bins with with \[O/Fe\]$`>0`$. Whereas the low mass, low initial density galaxy (Fig. 3), whose star formation stretched over long periods of time in several recurrent episodes, has the majority of stars in bins with \[O/Fe\]$`<0`$ ($``$44%).

The percentage of stars with $`\mathrm{N}([\mathrm{O}/\mathrm{Fe}])/\mathrm{N}_\mathrm{T}>0.5`$ are listed in column (9) of Table 1. Finally, in columns (3) and (4) and of Table 4 to better described below we list the maximum, $`\mathrm{Z}_{\mathrm{Max}}`$, and the mean metallicity, $`\mathrm{Z}`$, of our model galaxies.

## 4 What does determine the enhancement of $`\alpha `$-elements?

The entries of Table 1 show that the percentage of stars with \[O/Fe\]$`>`$0 decreases with the total duration of the star formation activity $`\mathrm{\Delta }\mathrm{T}_{\mathrm{SF}}`$. However all models but for the case 6B-LD have percentages of $`\alpha `$-enhanced stars higher than 50%. In order to understand under which circumstances (shape of the SFR and $`\mathrm{\Delta }\mathrm{T}_{\mathrm{SF}}`$) a galaxy would build a stellar population with a percentage of $`\alpha `$-enhanced stars less than 50%, we have performed the following experiment. We have taken one of the most unfavorable cases, i.e. a typical galaxy with $`\mathrm{M}_\mathrm{T}=10^{12}\mathrm{M}_{}`$ and hence $`\mathrm{M}_\mathrm{B}=10^{11}\mathrm{M}_{}`$, imposed that the same amount of mass in stars is build up at the end of the star forming activity, and varied both the shape and duration, $`\mathrm{\Delta }\mathrm{T}_{\mathrm{SF}}`$ of the SFR.

We adopt a simple analytical expression for the SFR as illustrated in Fig. 4. This is idealized to a trapezium of sides A, B and $`\mathrm{\Delta }\mathrm{T}_{\mathrm{SF}}`$. By varying the ratio B/A, the shape of the SFR changes from case (1) in Fig. 4 with B/A$`<<`$ 1 (SFR peaked in the past), to a constant SFR with B/A=1 of case (2) in Fig. 4, and finally to case (3) in Fig. 4 (SFR peaked toward the present). The total area of the trapezium (rectangle) yields the total mass $`\mathrm{M}_\mathrm{s}`$ in stars build up in the galaxy. The galaxy under consideration has a total mass in stars of about $`6.4\times 10^{10}\mathrm{M}_{}`$.

The following relationship holds between $`\mathrm{M}_\mathrm{s}`$ and SFR

$$\mathrm{M}_\mathrm{s}=0.5\times \mathrm{A}\left[1+\frac{\mathrm{B}}{\mathrm{A}}\right]\times \mathrm{\Delta }\mathrm{T}_{\mathrm{SF}}$$
(4)
from which we may easily derive an analytical expression for the SFR to be plugged into the chemical code. We have explored the following combinations of the parameters B/A=0.1, 1.0, and 10 and $`\mathrm{\Delta }\mathrm{T}_{\mathrm{SF}}`$=1, 2, 4, 6 and 10 Gyr. The results are displayed in Fig. 5, which plots the present-day percentage of $`\alpha `$-enhanced stars as a function of $`\mathrm{\Delta }\mathrm{T}_{\mathrm{SF}}`$ and the ratio B/A (the dotted line for B/A=0.1, the solid line for B/A=1, and the dashed line for B/A=10). It is soon evident that $`\mathrm{\Delta }\mathrm{T}_{\mathrm{SF}}`$ is the dominant parameter, whereas the ratio B/A plays a marginal role. With the prescription for Type Ia supernovae we have adopted, to fall below the border line of 50%, $`\mathrm{\Delta }\mathrm{T}_{\mathrm{SF}}`$ should be longer than about 6 Gyr and the stellar activity should be uniformly diluted. This is indeed the case of 6B-LD model which finds confirmation by the present analysis. A counter-example is the case of the 7B model with $`\mathrm{M}_\mathrm{T}=10^8\mathrm{M}_{}`$ whose $`\mathrm{\Delta }\mathrm{T}_{\mathrm{SF}}`$6.5 Gyr but the SFR is mostly confined between 2 to 4 Gyr with the peak at about 3 Gyr. The reason for that is the initial density which is slightly higher than that of case 6B-LD (see Table 1 and Chiosi & Carraro (2001) for all details). The total duration of the SFR and not its shape drives the level of enhancement.

## 5 Mean narrow band indices for galaxies

In this section we evaluate the narrow band indices for our model galaxies in order to check how far they can go in reproducing important relationships such as the $`\mathrm{Mg}_2`$ versus $`\sigma `$ (Burstein et al. 1988) and the $`\mathrm{H}_\beta `$ versus $`\mathrm{Mg}_2`$ and/or $`\mathrm{MgFe}`$, see for instance González (1993); Trager et al. (2000a, b).

In presence of enhancement in $`\alpha `$-elements one has to modify the relationship between the total metallicity $`Z`$ and the iron content \[Fe/H\] by suitably defining the enhancement factor $`\mathrm{\Gamma }`$. For more details see Tantalo (1998) and Tantalo & Chiosi (2002).

Let us split the metallicity $`Z`$ in the sum of two terms

$$\mathrm{Z}=\underset{\mathrm{j}}{}\mathrm{X}_\mathrm{j}+\mathrm{X}_{\mathrm{Fe}}$$
(5)
where $`\mathrm{X}_\mathrm{j}`$ are the abundances by mass of all heavy elements but Fe, and $`\mathrm{X}_{\mathrm{Fe}}`$ is the same for Fe. Recasting eq.(5) as

$$\mathrm{Z}=\frac{\mathrm{X}_{\mathrm{Fe}}}{\mathrm{X}_\mathrm{H}}\mathrm{X}_\mathrm{H}\left[1+\frac{_\mathrm{j}\mathrm{X}_\mathrm{j}}{\mathrm{X}_{\mathrm{Fe}}}\right]$$
(6)
and normalizing it to the solar values we get

$$\left(\frac{\mathrm{Z}}{\mathrm{Z}_{}}\right)=\left(\frac{\mathrm{X}_{\mathrm{Fe}}}{\mathrm{X}_\mathrm{H}}\right)\left(\frac{\mathrm{X}_\mathrm{H}}{\mathrm{X}_{\mathrm{Fe}}}\right)_{}\left(\frac{\mathrm{X}_\mathrm{H}}{\mathrm{X}_{\mathrm{H},}}\right)\frac{\left[1+\frac{_\mathrm{j}\mathrm{X}_\mathrm{j}}{\mathrm{X}_{\mathrm{Fe}}}\right]}{\left[1+\frac{_\mathrm{j}\mathrm{X}_\mathrm{j}}{\mathrm{X}_{\mathrm{Fe}}}\right]_{}}$$
(7)
from which we finally obtain

$$\left[\frac{\mathrm{Fe}}{\mathrm{H}}\right]=\mathrm{lg}\left(\frac{\mathrm{Z}}{\mathrm{Z}_{}}\right)\mathrm{lg}\left(\frac{\mathrm{X}}{\mathrm{X}_{}}\right)\mathrm{\Gamma }$$
(8)
which provides also the definition of $`\mathrm{\Gamma }`$. Obviously for solar metallicity $`\mathrm{Z}_{}`$ and solar-scaled mixture ($`\mathrm{\Gamma }=0.0`$) eq.(8) yields \[Fe/H\]$`=0.0`$. In general, when dealing with enhanced chemical mixtures, it suffices to determine $`\mathrm{\Gamma }`$ according to the degree of enhancement and simply re-scale the zero-point for the iron content of the Sun, i.e \[Fe/H\]=-$`\mathrm{\Gamma }`$). Remarkably, enhancing $`\alpha `$-elements simply reduces to decreasing \[Fe/H\].

Table 2 lists the initial chemical composition and the iron abundances corresponding to $`\mathrm{\Gamma }`$=0 and $`\mathrm{\Gamma }`$=0.356, the same as the stellar models and isochrones of Salasnich et al. (2000) on which our indices are ultimately based.

Once $`\mathrm{\Gamma }`$ is known for the stellar population under consideration, to proceed further we need a calibration expressing the indices as a function of it. To this purpose we make use of the calibration by Borges et al. (1995) for $`\mathrm{Mg}_2`$ in which the effect of enhancement is considered. According to Borges et al. (1995), the index $`\mathrm{Mg}_2`$ of a “star” with enhanced content of Mg is given by the relation

$$\mathrm{ln}\mathrm{Mg}_2=9.037+5.795\frac{5040}{\mathrm{T}_{\mathrm{eff}}}+0.398\mathrm{log}\mathrm{g}+$$

$$0.389\left[\frac{\mathrm{Fe}}{\mathrm{H}}\right]0.16\left[\frac{\mathrm{Fe}}{\mathrm{H}}\right]^2+0.981\left[\frac{\mathrm{Mg}}{\mathrm{Fe}}\right]$$
(9)
which holds for effective temperatures and gravities in the ranges $`3800<\mathrm{T}_{\mathrm{eff}}<6500`$ K and $`0.7<\mathrm{log}\mathrm{g}<4.5`$. To make use of eq.(9), we need \[Mg/Fe\] as function of $`\mathrm{\Gamma }`$. The correspondence between \[Mg/Fe\] and $`\mathrm{\Gamma }`$ can be obtained in the following way: from the pattern of abundances by Anders & Grevesse (1989); Grevesse (1991) and Grevesse & Noels (1993) we get the approximated relation

$$\mathrm{\Gamma }0.8\left[\frac{\mathrm{Mg}}{\mathrm{Fe}}\right]$$
(10)
According to which \[Mg/Fe\]=0.0 for $`\mathrm{\Gamma }`$=0 (solar-scaled abundances), and \[Mg/Fe\]=0.444 for $`\mathrm{\Gamma }`$=0.356.

For all remaining indices we use the Worthey et al. (1994) fitting functions in which the enhancement of $`\alpha `$-elements is taken into account only via the re-scaled value of \[Fe/H\] of eq.(8).

Looking at the response of indices to variations of the age, $`\mathrm{log}(\mathrm{t})`$ in Gyr, metallicity $`\mathrm{log}(\mathrm{Z}/\mathrm{Z}_{})`$, and degree of enhancement $`\mathrm{\Gamma }`$ of SSPs, we derive the following analytical relationships,

$$\begin{array}{cc}\mathrm{I}_\mathrm{i}=\left(\mathrm{a}_{1,\mathrm{i}}\mathrm{\Gamma }+\mathrm{a}_{0,\mathrm{i}}\right)\mathrm{log}(\mathrm{t})\hfill & \\ +\left(b_{1,i}\mathrm{\Gamma }+b_{0,i}\right)\mathrm{log}\left(\frac{Z}{Z_{}}\right)+c_{0,i}\mathrm{\Gamma }+\delta _i\hfill & \end{array}$$
(11)
where $`\mathrm{I}_\mathrm{i}`$ for $`i`$=1 to 4 stand for the $`\mathrm{H}_\beta `$, $`\mathrm{Mg}_2`$, $`\mathrm{Fe}`$ and $`\mathrm{Mg}_\mathrm{b}`$, respectively. The coefficients are tabulated in Table 3 for the Borges et al. (1995) calibration. These analytical fits are accurate enough for the purposes of the present study.

Given these premises, we derive the indices for our model galaxies as follows:

(i) With the aid of the cumulative distributions of living stars per \[O/Fe\] bins (see the bottom panels of Figs. 3 and 2 and those for all remaining models not shown here for the sake of brevity), we estimate the mean value of \[O/Fe\] weighed on the number of stars per bin and consider it as an estimate of the enhancement factor $`\mathrm{\Gamma }`$. This is found to weakly correlate with the velocity dispersion of the galaxy $`\mathrm{\Gamma }=0.072log\sigma +0.134`$.

(ii) From the history of star formation (top panels of the same Figures.) we estimate the mean age $`\mathrm{T}_{\mathrm{BS}}`$ of the bulk stars. To do so one has to convert the rest-frame ages of the galaxy models into an absolute age by means of a cosmological model of the Universe. To make things simple we adopt here the Freedman model with $`\mathrm{H}_0`$=50 km s-1 Mpc-1, $`\mathrm{q}_0`$=0 and red-shift of galaxy formation $`\mathrm{z}_{\mathrm{for}}=5`$, which yield the maximum age for galaxies $`\mathrm{T}_\mathrm{G}`$=16.5 Gyr. The age $`\mathrm{T}_{\mathrm{BS}}`$ is therefore $`\mathrm{T}_{\mathrm{BS}}=\mathrm{T}_\mathrm{G}\mathrm{T}_{\mathrm{PSF}}`$, where $`\mathrm{T}_{\mathrm{PSF}}`$ is the time at which the peak of star formation occurred. These quantities are listed in columns (5) and (6) of Table 4.

(iii) Finally, we make use of relations (11) to derive the indices $`\mathrm{H}_\beta `$, $`\mathrm{Mg}_2`$, $`\mathrm{Fe}`$, and $`\mathrm{Mg}_\mathrm{b}`$ for both the maximum and mean metallicity (age $`\mathrm{T}_{\mathrm{BS}}`$ and $`\mathrm{\Gamma }`$ are the same). The results are listed in Table 4.

The $`\mathrm{Mg}_2`$ vs $`\sigma `$ relation. This is shown in Fig. 6 and compared with the observational data by Trager et al. (2000a, b) for a sample of local EGs. The least squares fit of the observational data yields the relation $`\mathrm{Mg}_2=0.186\mathrm{log}\sigma 0.125`$ (where $`\sigma `$ is in km s-1). The theoretical results are plotted separately for each of the metallicities in use: the filled triangles are for $`\mathrm{Z}`$ (a more realistic approximation of reality) and the open triangles are for $`\mathrm{Z}_{\mathrm{Max}}`$. The least squares fit of the data yields $`\mathrm{Mg}_2=0.063\mathrm{log}\sigma +0.001`$ for $`Z`$ and $`\mathrm{Mg}_2=0.094\mathrm{log}\sigma +0.063`$ for $`Z_{Max}`$. The slope is about half of what observed. Where is the source of disagreement? There are several candidates: (i) the velocity dispersion, (ii) the calibration in use, (iii) the set of parameters $`\mathrm{T}_{\mathrm{BS}}`$, $`\mathrm{Z}`$ and $`\mathrm{\Gamma }`$ assigned to the model galaxies. Limiting the discussion to the range of massive models for which data exist, we note:

(i) The velocity dispersion is not likely because it should be decreased by a factor of two and in any case the relation would simply horizontally shift without changing its slope.

(ii) The calibration is more uncertain, even though the one in use here allowing for the maximum effect. In this context one should, however, keep in mind that, contrary to what one may guess, increasing $`\mathrm{\Gamma }`$ (higher enhancement) does not imply a stronger $`\mathrm{Mg}_2`$ because of the counter-reaction by \[Fe/H\].

(iii) Assuming that our velocity dispersions and calibration are correct, inspection of the weight of $`T_{BS}`$, $`\mathrm{\Gamma }`$, and $`\mathrm{Z}`$ on building up the final result reveals that the metallicity dominates. The three empty stars show the effect of arbitrarily increasing $`\mathrm{Z}`$ by a factor of 3, limited to the three most massive galaxies. This implies that the $`\mathrm{Mg}_2`$ vs $`\sigma `$ relation reflects the mass-metallicity relation, rather a sequence in the enhancement level.

Are the mean (and maximum) metallicity reasonable? Is there any plausible effect making a massive galaxy more efficient in building metals? With the IMF adopted here and in the dynamical models, the production of metals is scarcely efficient because of the low $`\zeta `$, Increasing this for all galaxies would simply shift vertically the $`\mathrm{Mg}_2`$ vs $`\sigma `$ relation. However, if $`\zeta `$ increases with the galaxy mass (top-heavy IMF), theory and data could be reconciled. Although this possibility has often been advanced by several authors in different astrophysical contexts – see the arguments given by Chiosi et al. (1998) and Larson (1998) – further exploring the subject is beyond the scope of this paper. Suffice to conclude that dynamical models almost match the observational $`\mathrm{Mg}_2`$ vs $`\sigma `$ relation provided some fine tuning of the chemical yields per stellar generation is applied. Just for the sake of curiosity, we derive the relationship one would obtain by letting the mean metallicity to increase from the present values to $`3\times \mathrm{Z}`$ as the galaxy mass mass increases from $`10^9`$ to $`5\times 10^{12}\mathrm{M}_{}`$. We get $`\mathrm{Mg}_2=0.149\mathrm{log}\sigma 0.070`$, which is about right.

The $`\mathrm{H}_\beta `$ versus $`\mathrm{Mg}_2`$ and $`\mathrm{MgFe}`$ planes. The index-index planes are customarily used to infer age and metallicity of early-type galaxies (González 1993; Bressan et al. 1996; Trager et al. 2000a, b; Kuntschner 2000, 2001). In Fig. 7 we present the $`\mathrm{H}_\beta `$ vs $`\mathrm{Mg}_2`$ and $`\mathrm{H}_\beta `$ vs $`\mathrm{MgFe}`$ planes left and right panels respectively, for our models, single stellar populations of similar metallicity, and observational data for a small sample of EGs in the Fornax cluster by Kuntschner (2000) and the wider sample of nearby galaxies by Trager et al. (2000a). The following remarks can be made:

(i) Our models lie on a strip which roughly correspond to that covered by SSP of different metallicity and age in the range say 15 to 8-10 Gyr. This expected because all models have terminated their star formation far in the past. Even the most extreme case has ceased star formation 6 Gyr ago.

(ii) Also in these planes we notice the same kind of difficulty encountered with the $`\mathrm{Mg}_2`$ vs $`\sigma `$ relation, i.e. our model galaxies do not make enough metals so that even the most massive ones only marginally reach the region crowed by the data. The same kind of arguments brought about to cure the $`\mathrm{Mg}_2`$ vs $`\sigma `$ relation apply here. Increasing the metal content of the model galaxies is less of a problem.

## 6 Concluding remarks

In this paper we have derived the present-day level of enhancement in $`\alpha `$-elements in the stellar content of the N-Body Tree-SPH models of EGs by Chiosi & Carraro (2001), cast light on the main cause of it, derived the mean spectral indices of the model galaxies, and compared them to the observations with particular attention to the $`\mathrm{Mg}_2`$-$`\sigma `$ relation and the $`\mathrm{H}_\beta `$ \- $`\mathrm{Mg}_2`$–\[MgFe\] planes.

Combining the results of Chiosi & Carraro (2001) summarized in Section 3 and those obtained here the following remarks can be made:

(i) The star formation history of EGs depends on their total mass and initial density, i.e. the physical conditions of the environment in which they have formed. High-density systems, owing to their very early and short star forming activity would now appear as old galaxies, whose stellar content is enhanced in $`\alpha `$-elements. Low-density systems of high mass would resemble the previous ones, their stellar content should, however, span a broader range of ages. Most of the component stars should be enhanced in $`\alpha `$-elements but with a broader distribution of abundance ratios \[$`\alpha `$/Fe\]. At decreasing mass the star formation activity should become more diluted in time, the stellar populations should span wider and wider age ranges, and the level of enhancement should progressively tend to solar. The extreme case is reached with the low-mass, low-density systems which could undergo many episodes of star formation up to very recent times, thus resembling young objects with solar or sub-solar partitions of elements. It is worth recalling here that $`\mathrm{\Gamma }`$ or equivalently \[O/Fe\] is found to weakly increase with the velocity dispersion (mass) of the galaxy.

(ii) The scenario sketched above is perhaps confirmed by the many observational studies based on different kinds of diagnostic. For instance Kuntschner & Davies (1998); Kuntschner (2000, 2001) studying early-type galaxies in the Fornax cluster find the EGs to be coeval and to have metallicities varying from solar to three times solar, and the ratio \[Mg/Fe\] –inferred from the indices $`\mathrm{Mg}_2`$ and Fe3– to weakly correlate with the velocity dispersion. In contrast the data for nearby field galaxies by González (1993); Bressan et al. (1996), and Trager et al. (2000a, b), seem to indicate a much broad range of ages as perhaps compatible with our model galaxies of low initial density. Finally, we call attention on the recent result by Poggianti et al. (2001) for dwarf galaxies confirming the suggestion made by Bressan et al. (1996) and Trager et al. (2000a, b) of a recent star formation in these systems.

(iii) The $`\mathrm{Mg}_2`$ vs $`\sigma `$ relation mainly reflects the metallicity-$`\sigma `$-mass sequence rather than enhancement in $`\alpha `$-elements and/or age.

(iv) Although the dynamical models and companion chemical models we have discussed do not fully match the high metallicities inferred from other diagnostics and some improvement is required, they however seem to offer a coherent scenario accounting for many observed properties of EGs, thus perhaps indicating that we are on the right track.

###### Acknowledgements.

We like to thank Dr. F. Ferrini for several constructive suggestions in his referee report. C.C. likes to acknowledge the friendly hospitality and stimulating environment provided by ESO in Garching where this paper has been written during leave of absence from the Astronomy Department of the Padua University. This study has been financed by the Italian Ministry of Education, University, and Research (MIUR).


text:
# Dynamics of folding in Semiflexible filaments

## Abstract

We investigate the dynamics of a single semiflexible filament, under the action of a compressing force, using numerical simulations and scaling arguments. The force is applied along the end to end vector at one extremity of the filament, while the other end is held fixed. We find that, unlike in elastic rods the filament folds asymmetrically with a folding length which depends only on the bending stiffness $`\kappa `$ and the applied force $`f`$. It is shown that this behavior can be attributed to the exponentially falling tension profile in the filament. While the folding time $`\tau _0`$ depends on the initial configuration, the distance moved by the terminal point of the filament and the length of the fold scales as $`\tau ^{1/2}`$ at time $`\tau >>\tau _0`$ and is independent of the initial configuration.

Recent developments in micro manipulation techniques using optical tweezers, micro pipette and high resolution optical imaging have made the study of single filament dynamics experimentally feasible . Today it is possible to apply forces of the order of pico Newton with precision. This has spurred lot of activity in the study of single bio-polymers. These techniques have been used extensively in the study of stretching elasticity and force extension relation of bio-polymers . Some of these techniques can also be applied to study the effect of compressing forces on single polymer filaments .

In a physical description of their statistical properties, the chemical details of polymers can be summarized using a set of parameters. These parameters characterize the individual features of each polymer species. Two of the most important parameters are persistence length and bending stiffness. Persistence length, $`l_p`$, is the distance along the contour of the filament over which the tangent-tangent correlation is destroyed by thermal fluctuations. Semiflexible polymers have contour length comparable to its persistence length. This means that unlike in the case of flexible filaments or rigid rods the dynamical properties of semiflexible filaments are determined by both conformational entropy and bending stiffness. Thermal fluctuations also play an important role in determining the viscoelastic behavior of semidilute solutions of semiflexible filaments .

From the biological point of view, a study of the mechanical and statistical mechanical properties of semiflexible filaments are important to understand the dynamical behavior of cytoskeletal network. Cytoskeleton, which controls the shape recovery and dynamic cellular reorganization of living cells is a complex network of polymers . At the length scale of a cell most of the polymers which build up this network are semiflexible . Microtubules, $`l_p=5.2mm`$ , and Actin filaments, $`l_p=17\mu m`$  are some of the examples.

For length scales much smaller than the persistence length the filaments are like rigid rods. Under the action of longitudinal compressing forces these filaments exhibit the classical Euler buckling . The dynamics of this buckling has been studied recently . In the case of bio-polymers like microtubules, such buckling is also observed as a result of the force generated by polymerization  and by the movement, on them, of kinesin motors . On the other hand in flexible polymers with persistence length much smaller than the contour length, the compression forces applied at the extremities do not propagate along the polymer. This naturally raises the question, what happens when we apply a longitudinal compressing force on a semiflexible filament? Here we investigate the dynamics of a single semiflexible filament under the action of an external compressing force.

In this paper we look at the folding of a semiflexible filament under the action of a compressing force acting along the end to end vector. Previous calculations of folding in semiflexible polymers assumes uniform tension in the filament . We show that, at finite temperature, the tension profile in a semiflexible filament in equilibrium decreases exponentially along the contour. This tension profile is shown to have non trivial effects on the folding dynamics of filaments. The filament is described by the position vector $`𝐫(s)`$, parameterized by the contour length $`s`$ of the filament, where $`0sL`$. Variety of boundary conditions are possible. Throughout this paper we will discuss only the case wherein the curvature of the filament at the two extreme ends are zero (hinged boundary). This implies free rotation of the tangent vector at the ends. The $`s=0`$ tip of the filament is allowed to move only along the direction of the end to end vector which is taken to be along the $`x`$axis. The position of the $`s=L`$ tip, $`𝐫(L)`$, is fixed. An external force $`𝐟`$ is applied at the $`s=0`$ tip along the $`x`$direction. A detailed discussion of other possible boundary conditions will be published elsewhere . The elastic coefficients for the transverse and longitudinal fluctuations of semiflexible filaments are not the same. This implies that, for the above boundary conditions, a compressing force could lead to either folding or buckling. However, the transverse spring coefficient in the stiff limit is proportional to stiffness $`\kappa `$ and the longitudinal spring coefficient proportional to $`\kappa ^2/T`$ . Thus an applied force on a filament with hinged ends will always lead to folding.

Actin filaments exhibit a persistence length of about $`17\mu m`$ at room temperature . Critical force for buckling is then of the order of a few pico Newton. In a laser tweezers experiment we can apply forces of the order of pico Newton with precision. Typically, in such experiments the filament is held by a bead of size $`0.5\mu m`$ to a precision of $`\pm 50nm`$ . In view of these, it should be possible to observe folding dynamics in laser tweezers experiments using actin filaments of length $`<10\mu m`$ or microtubules of length $`>100\mu m`$ .

In our numerical simulations, semiflexible filament are modeled as $`N`$ rods of equal length hinged together linearly. Position of each hinge is defined by the Cartesian coordinates $`𝐫_i=(x_i,y_i)`$. The total energy of the filament, with the compressing force $`𝐟`$ acting at $`𝐫_1`$, is ,

$$F=\frac{\kappa }{2}\underset{i=1}{\overset{N2}{}}(1𝐭_i𝐭_{i+1})+\underset{i=1}{\overset{N1}{}}g_i𝐭_i^2x_1\widehat{𝐞}_𝐱𝐟,$$
(1)
where $`𝐭_𝐢`$ is the unit vector pointing from $`i`$ to $`i+1`$. The first term is the contribution from the curvature of the polymer, $`\kappa `$ being the elastic constant, and the second term is to enforce the local inextensibility condition, $`𝐭_{𝐢}^{}{}_{}{}^{2}=1`$, of the rods. We work in the over damped limit where the dynamics of the filament is described by the Langevin equations:

$`\mathrm{\Gamma }^1{\displaystyle \frac{x_i}{t}}`$ $`=`$ $`{\displaystyle \frac{F}{x_i}}+\eta _{x_i}(t)`$ (2)
$`\mathrm{\Gamma }^1{\displaystyle \frac{y_i}{t}}`$ $`=`$ $`{\displaystyle \frac{F}{y_i}}+\eta _{y_i}(t).`$ (3)
Here the effect of the solvent enters only through the inverse friction coefficient $`\mathrm{\Gamma }`$. Equations 3 are solved along with the constraint $`𝐭_i^2=1`$. The noise $`\eta `$ obeys the fluctuation dissipation theorem which in the case of filaments with equal monomer mass leads to

$$\eta _{x_i}(t)\eta _{y_j}(t^{})=2\mathrm{\Gamma }^1\delta _{xy}\delta _{ij}\delta (tt^{}).$$

Integrating equations 3 we get the increments $`\mathrm{\Delta }x_i`$ and $`\mathrm{\Delta }y_i`$, which are functions of the unknown Lagrange multipliers $`g_i`$. Increment in time is represented by the dimensionless variable $`\mathrm{\Delta }\tau =\mathrm{\Gamma }\mathrm{\Delta }t`$ . The coordinates of the vertices at time $`\tau +\mathrm{\Delta }\tau `$ obtained from the above equations have to satisfy the constraint $`𝐭_i^2(\tau +\mathrm{\Delta }\tau )=1`$. These $`N1`$ simultaneous nonlinear equations are solved to determine the Lagrange multipliers $`g_i`$. We choose $`\mathrm{\Delta }x_N=0`$, $`\mathrm{\Delta }y_N=0`$ and $`\mathrm{\Delta }y_0=0`$ at all times. There is no constraint on the angle between the filament and the force direction at both ends.

We have carried out extensive numerical simulation of the above model using a chain with length $`N=150`$ for different values of $`\kappa `$ and $`f`$. One run in the simulation correspond to $`\mathrm{\Delta }\tau =\mathrm{\hspace{0.17em}1}`$. All results are averaged over 10 initial configurations. The initial conformation of the filament $`(\tau =0)`$ is chosen from an equilibrium distribution which satisfy the boundary condition $`y_1=0`$ and $`y_N=0`$ . The snapshots from our simulation for a constant applied force $`f`$ is shown in figure $`1`$. We observe that the filament first develop a fold at a distance $`l_f`$ from the end at which the force is applied. $`l_f`$ is the position of the first maximum in curvature along the contour when the free end is perpendicular to the end to end vector. For a given $`\kappa `$ and force $`f`$, this ”folding length”$`(l_f)`$ is independent of the initial configuration and is found to depend on bending stiffness as $`l_f\sqrt{\kappa }`$. This is shown in figure $`2`$.

FIG. $`1`$ The snapshots of conformations of a semiflexible filament with $`\kappa =500`$ and $`N=150`$, taken at equal intervals of time $`\mathrm{\Delta }\tau =0.2`$ million for an applied force $`f=30`$.

FIG.$`2`$ A plot of $`\mathrm{log}_{10}\kappa `$ vs $`\mathrm{log}_{10}l_f`$. The filled squares are from the simulation. The dotted line is a fit which shows $`l_f\sqrt{\kappa }`$.

A better understanding of the dependence of folding length on $`\kappa `$ is obtained from the continuum equations for the worm like chain model  . The equation of motion in the over damped limit for the boundary conditions, $`𝐫(L)=`$a constant and $`y(0)=0`$, is then

$`{\displaystyle \frac{1}{\mathrm{\Gamma }}}{\displaystyle \frac{𝐫}{t}}`$ $`=`$ $`{\displaystyle \frac{\delta F}{\delta 𝐫}}=l_p_s^4𝐫+2_s(g\dot{𝐫})`$ (5)
$`+\left(f+2g(s)\dot{𝐫}(s)l_p_s^3𝐫\right)\delta (s)+\eta (s,t),`$
here dots denote the derivatives with respect to arc length. $`\eta (s,t)`$ is the thermal noise and $`\mathrm{\Gamma }`$ is a kinetic coefficient. Using the constraint $`\dot{𝐫}(s)^2=1`$, and equation 5 we get,

$$\frac{1}{2\mathrm{\Gamma }}\frac{\dot{𝐫}^2}{t}=2\ddot{g}2g\ddot{r}^2l_p\dot{𝐫}_s^5𝐫=0.$$
(6)
Comparing the curvature and force terms in equation 5, we get a length which scales as $`l\sqrt{(}\kappa /f)`$. Similar length scale is also obtained from the constraint equation.

Let us try to get some insight into the tension profile along the filament. Dropping higher order derivative terms from the constraint equation 6 leads to

$$\frac{1}{2\mathrm{\Gamma }}\frac{\dot{𝐫}^2}{t}=\ddot{g}^fg^f\ddot{r}^2=0.$$
(7)
We use $`g^f`$ to distinguish this tension from that obtained with $`\kappa `$. Numerical solution of the above equation for time $`t=0`$ is shown in figure $`3`$. Here we have replaced $`\ddot{𝐫}(s)^2`$ by its approximate ensemble averaged value

$$\ddot{r}(s)^2=\left(\frac{1}{l_pL}\right)\frac{_n\mathrm{sin}^2\left(\frac{n\pi s}{L}\right)}{1\frac{L}{\kappa \pi }_n\frac{1}{n^2}\mathrm{sin}^2\left(\frac{n\pi s}{L}\right)}$$
(8)
FIG.$`3`$ $`g(s)`$ obtained from equation 7 for $`\kappa =500,1000,1500,2000,2500`$ and n=150.

We see that the resulting tension decays exponentially. We find the length scale to vary as $`\sqrt{\kappa }`$. The tension profile is found to depend on the number of modes taken in equation 8 and is exponential only in the limit when the number of modes are large. The lower cut off in length being set by the segment length.

We now compare this tension profile with that obtained from the simulations. In the simulations, the tension $`g_i`$ has contributions from three sources, namely, curvature, noise and the applied force. Like in the case of continuum equations we can calculate a tension $`g_i^f`$ by setting the noise and curvature terms to zero. Since the equations of motion are nonlinear one may doubt the validity of a comparison between the tension $`g_i^f`$ obtained from this calculation with $`g^f(s)`$ obtained from the continuum equations. To dispel this doubt we linearize the equations for the tension profile by neglecting terms of the order $`(\mathrm{\Delta }x_i)^2`$ and $`(\mathrm{\Delta }y_i)^2`$, in equation $`𝐭_i^2=1`$. The results obtained from such a linear calculation agree well with the full non-linear results (data not shown). However, this linearization also means deviation from the constraint $`𝐭_i^2=1`$. Even for small $`\mathrm{\Delta }\tau `$ these errors in segment length can become very large in a short time. Therefore we have to use the full non-linear equations for updating the filament configuration.

In figure $`4`$ we show this tension profile $`g_i^f`$ as a function of the vertex index $`i`$. The tension $`g_i^f`$ goes to zero exponentially as predicted earlier. The distance $`\xi `$ over which the tension falls to half its maximum value is found to increase as $`\sqrt{\kappa }`$ for a fixed value of $`f`$.

FIG.4 $`g^f/f`$ as a function of the contour length $`i`$ from the point at which force is applied, for $`\kappa =300,1000,2000,3000`$. Inset: $`\mathrm{log}(\kappa )`$ is plotted against the $`\mathrm{log}\xi `$, showing $`\xi \kappa ^{1/2}`$

We now look at how the buckling length $`l_f`$ depends on the applied force $`f`$. We have seen in the previous paragraph that the tension in a semiflexible filament falls exponentially along its contour. This exponential tension profile differentiates between semiflexible filaments and stiff elastic rods wherein the tension profile is uniform. Thus it is instructive to extend the stability calculations for a rod with uniform tension  to the present case of exponential tension profile of the form $`\mathrm{exp}(\alpha s/\sqrt{\kappa })+1`$. The equation of equilibrium is,

$$\kappa \frac{^4y}{s^4}+\frac{f}{2}\frac{}{s}\left[\frac{y}{s}\left[\mathrm{exp}(\alpha s/\sqrt{\kappa })+1\right]\right]=0.$$
(9)
These equations are solved numerically with the boundary conditions $`\ddot{y}|_{s=L}=0`$, $`\ddot{y}|_{s=0}=0`$, $`y(0)=0`$, $`y(L)=0`$. The lowest energy solutions of this equation show that the value $`s=s_{max}`$, of the first maximum in curvature, shifts to smaller values of $`s`$ with increasing $`f`$. This is indeed what we see in our simulations of the semiflexible filament.

FIG.$`5`$ Folding length ($`l_f`$) as a function of the applied force($`f`$); The dotted line represents the buckling point obtained from the numerical solution of equation 9 with $`\alpha =10`$. The straight line corresponds to $`l_f\frac{1}{\sqrt{f}}`$

In figure $`5`$ we compare the buckling length $`l_f`$ obtained from the simulation with $`s_{max}`$ from the solution of equation 9. $`\alpha `$ is the only adjustable parameter used to get the agreement between the two results.

We now come to the dynamical behavior of the fold. The length of the segment before the fold, the end to end distance and tension $`g_i^f`$ all depend on time. As can be seen from the configuration snap shots in figure $`1`$, once the filament has folded the applied force acts to pull the filament. This implies a change in sign of tension $`g_i^f`$. At $`\tau =0`$ we have $`g_1^f=f`$ . As the filament folds $`g_1^f`$ changes sign saturating to $`g_1^f=f`$ . The time $`\tau _0`$ at which $`g_1^f(\tau _0)=0`$ can be designated as the folding time. We find this value of the folding time to depend on the initial configuration of the filament. The end of the filament at which the force is applied, moves by a distance $`x_1(\tau )`$ at time $`\tau `$. This is depicted in figure $`6`$ along with the length of the fold.

For $`\tau >>\tau _0`$ all the distances scale as $`\tau ^{1/2}`$. This shows that in this regime the problem is similar to that of stretching a tense string . For large $`f`$ we can then neglect the contribution from the curvature energy and noise. The dynamics is then determined by the straight segment up to the fold moving through the fluid. This results in a linear tension profile along the filament , with the tension falling to zero approximately at the fold. At $`s=0,g^f(s=0)=f`$, by neglecting the curvature and noise terms on the right hand side of the equation 5, using dimensional analysis, we get, $`x(s=0,\tau )(\tau )^{1/2}`$. This is shown in figure $`6`$.

FIG. $`6`$ The displacement of $`x_1(\tau )`$ after folding of the filament (Filled squares). Crosses shows the position of the fold along the contour. Open circles represent the arc length up to which $`g^f`$ is linear. $`(\tau \tau _0)^{1/2}`$ line (solid line) is given as a guide to the eye

As the end moves, the position of the fold moves along the contour ( see figure $`1`$). Hence the fold length and the point at which tension falls to zero too follow similar behavior as that of $`x(s=0,\tau )`$.

To conclude, we have studied the dynamics of a single semiflexible filament under the action of a compressing force at one end which lead to folding. We have shown that the folding length($`l_f`$) of the filament is proportional to $`\sqrt{\kappa }`$ and decreases monotonically with the applied force $`f`$. This follows from the fact that the tension in semiflexible filaments decay exponentially along the contour. It is shown that buckling of a homogeneous elastic rod with an exponential tension profile is qualitatively similar to that of the folding of a semiflexible filament. The folding length for a given $`\kappa `$ and force $`𝐟`$ does depend on the strength of the noise $`\eta `$. Thus once the initial configuration is chosen, a ”zero temperature” calculation should give the same results. As the filament folds the tension changes sign and we move into the stretching regime. We have shown, using numerical simulations and scaling arguments, that in this regime the fold length, the end point $`x_1(\tau )`$ and the distance along the contour to which tension has propagated, scales as $`\tau ^{1/2}`$.

Acknowledgments
We would like to thank S. K Srivatsa and Madan Rao for many discussions and comments. P.B.S thanks U. Seifert for discussions during a stay at MPIKG, Golm. We thank CSIR, India for financial support.