text:
# Adsorption of Oppositely Charged Polyelectrolytes onto a Charged Rod

## I Introduction

The adsorption of charged polymer chains \[polyelectrolytes (PEs)\] on oppositely charged surfaces and particles has been the subject of intense theoretical and experimental research work. Due to the strong mutual electrostatic attraction, the resulting complex is highly stable and can exhibit new and interesting physico-chemical properties. For instance, such a process is employed to stabilize charged colloidal suspensions. Another more elaborated PE adsorption technique, consists of successive deposition of polycations (PAs) and polyanions (PAs) leading to well controlled PE multilayer structures. Such materials present enormous technological applications and this deposition method has been especially exploited with planar Decher\_1997  and spherical Caruso\_Science\_1998  substrates. In principle this technique could also be applied to (ultra-thin) cylindrical substrates. However due to the very low dimensionality there, a quantitative experimental characterization of the PE multilayering is more difficult.

On the theoretical side, there exist a few analytical works about PE multilayers on charged planar surfaces,Solis\_JCP\_1999 ; Netz\_Macromol\_1999b ; Castelnovo\_2000  but no theory is available for thin cylindrical substrates. Recent MC simulations were devoted for the PE multilayering on planar Messina\_Macromol\_2003  and spherical Messina\_Langmuir\_2003  charged surfaces. The aim of this paper is to fill this gap and to investigate by means of MC simulations the structure of PCs and PAs adsorbed onto a charged cylindrical surface. The relevant role of the entropy is especially pointed out.

Our paper is organized as follows: Sec. II is devoted to the description of our MC simulation model. The measured quantities are specified in Sec. III. The PC adsorption is studied in Sec. IV. Then the adsorption of both PC and PA chains is addressed in Sec. V. Finally, Sec. VI provides brief concluding remarks.

## II Simulation model

The setup of the system under consideration is similar to those recently investigated with spherical Messina\_Langmuir\_2003  and planar Messina\_Macromol\_2003  substrates. Within the framework of the primitive model we consider a PE solution near a charged hard rod with an implicit solvent (water) of relative dielectric permittivity $`ϵ_r80`$. This fixed charged rod (whose axis is located at $`r=\sqrt{x^2+y^2}=0`$) is characterized by a negative linear bare charge density $`\lambda _0e`$, where $`e`$ is the (positive) elementary charge and $`\lambda _0>0`$ is the number of charges per unit length. Electroneutrality is always ensured by the presence of explicit monovalent ($`Z_c=1`$) rod’s counterions of diameter $`a`$. PE chains ($`N_+`$ PCs and $`N_{}`$ PAs) are made up of $`N_m`$ monovalent monomers ($`Z_m=1`$) of diameter $`a`$. Hence, all microions are monovalent: $`Z=Z_c=Z_m=1`$ and have the same diameter $`a`$. For the sake of simplicity, we only consider here symmetrical complexes where PC and PA chains have the same length and carry the same charge in absolute value. Here the radius $`r_{rod}`$ of the charged rod is also $`a/2`$. To keep further simplicity in our model, we also assume the same linear charge density for the rod $`\lambda _0e=1e/a`$ as the flexible polyelectrolytes. Hence the only differences between the rod and the PE chains are the chain stiffness and and the chain length, where both of them are infinite for the rod.

All these particles making up the system are confined in a $`R\times L`$ cylindrical box. Periodic conditions are applied in the $`z`$ direction, whereas hard coaxial cylinders are present at $`r=a/2`$ (location of the cylindrical charged wall) and at $`r=R`$ (location of an uncharged cylindrical wall). To avoid the appearance of image charges,Fixman\_Macromol\_1978 ; Messina\_image\_2002  we assume that on both parts of the charged rod (at $`r=a/2`$) the dielectric constants are the same.

The total energy of interaction of the system can be written as

$`U_{tot}`$ $`=`$ $`{\displaystyle \underset{i}{}}\left[U_{hs}^{(rod)}(r_i)+U_{coul}^{(rod)}(r_i)+U_{vdw}^{(rod)}(r_i)\right]+`$
$`{\displaystyle \underset{i,i<j}{}}\left[U_{hs}(R_{ij})+U_{coul}(R_{ij})+U_{FENE}(R_{ij})+U_{LJ}(R_{ij})\right],`$
where the first (single) sum stems from the interaction between an ion $`i`$ (located at $`r=r_i`$) and the charged plate, and the second (double) sum stems from the pair interaction between ions $`i`$ and $`j`$ with $`R_{ij}=|𝐑_i𝐑_j|`$ where $`𝐑_i=x_i𝐞_x+y_i𝐞_y+z_i𝐞_z`$. All these contributions to $`U_{tot}`$ in Eq. (II) are described in detail below.

Excluded volume interactions are modeled via a hardcore potential defined as follows

$$U_{hs}(R_{ij})=\{\begin{array}{cc}\mathrm{},\hfill & \mathrm{for}R_{ij}<a\hfill \\ 0,\hfill & \mathrm{for}R_{ij}a\hfill \end{array}$$
(2)
for the microion-microion one, except for the monomer-monomer one,note\_HS  and

$$U_{hs}^{(rod)}(z_i)=\{\begin{array}{cc}\mathrm{},\hfill & \mathrm{for}r_i<r_{rod}+a/2\hfill \\ \mathrm{},\hfill & \mathrm{for}r_i>Ra/2\hfill \\ 0,\hfill & \mathrm{for}r_{rod}+a/2r_iRa/2\hfill \end{array}$$
(3)
for the rod-microion one. Hence the minimal radius, $`r_0`$, of closest approach between the rod-axis and a microion is $`r_0=r_{rod}+a/2=a`$.

The electrostatic energy of interaction between two ions $`i`$ and $`j`$ reads

$$\frac{U_{coul}(R_{ij})}{k_BT}=\pm \frac{l_B}{R_{ij}},$$
(4)
where +(-) applies to charges of the same (opposite) sign, and $`l_B=e^2/(4\pi ϵ_0ϵ_rk_BT)`$ is the Bjerrum length corresponding to the distance at which two monovalent ions interact with $`k_BT`$. The electrostatic energy of interaction between an ion $`i`$ and the (uniformly) charged rod reads

$$\frac{U_{coul}^{(rod)}(r_i)}{k_BT}=\pm 2l_B\lambda _0\mathrm{ln}(r_i/r_{rod})$$
(5)
where +(-) applies to positively (negatively) charged ions. An appropriate and efficient modified Lekner sum was utilized to compute the electrostatic interactions with periodicity in one direction. Brodka\_MolPhys\_2002b  To link our simulation parameters to experimental units and room temperature ($`T=298`$K) we choose $`a=4.25`$ Å leading to the Bjerrum length of water $`l_B=1.68a=7.14`$ Å and to the so-called Manning condensation parameter Manning\_JCP\_1969 ; Manning\_BBPC\_1996  $`\xi _M=l_B/a=1.68`$.

The polyelectrolyte chain connectivity is modeled by employing a standard FENE potential in good solvent (see, e.g., Ref. Kremer\_FENE\_1993 ), which reads

$$U_{FENE}(R_{ij})=\{\begin{array}{cc}\frac{1}{2}\kappa R_0^2\mathrm{ln}\left[1\frac{R_{ij}^2}{R_0^2}\right],\hfill & \text{for}R_{ij}<R_0\hfill \\ & \\ \mathrm{},\hfill & \text{for}R_{ij}R_0\hfill \end{array}$$
(6)
with $`\kappa =27k_BT/a^2`$, $`R_0=1.5a`$, and $`i`$ and $`j`$ stand for two consecutive monomers along a chain. The excluded volume interaction between chain monomers is taken into account via a purely repulsive Lennard-Jones (LJ) potential given by

$$U_{LJ}(R_{ij})=\{\begin{array}{cc}4ϵ\left[\left(\frac{a}{R_{ij}}\right)^{12}\left(\frac{a}{R_{ij}}\right)^6\right]+ϵ,\hfill & \text{for}R_{ij}2^{1/6}a\hfill \\ & \\ 0,\hfill & \text{for}R_{ij}>2^{1/6}a\hfill \end{array}$$
(7)
where $`ϵ=k_BT`$. These parameter values lead to an equilibrium bond length $`l=0.98a`$.

An important interaction in PE multilayering is the non-electrostatic short ranged attraction, $`U_{vdw}^{(rod)}`$, between the cylindrical macroion and the PC chain. To include this kind of interaction, we choose without loss of generality a (microscopic site-site) van der Waals (VDW) potential of interaction between the cylindrical macroion and a PC monomer that is given by

$$U_{vdw}^{(rod)}(r_i)=ϵ\chi _{vdw}\left(\frac{a}{r_ir_0+a}\right)^6\text{for}r_ir_0,$$
(8)
where $`\chi _{vdw}`$ is a positive dimensionless parameter describing the strength of this attraction. Thereby, at contact (i.e., $`r=r_0`$), the magnitude of the attraction is $`\chi _{vdw}ϵ=\chi _{vdw}k_BT`$ which is, in fact, the relevant characteristic of this potential.note\_VDW  Since it is not straightforward to directly link this strength of adsorption to experimental values, we choose $`\chi _{vdw}=5`$ as considered previously.Messina\_Langmuir\_2003 ; Messina\_Macromol\_2003

All the simulation parameters are gathered in Table 1. The set of simulated systems can be found in Table 2. The equilibrium properties of our model system were obtained by using standard canonical MC simulations following the Metropolis scheme.Metropolis\_JCP\_1953 ; Allen\_book\_1987  Single-particle moves were considered with an acceptance ratio of $`30\%`$ for the monomers and $`50\%`$ for the counterions. Typically, about $`10^5`$ to $`10^6`$ MC steps per particle were required for equilibration, and no less than $`5\times 10^5`$ subsequent MC steps were used to perform measurements. To improve the computational efficiency, we omitted the presence of PE counterions when $`N_+=N_{}`$ so that the system is still globally electroneutral. We have systematically checked for $`N_+=N_{}=6`$ (system $`B`$) that the (average) PE configurations (especially the monomer distribution) are qualitatively the same of those where PE counterions are explicitly taken into account, as it should be.

## III Measured quantities

We briefly describe the different observables that are going to be measured. In order to characterize the PE adsorption, we compute the monomer density $`n_\pm (r)`$ that is normalized as follows

$$_{r_0}^{Ra/2}n_\pm (r)2\pi rL𝑑r=N_\pm N_m$$
(9)
where $`+()`$ applies to PCs (PAs). This quantity is of special interest to characterize the degree of ordering in the vicinity of the cylindrical macroion surface.

The total number of accumulated monomers $`\overline{N}_\pm (r)`$ within a distance $`r`$ from the axis (located at $`r=0`$) of the cylindrical macroion is given by

$$\overline{N}_\pm (r)=_{r_0}^rn_\pm (r^{})2\pi r^{}L𝑑r^{}$$
(10)
where $`+()`$ applies to PCs (PAs). This observable will be addressed in the study of PC adsorption. (Sec. IV).

Another relevant quantity is the global net fluid charge $`\lambda (r)`$ which is defined as follows

$$\lambda (r)=_{r_0}^r\left[\stackrel{~}{n}_+(r^{})\stackrel{~}{n}_{}(r^{})\right]2\pi r^{}𝑑r^{},$$
(11)
where $`\stackrel{~}{n}_+`$ ($`\stackrel{~}{n}_{}`$) stand for the density of all the positive (negative) microions (i.e., monomers and counterions). Thus, $`\lambda (r)`$ corresponds to the (reduced) net fluid charge per unit length (omitting the bare macroion linear-charge $`\lambda _0`$) within a distance $`r`$ from the axis of the charged cylindrical wall. At the uncharged wall, electroneutrality imposes $`\lambda (r=Ra/2)=\lambda _0`$. By simple application of the Gauss’ law, $`\left[\lambda (r)\lambda _0\right]/r`$ is directly proportional to the mean electric field at $`r`$. Therefore $`\lambda (r)`$ can measure the screening strength of the cylindrical macroion-charge by the neighboring solute charged species.

## IV Adsorption of polycations

In this part, we study the adsorption of PC chains (system $`A`$) for the two couplings $`\chi _{vdw}=0`$ and $`\chi _{vdw}=5`$. The case $`\chi _{vdw}=0`$, corresponding to a purely electrostatic regime, is well understood for ”normal spherical” electrolyte ions near a charged rod. Manning\_JCP\_1969 ; Belloni\_CollSurf\_1998 ; Deserno\_Macromol\_2000  Nevertheless, the situation is much more delicate for PE adsorption, and that has not been the subject of numerical simulations.note\_dias  This is a decisive step to elucidate the even more complex PE multilayer structures where additionally PAs are also present.

Here, where $`N_{}=0`$ (i.e., no polyanions), global electroneutrality is ensured by the presence of explicit PC’s counterions (i.e., monovalent anions) and the macroion-rod’s counterions (i.e., monovalent cations).

The profiles of the monomer density $`n_+(r)`$ are depicted in Fig. 1. Near contact (i.e., $`rr_00`$), the density $`n_+(r)`$ at $`\chi _{vdw}=5`$ is about one order of magnitude larger than at $`\chi _{vdw}=0`$. To further estimate the degree of adsorption of PC monomers we have plotted the fraction $`\overline{N}_+(z)/(N_+N_m)`$ of adsorbed monomers in Fig. 2. At a radial distance $`rr_0=a`$ (corresponding to a width of two monomers), almost all the the monomers (more than 99%) are adsorbed for $`\chi _{vdw}=5`$ against only $`68\%`$ for $`\chi _{vdw}=0`$. This strong observed PC adsorption at $`\chi _{vdw}=0`$ is of course due to the high valence of the polycation (here $`N_m=20`$) and the intense attractive electric field near the thin rod. On the other hand, at high enough $`\chi _{vdw}`$ and (even) without rod-monomer electrostatic interaction \[i.e., for a neutral rod ($`\lambda _0=0`$)\], the PC adsorption would be very similar to that obtained at finite $`\lambda _0`$.

The (global) net fluid charge $`\lambda (r)`$ is reported in Fig. 3. In all cases we detect a macroion-surface charge reversal (i.e., $`\lambda (r)/\lambda _0>1`$). The position $`r=r^{}`$ at which $`\lambda (r^{})`$ gets its maximal value decreases with $`\chi _{vdw}`$, due to the $`\chi _{vdw}`$-enhanced adsorption of the PCs. Concomitantly, this overcharging increases with $`\chi _{vdw}`$, since the (extra) gain in energy by macroion-monomer VDW interactions can better overcome (the higher $`\chi _{vdw}`$) the cost of the self-energy stemming from the adsorbed excess charge. Messina\_Langmuir\_2003 ; Messina\_Macromol\_2003  More quantitatively, we have $`\lambda (r^{}r_01.1a)/\lambda _02.96`$ at $`\chi _{vdw}=5`$ against $`\lambda (r^{}r_03.4a)/\lambda _02.67`$ at $`\chi _{vdw}=0`$. Note that the maximal value of charge reversal of $`(6\times 2040)/40=200\%`$ (i.e., $`\lambda (r^{})/\lambda _0=3`$) allowed by the total charge of PCs can not be exactly reached due to a slight accumulation of microanions. However at $`\chi _{vdw}=5`$, the presence of microanions near the rod is very marginal, leading in practice to an overcharging that is fully governed by the PCs. Again this giant overcharging found at any $`\chi _{vdw}`$ is due to the high (intrinsic) electric field generated by the macroion rod. Indeed, for the same chain parameter, it was found in planar geometry (with a surface charge density $`\sigma _0e=0.165\text{C/m}^2`$) that only a charge reversal of about 70% \[i. e., $`\sigma (r^{})/\sigma _0=1.7`$\] could be reached at $`\chi _{vdw}=5`$.Messina\_Macromol\_2003  For cylindrical substrates, it turns out that the equivalent surface charge density \[here $`\lambda _0e\times L/(2\pi r_{rod}L)0.282\text{C/m}^2`$\] can be much higher in experimental conditions. Those (locally) overcharged states should be the driving force for the building of subsequent PE bilayers when PA chains are added.

Typical equilibrium configurations can be seen in Fig. 4. The qualitative difference between $`\chi _{vdw}=0`$ \[Fig. 4(a)\] and $`\chi _{vdw}=5`$ \[Fig. 4(b)\] is remarkable. Without additional VDW attraction ($`\chi _{vdw}=0`$) the adsorption is more diffuse than at $`\chi _{vdw}=5`$, where in the latter situation the $`r`$-fluctuation is very weak. Interestingly, one can observe that the PC chains tend to adopt a helical structure. This kind of conformation is the non-trivial result of two antagonistic forces: (i) electrostatic correlations and (ii) entropy. To clarify this point, one has to consider two typical length scales:

* If the cylindrical macroion is larger or about the PE chain (i.e., $`LN_ml`$ as it is presently the case), then the ground state note\_helix  corresponds to straight PE chains sticked to the surface of the rod. At finite temperature one obtains PC chains spiraling around as the result of the coupled effects of entropy and electrostatic correlations.
* If $`LN_ml`$, then the ground state (already) corresponds to perfectly ordered helical structures (so as to overscreen the rod) which locally resemble a one-dimensional Wigner crystal. Nguyen\_JCP\_2000  At finite temperature and strong Coulomb coupling (as it is the case with highly charged PEs and macroion-substrates), this high helical ordering persists at shorter range.

Those reasonings above should at least hold when the equivalent surface charge densities of the PE and that of the rod are similar, as it is the case in this study.

## V Adsorption of polycations and polyanions

We now consider the systems $`BD`$ where additionally PA chains are present, so that we have a (globally) neutral polyelectrolyte complex (i.e., $`N_+=N_{}=N_{PE}/2`$). Global electroneutrality is ensured by the counterions of the cylindrical macroion as usual. We stress the fact that this process is fully reversible for the parameters under investigation. In particular, we checked that the same final equilibrium configuration is obtained either by (i) starting from system with isolated PCs and then adding PAs or (ii) starting directly with the oppositely charged PEs.

### V.1 Monomer distribution

The profiles of the monomer density $`n_\pm (r)`$ at $`\chi _{vdw}=0`$ and $`\chi _{vdw}=5`$ are depicted in Fig. 5(a) and Fig. 5(b), respectively. The corresponding microstructures are sketched in Fig. 6.

#### V.1.1 Case of zero-$`\chi _{vdw}`$

At $`\chi _{vdw}=0`$, a comparison between systems $`A`$ (see Fig. 1) and $`B`$ indicates that the adsorption of PC monomers is weaker when additional PAs are present. This effect was already observed with spherical Messina\_Langmuir\_2003  and planar Messina\_Macromol\_2003  substrates, and the same mechanisms apply here to cylindrical substrates. More explicitly, the PC tends to build up a globular state (reminiscent of the classical bulk PE collapse) by getting complexed to the PA chain, as perfectly illustrated in Fig. 6(a). Thereby, the mean number of nearest monomer-neighbors gets higher which is both (i) entropically and (ii) energetically (at least from the PE complex viewpoint) favorable. It turns out that for macroion-rods with large curvature, this PC desorption is counter-intuitively much higher than in the planar case, Messina\_Macromol\_2003  although a much weaker charge reversal (without PAs) was found for the latter. Indeed, we just saw that the overcharging by PCs alone can be quite huge (of the order of 200 % even at $`\chi _{vdw}=0`$ \- see Fig. 3), and thus, one would naively expect a new wave of charge reversal in system $`B`$ (where $`N_+=N_{}=6`$) of the same order (see also Sec. V.2 for a more quantitative description) which in turn should lead to a (very) stable bilayering. But Figs. 5(a) and 6(a) demonstrate that it is not the case. Furthermore, an inspection of the density profile $`n_{}(r)`$ of the PA monomers reveals that there is no peak for $`N_{PE}=12`$. Those interesting findings allow us to already draw a partial conclusion:

* At large enough rod-curvature and small enough PE concentration, even if a giant overcharging would occur in the absence of PAs, it is not sufficient to promote bilayering.

Hence, this spectacular PA-induced PC desorption is due to entropic effects. More precisely, the adsorption onto quasi one-dimensional objects is extremely costly from an entropy viewpoint, and if the number of available dipoles (i.e., PC-PA ion pairs) is too small then ion-dipole correlations are not strong enough to overcome the “entropic barrier”.

Upon increasing $`N_{PE}`$ (i.e., the PE concentration) one sees that the density $`n_+(r)`$ at contact increases \[also identifiable on Figs. 6(a), (c) and (e)\] and gets gradually less sensitive to $`N_{PE}`$, showing that a ”steady” regime should be reached at large enough $`N_{PE}`$. Concomitantly, by increasing $`N_{PE}`$, Fig. 5(a) shows that the adsorption of PA monomers increases \[also identifiable on Figs. 6(a), (c) and (e)\]. This $`N_{PE}`$-enhanced adsorption of monomers is essentially governed by two concomitant mechanisms:

* At higher $`N_{PE}`$, the global ion-dipole correlation gets enhanced (due to an increase of PC-PA ion pair concentration) which favors a stronger monomer adsorption.
* Concomitantly, at higher $`N_{PE}`$, where the available volume becomes lower, the entropy loss undergone upon adsorption gets reduced, which also favors a stronger monomer adsorption.

Nevertheless, only a very broadened peak in $`n_{}(r)`$ (at $`N_{PE}=48`$) is seen and no second peak in $`n_+(r)`$ appears, indicating a weakly stable bilayer. This feature strongly contrasts with our previous findings for planar substrates Messina\_Macromol\_2003  where stable trilayers were reported at $`\chi _{vdw}=0`$ (also with $`N_m=20`$). Therefore, we can state:

* PE multilayering in the electrostatic regime is highly inhibited at large rod-curvature, mainly due to entropic effects. It is to say that the (conformational) PE entropy along the rod is much smaller than that in the bulk and also smaller than that on a planar surface.

#### V.1.2 Case of finite $`\chi _{vdw}`$

In general all our mechanisms previously described at $`\chi _{vdw}=0`$ are also present here at $`\chi _{vdw}=5`$, but one has now to take into account the possible dominance of the VDW interaction near the surface of the cylindrical macroion.

At $`\chi _{vdw}=5`$, Fig. 5(b) shows that the first PC layer is always very stable (despite of the presence of PAs) due to the the sufficiently strong rod-PC VDW attraction. In parallel, $`n_+(r)`$ and in general the ordering increase with $`N_{PE}`$, as it should be. At large $`N_{PE}=48`$, one can detect a second broadened peak in $`n_+(r)`$ which could be the signature of the onset of a third PE layer. However, the height of this peak is still significantly smaller than the corresponding decaying PA density $`n_{}(r)`$ in this region. A visual inspection of Fig. 6(b), (d) and (f) confirms all those features. Hence, even with the help of a specific VDW attractive interaction, it seems that the formation of stable structures beyond bilayers is highly inhibited at strong confinement. This again proves how important is the role of entropy in PE multilayering with low dimensional substrates.

### V.2 Charge distribution

The net fluid charge $`\lambda (r)`$ is reported in Fig. 7. In all cases, the cylindrical macroion gets (locally) overcharged. Surprisingly, beyond the PC overcharging layer the fluid charge does not reenter below unity (at least for $`r<6a`$), proving that there are no charge oscillations in the vicinity of the rigid rod. Even with an overcharging of $`\lambda (r^{})/\lambda _0=3.5`$ for $`\chi _{vdw}=5`$ and $`N_{PE}=48`$, no charge oscillation is generated. This is, nevertheless, consistent with our previous results of Fig. 5 where no multilayering was observed. This phenomenon is due to the fact that our PE concentration is still (very) small even at $`N_{PE}=48`$, which is also the experimental situation in the process of PE multilayering where a rinsing is applied at each successive layer deposition Decher\_1997 . It is well known that at high finite electrolyte concentration, charge oscillations take place near a charged cylindrical wall. Marcelo\_JCP\_1985 ; Deserno\_Macromol\_2000 ; Deserno\_JBCB\_2001

Our results clearly show that the degree of overcharging is systematically larger at high $`\chi _{vdw}=5`$, as also observed without PAs (see Fig. 3). The values of $`r^{}`$ are also systematically shifted to smaller ones as one increases $`\chi _{vdw}`$. However, the position of the peak ($`r^{}r_0a`$) for $`\chi _{vdw}=5`$ remains nearly independent of $`N_{PE}`$ in contrast to what happened at $`\chi _{vdw}=0`$, where $`r^{}`$ is systematically shifted to smaller radial distances upon increasing $`N_{PE}`$. All those observations are consistent with our mechanisms presented so far.

Nonetheless, as soon as oppositely charged polyions can interact, there is a subtle interplay between clustering and the lateral correlations of polyions that governs the degree of overcharging near the macroion surface, as it was already demonstrated in Refs. Messina\_Langmuir\_2003 ; Messina\_Macromol\_2003 . Additionally, at strong rod curvature, entropic effects become also considerable. Consequently, the physics of overcharging involved in the adsorption of oppositely charged PEs onto one-dimensional substrates is highly complicated.

## VI Concluding remarks

In summary, we have carried out MC simulations to elucidate the adsorption of flexible highly charged polycations and polyanions onto a charged cylindrical substrate. Within the dilute regime, we considered the influence of PE concentration. In order to enhance the possible formation of PE multilayers we have also considered an extra non-electrostatic short-range attraction (characterized here by $`\chi _{vdw}`$) between the cylindrical macroion surface and the polycations.

As far as the PC adsorption is concerned (in the absence of PAs), we demonstrated that huge macroion charge reversal occurs even in a purely electrostatic regime with $`\chi _{vdw}=0`$. By adding exactly the same amount of PAs, we surprisingly observe a (relatively) marginal overcharging which is due to (i) PC-PA clustering and (ii) above all to entropic effects.

At higher number of PEs, our results show that true bilayering (i.e.; flat and dense PE layers) can only occur at finite $`\chi _{vdw}`$, in contrast to what was recently found with planar substrates. Messina\_Macromol\_2003  Even at finite $`\chi _{vdw}`$, our simulation data demonstrate that stable multilayering (beyond bilayering) is hard to reach at large macroion rod-curvature, due to the high entropy loss there. This latter in turn inhibits the appearance of charge oscillations.

However, there must be a certain regime of curvature, typically when the rod radius is larger or about the chain size (i.e., $`r_{rod}N_ml`$), where such entropy effects become less relevant and would then lead to similar behaviors to those occurring with planar substrates. In this respect, a future work could include the effect of rod-curvature on the formation of PE multilayers. We hope that this contribution will trigger new theoretical as well as experimental investigations.

###### Acknowledgements.

The author thanks C. N. Likos and H. Löwen for helpful discussions and a critical reading of the manuscript.


format: text

text:
# Ferromagnetism in Fe-doped Ba6Ge25 Chiral Clathrate

## Abstract

We have successfully synthesized a Ba6Ge25 clathrate, substituting 3 Fe per formula unit by Ge. This chiral clathrate has Ge sites forming a framework of closed cages and helical tunnel networks. Fe atoms randomly occupy these sites, and exhibit high-spin magnetic moments. A ferromagnetic transition is observed with $`T_c`$ = 170 K, the highest observed $`T_c`$ for a magnetic clathrate. However, the magnetic phase is significantly disordered, and exhibits a transformation to a re-entrant spin glass phase. This system has a number of features in common with other dilute magnetic semiconductors.

Si, Ge and Sn clathrates have received increasing attention over the past few years. These materials have a semiconducting framework into which metal atoms can be substituted, providing a number of possibilities for new materials,Roy et al. (1992) for example new thermoelectrics,Nolas et al. (1998); Sales et al. (2001) and superconductors.Kawaji et al. (1995) Ferromagnetism has been observed in Mn-doped Ge clathrate, with $`T_c`$ = 10 K,Kawaguchi et al. (2000) and clathrates having Eu encapsulated in the cages, with $`T_c`$ up to 37 K.Sales et al. (2001) Currently there is strong interest in ferromagnetic semiconductors for spintronic applications,Ohno et al. (2000) and much effort has been devoted to the magnetic doping of semiconductors using low-temperature epitaxy.Ohno (2003) Clathrates offer the possibility of semiconducting phases containing magnetic atoms at equilibrium, which might be compatible with conventional substrates.Munetoh et al. (2001) Here, we report the synthesis of a stable Fe-doped Ge clathrate, and the observation of a ferromagnetic transition with $`T_c`$ = 170 K.

Although several transition metals substitute for Ge in the type-I clathrate structure,Cordier and Woll (1991) attempts at Fe substitution have failed due to competing Fe-Ge phases. However, we find that higher temperature synthesis produces a Fe-containing clathrate of the recently discovered chiral structure type, Ba6Ge25.Fukuoka et al. (2000); Carrillo-Cabrera et al. (2000); Kim et al. (2000) This clathrate contains an open Ge-Ge bonded framework, threaded by closed dodecahedral cages. The material actually forms as Ba6Ge24, with vacancies occurring preferentially by a Zintl mechanism.Kim et al. (2000)

Our method followed several steps. Elemental mixtures, with excess Ba, were placed in BN crucibles, sealed under Ar, and pre-reacted by rf induction. The resulting mixtures included type-I clathrates and FeGe2. Characterization was by x-ray diffraction (Bruker D8 Advance Powder, Cu $`K\alpha `$ radiation), and refinement using GSAS software.Larson and von Dreele (2000) To complete the reaction and eliminate FeGe2, the mixtures were arc-melted several times in a water-cooled copper crucible. The ingots were then sealed in quartz ampoules and annealed at 700 C for 100 hours. By this method we obtained single-phase or nearly single-phase Ba6Ge25-xFex with $`x`$ up to 3. However, furnace melting and slower cooling through the melting point produced a mixture of FeGe2 and clathrate phases. An x-ray analysis for $`x`$ = 3 is shown in Fig. 1. All observed peaks were indexed according to the cubic space group P4132 (#213), as previously reported for Ba6Ge25.Fukuoka et al. (2000); Carrillo-Cabrera et al. (2000); Kim et al. (2000) The lattice constant for $`x`$ = 3 is $`a`$ = 1.45520 nm, slightly larger than reported for Ba6Ge24 ($`a`$ = 1.45483 nm).Kim et al. (2000)

Energy-dispersive x-ray spectroscopy (EDS) and electron diffraction measurements (JEOL JEM-2010 electron microscope, 100 kV) for the $`x`$ = 3 sample showed that Fe was indeed incorporated in the clathrate phase. These measurements also showed trace amounts of Ge and FeGe2 to be present, but below the x-ray detectability limit (a few percent). Thus a minor fraction of the Fe atoms are present as FeGe2. No other phases were observed. Analysis gave an average clathrate composition Ba6Ge21.9Fe3.1 for the $`x`$ = 3 sample, a Ge:Fe ratio close to the starting composition. The Fe content per formula unit was 3.1 $`\pm `$ 1, with variations observed between crystallites, however the crystallites appeared locally homogeneous, with no apparent clustering.

The chiral clathrate structure (Fig. 1) has Ba situated in open channels as well as closed cages. In the Ge25 framework, two of the six distinct sites have three Ge-Ge bonds, the others having four. By contrast, type I and II clathrates have only 4-bonded Ge sites in approximately tetrahedral bonding configurations, and closed cages. In our x-ray analysis we placed Fe predominantly on three of the six sites, similar to the occupation parameters for In-doped Ge25 clathrate,Kim et al. (2000) and the refinement gave a composition Ba6Ge21.8Fe3.0, in excellent agreement with the EDS result (Ba6Ge21.9Fe3.1). However, since the Fe and Ge x-ray structure factors are similar, we found that Fe could be redistributed among framework sites to give a similar composition and nearly identical goodness of fit ($`R`$ values: $`R_{wp}`$ = 0.077 and $`R_p`$ = 0.053). The specific Fe locations therefore are uncertain, although it is clear that they are distributed in the clathrate framework.

Magnetic measurements were made on a Ba6Ge22Fe3 sample using a SQUID magnetometer (Quantum Design MPMS-XL). Zero-field-cooled (ZFC) magnetization was measured on warming in fixed field, while field-cooled (FC) magnetization was measured vs. decreasing temperature in constant field. Fig. 2 shows data for 0.1 T applied field. The high-temperature behavior was fit to a Curie-Weiss law, $`\chi =\frac{C}{TT_c}+\chi _{dia}`$, where $`\chi _{dia}`$ denotes a diamagnetic background. The value $`C`$ = 5.6 $`\times 10^5`$ K m3/kg indicates an effective moment peff = 5.5 $`\mu _B`$ per Fe, close to the free-ion value, 5.9 $`\mu _B`$, for high-spin Fe (peff = $`g[J(J+1)]^{1/2}`$). $`T_c`$ obtained from the fit is 180 K. The diamagnetic term, $`\chi _{dia}`$ = $`8.2\times 10^8`$ m3/kg, is large and may be an overestimation due to rounding of the transition. Some distribution of $`T_c`$ seems likely because of the observed Fe concentration distribution. However, low temperature magnetization measurements agree with the Curie fit; $`T_c`$ = 170 K is estimated from a modified Arrott plot (inset to Fig. 3).

A divergence of $`\chi _{FC}`$ and $`\chi _{ZFC}`$ starts near 170 K, and becomes much more pronounced below 110 K (Fig. 2). The remanence was also measured, warming the sample after removal of the field, and found to decrease with temperature, and vanish near 110 K. The magnetization at 2 K (Fig. 3) is polarized to 0.75 $`\mu _B`$ per Fe in a small field. However, full saturation is not reached even at 7 T, where the magnetization corresponds to 1 $`\mu _B`$ per Fe. This value is considerably lower than expected, given the value peff = 5.5 $`\mu _B`$ obtained from the Curie fit. This suggests a canted or noncollinear spin configuration.

The observed magnetic behavior cannot be due to the small amount of FeGe2, which has a spiral spin structure below 289 K, becoming commensurate at 263 K.Corliss et al. (1985) These do not correspond to the observed transitions, and the observed moments clearly could not result from inclusions of this phase at the few percent level.

To further understand the magnetic behavior, in-phase ($`\chi ^{}`$) and out-of-phase ($`\chi ^{\prime \prime }`$) ac susceptibilities were measured under a 0.6 mT ac field, with the results shown in Fig. 4. The higher-temperature data could be fit to a Curie law with $`T_c`$ = 180 K, in accord with the dc results. $`\chi ^{}`$ exhibits a strong rise below 200 K, corresponding to the ferromagnetic transition. However, below a maximum, $`\chi ^{}`$ decreases at lower temperatures and becomes frequency-dependent, while $`\chi ^{\prime \prime }`$ also becomes frequency-dependent. The frequency dependence, combined with the decrease of ZFC vs. FC magnetization below 110 K, are characteristic of reentrant spin glass systems.Binder and Young (1986); Mydosh (1993) We conclude that the transition near 170 K is ferromagnetic, while the 110 K irreversibility point corresponds to a spin glass freezing transition.

Frequency dependent $`\chi _{AC}`$ is found in spin glasses, while for ordinary ferromagnets the time-scale is much too small, even for the critical slowing down region near $`T_c`$.Binder and Young (1986); Mydosh (1993) Canted ferromagnets with a low-temperature antiferromagnetism have exhibited ac behavior similar to that of a reentrant spin glass,Mukherjee et al. (1994) however $`\chi ^{\prime \prime }`$ in that case shows only an increase with frequency, and no shift in canting transition temperature. The variation of $`\chi ^{\prime \prime }`$ peak position, $`T_p`$, per decade of frequency, $`[dT_p/d(\mathrm{log}f)/T_p]`$, is 0.042 for Ba6Ge22Fe3, which is in the range of values observed for spin glasses, 0.004 to 0.08.Mydosh (1993)

Spin-glass reentrance is attributed to the freezing of transverse degrees of freedom in a ferromagnet containing disorder.Gabay and Toulouse (1981) Alternatively, a cluster-glass mechanism has been argued,Rakers and Beck (1987) and it is clear that ferromagnetic clusters and considerable disorder in the ferromagnetic stateBao et al. (1999) do play a role in real materials. However, a cluster glass alone does not explain the present data, particularly the low-temperature saturation. A remanence increase on cooling in the glassy state is also typical of such spin glasses.Wynn et al. (1998)

Fe occupancy of framework sites is below the percolation threshold, even for second-neighbor bonds, so it is likely that the mechanism for magnetic ordering is conduction-electron mediation, as in diluted magnetic semiconductors. For comparison, in (Ga1-xMnx)As, Mn carries about 5 $`\mu _B`$, and $`T_c`$ scales linearly with small values of $`x`$, reaching $`T_c`$ = 110 K for $`x`$ = 0.05.Ohno (2003) The present system is quite similar in $`T_c`$ vs. $`x`$. Band-structure calculationsZerec et al. (2002) indicate a gap below the Fermi level in Ba6Ge25 which may quite possibly be the EF location in the Zintl composition Ba6Ge24. The Zintl conceptKauzlarich (1996) involves valence counting to give a closed-shell bonding configuration. Assuming valence 3 for Fe, Ba6Ge$`_{21\frac{3}{4}}`$Fe3 is also a Zintl phase, and very close to the average composition found here. Thus it seems likely that EF is near a band edge in this material, and the semiconductor model is appropriate.

There is considerable interest in the effects of disorder on the magnetic behavior of magnetic semiconductors;Schliemann and MacDonald (2002); Zaránd and Jankó (2002); das Sarma et al. (2003) the conduction-electron mediated interaction provides frustration, giving a non-collinear ordered state.Schliemann and MacDonald (2002) For the case of (Ga,Mn)N, the ground state appears to be a spin glass with no ferromagnetism.Dhar et al. (2003) Here we observe a ferromagnetic transition and spin-glass state in the same material.

In summary, a Ba6Ge22Fe3 clathrate was synthesized, exhibiting a ferromagnetic $`T_c`$ = 170 K. Below 110 K, a transition to a reentrant spin glass was evidenced by the onset of frequency-dependent ac susceptibility and a characteristic field-dependent magnetization.

###### Acknowledgements.

This work was supported by the Robert A. Welch Foundation, Grant No. A-1526, and by the National Science Foundation (DMR-0103455).


text:
# Model for Relaxation Oscillations of Luminous Accretion Disk in GRS 1915+105: Variable Inner Edge

## 1 Introduction

Microquasars in our Galaxy are enigmatic objects on account of their curious time variability, extremely large brightness (Belloni et al. 1997a, 1997b; Mirabel et al. 1994; Rao, Yadav, & Paul 2000), and energetic super-luminal jets (Mirabel & Rodriguez 1998). Large numbers of broad-band observations have been so far performed, revealing that they are low mass X-ray binaries including black holes and that the accretion flow around the black hole illuminates via the release of the gravitational energy. Their long-term bursts ($``$ months) could be triggered by a thermal instability in the disk (e.g. Mineshige & Wheeler 1989), however, the mechanism causing their short-term peculiar behavior is not well understood yet (Mirabel & Rodriguez 1999; Belloni et al. 2000; Naik et al. 2002). From the theoretical viewpoints, it is interesting to investigate how their time-dependent properties of accretion disks can be understood in the frame work of the disk models.

Among plenty of Galactic black-hole candidates, GRS 1915+105, most famous microquasar, is unique in exhibiting quasi-periodic luminosity variation (Belloni et al. 1997a, 1997b). Yamaoka, Ueda & Inoue (2001) analyzed the X-ray data of GRS 1915+105, taken by ASCA and RXTE during 1998-2000 and found that it exhibits mainly two branches in its high luminosity state: the high temperature branch (HTB), $`kT2`$ keV, and the low temperature branch (LTB), $`kT1`$ keV (see also Belloni 2002). The transition time between them is $`1`$ s. Interestingly, both states can be fitted by the multi-color disk model (hereafter, MCD; Mitsuda et al. 1984), however, what is puzzling in their HTB data is too high temperature and too small inner radius ($`R_{\mathrm{in}}20`$ km) to account for within the framework of the standard disk model for black-hole candidates (BHCs). It is interesting to note, in this respect, that a similar puzzle is known in ultra-luminous compact X-ray sources (ULXs, see Makishima et al. 2000, Mizuno, Kubota & Makishima 2001). In this case, several possibilities have been discussed; rotating black-hole hypothesis (Makishima et al. 2000), moderate beaming effects (King et al. 2001), and the shift of the inner edge of the disk in super-critical accretion regimes (Watarai, Mizuno, & Mineshige 2001, hereafter WMM01).

Then, why is GRS1915+105 so unique? What parameter is peculiar to GRS 1915+105? The fact that its luminosity is very large, comparable to the Eddington luminosity, indicates that this source seems to suffer very high mass accretion flow. The accretion disk theory predicts that an accretion disk becomes secularly and thermally unstable, when the mass accretion rate is very high so that the disk is radiation pressure-dominated (Lightman & Eardley 1974; Shibazaki & H$`\overline{\mathrm{o}}`$shi 1975; Shakura & Sunyaev 1976). What are then the observable consequences of such an unstable accretion disk?

As for the solution of accretion flow at accretion rates exceeding the critical rate ($`\dot{M}_{\mathrm{crit}}=L_\mathrm{E}/c^2`$), Abramowicz et al. (1988) found a self-consistent, steady-state solution, which is now known as the “slim disk.” It was then suggested that a near-critical disk may undergo limit-cycle oscillations, like the case of dwarf-nova outbursts, between the gas pressure-dominated standard disk and the slim disk (Kato 1983; Abramowicz, Lasota, & Xu 1986). Honma, Matsumoto & Kato (1991, hereafter HMK91) were the first to numerically demonstrate limit-cycle oscillations between the two different states based on the slim disk model. They show that an existence of the enough radiation pressure dominated region is essential for developing the thermal instability. Szuszkiewicz & Miller (1997, 1998) confirmed their results for wider $`\alpha `$ parameter ranges. The X-ray spectral variability of these oscillations was calculated by Zampieri, Turolla, & Szuszkiewicz (2001), who discussed the spectral behavior in each state and the corresponding hardening factor. More recently, Szuszkiewicz & Miller (2001) calculated the non-linear evolution of the disk with diffusion type viscosity, which is a more physical, and could also reproduce limit-cycle oscillation as in previous papers.

In the present study, we pay special attention to the behavior of the innermost region. The steady calculation of the slim disk indicates that a super-critical accretion flow can emit substantially from the part even inside $`3r_\mathrm{g}`$, since large amount of gas exists inside $`3r_\mathrm{g}`$ because of high mass accretion flow. Therefore, it has been anticipated to see a temporary reduction of the radius of the inner edge of the disk during a burst (Watarai et al. 2000; Mineshige et al. 2000). To examine if this is the case, we carefully calculate the time evolution of the innermost part of a super-critical flow, compare the theoretical results with observations, and discuss the possibility of obtaining a simpler proof for the existence of the super-critical accretion and thus probing the disk-instability to be a cause of the light variations of GRS 1915+105. Note that previous calculations were all performed down only to the radius of $`2.7r_\mathrm{g}`$ (as in HMK91) or $`2.5r_\mathrm{g}`$ (Szuszkiewicz & Miller 1997, 1998, 2001; Zampieri et al. 2001) so that an important fraction of the total radiation release might have been missed. In the present study, we calculate the disk structure down to $``$ 2.2 – 2.4$`r_\mathrm{g}`$, well inside the sonic radius, to see an observable feature in terms of the fitting parameters more clearly. This smaller boundary is critical to test the disk-instability scenario for GRS 1915+105.

In the next section we present the basic equations and our model assumptions. We show the results of the time-dependent calculations of the flow in section 3. The comparison with the X-ray observation of the microquasar GRS 1915+105 is discussed in section 4. The last section is devoted to conclusions.

## 2 Basic Equations and Numerical Methods

The methods of calculation are nearly the same as those by HMK91 (see Kato, Fukue, and Mineshige 1998 for the detailed discussion on the basic equations). The basic equations are written by using the cylindrical coordinates ($`r,\phi ,z`$). The assumptions of our calculation are; (i) a disk is axisymmetric, (ii) pseudo-Newtonian potential (Paczyńsky & Wiita 1980) is adopted; that is, $`\psi =GM/(Rr_\mathrm{g})`$ with $`R\sqrt{r^2+z^2}`$, and (iii) we employ one zone approximation toward the vertical direction. Here, $`r_\mathrm{g}`$ is the Schwarzschild radius defined by $`2GM/c^2=3\times 10^6(M/10M_{})`$ cm. We use the height integrated quantities, such as

$$\mathrm{\Sigma }=_H^H\rho 𝑑z=2I_\mathrm{N}\rho H,$$
(1)
and

$$\mathrm{\Pi }=_H^Hp𝑑z=2I_{\mathrm{N}+1}pH.$$
(2)
Here, $`\mathrm{\Sigma }`$, $`\mathrm{\Pi }`$, $`\rho `$, $`p`$, and $`H`$ are surface density, height integrated pressure, mass density, total pressure, and scale height, respectively. The coefficients, $`I_N`$ and $`I_{N+1}`$, are introduced by H$`\overline{\mathrm{o}}`$shi (1977). The density and pressure are related to each other by the polytropic relation, $`p\rho ^{1+1/N}`$. We assign $`N=3`$ in the entire calculations (i.e., $`I_3=16/35`$ and $`I_4=128/315`$). The equation of state is

$$p=p_{\mathrm{rad}}+p_{\mathrm{gas}}=\frac{a}{3}T_c^4+\frac{k_\mathrm{B}}{\overline{\mu }m_\mathrm{H}}\rho T_c,$$
(3)
where the first term on the right hand side represents the radiation pressure ($`a`$ and $`T_c`$ are the radiation constant and the temperature on the equatorial plane, respectively) and the second term represents the gas pressure ($`k_\mathrm{B}`$ is the Stefan-Boltzmann constant, $`\overline{\mu }=0.5`$ is the average molecular weight, and $`m_\mathrm{H}`$ is the hydrogen mass). We assume hydrostatic equilibrium (H$`\overline{\mathrm{o}}`$shi 1977),

$$(2N+3)\frac{\mathrm{\Pi }}{\mathrm{\Sigma }}=H^2\mathrm{\Omega }_\mathrm{K}^2H^2\frac{GM}{r(rr_\mathrm{g})^2}.$$
(4)
Here, $`\mathrm{\Omega }_\mathrm{K}`$ denotes the Keplerian angular frequency under the pseudo-Newtonian potential.

Mass conservation, the radial component of the momentum conservation, and angular momentum conservation are written as follows:

$$\frac{}{t}(r\mathrm{\Sigma })+\frac{}{r}(r\mathrm{\Sigma }v_r)=0,$$
(5)
$$\frac{}{t}(r\mathrm{\Sigma }v_r)+\frac{}{r}(r\mathrm{\Sigma }v_r^2)+r\frac{\mathrm{\Pi }}{r}+\mathrm{\Pi }\frac{\mathrm{ln}\mathrm{\Omega }_\mathrm{K}}{\mathrm{ln}r}=\mathrm{\Sigma }(l^2l_\mathrm{K}^2),$$
(6)
and

$$\frac{}{t}(r^2\mathrm{\Sigma }v_\phi )+\frac{}{r}(r^2\mathrm{\Sigma }v_rv_\phi )=\frac{}{r}(r^2T_{r\phi }),$$
(7)
Here, the velocity is expressed by $`v_r`$ and $`v_\phi `$ for the radial and azimuthal components, respectively, and the angular momentum of gas is expressed by $`l`$ and $`l_\mathrm{K}`$, which are defined by $`lrv_\phi `$ and $`l_\mathrm{K}r^2\mathrm{\Omega }_\mathrm{K}`$, respectively.

The $`r`$-$`\phi `$ component of viscous stress tensor in equation (7) is prescribed as

$$T_{r\phi }_H^Ht_{r\phi }𝑑z=\alpha \beta ^\mu \mathrm{\Pi },$$
(8)
where $`\mu `$ is a parameter ($`0\mu 1`$) and $`\beta p_{\mathrm{gas}}/p`$ (Szuszkiewicz 1990; HMK91; Watarai & Mineshige 2001). In the limit of vanishing $`\mu `$, we recover the ordinary $`\alpha `$ viscosity prescription (Shakura & Sunyaev 1973). If we take the limit of $`\mu 1`$, on the other hand, the viscosity depends on the gas pressure only (Sakimoto & Coroniti 1981).

Finally, the energy equation is

$$\frac{}{t}(r\mathrm{\Sigma }ϵ_{\mathrm{tot}})+\frac{}{r}[r(\mathrm{\Sigma }ϵ_{\mathrm{tot}}+\mathrm{\Pi })v_rrT_{r\phi }v_\phi ]=2rF^{},$$
(9)
in which advective cooling, viscous, and radiative cooling are considered. The explicit form of the total energy ($`ϵ_{\mathrm{tot}}`$) is

$$ϵ_{\mathrm{tot}}=\left[3(1\beta )+\frac{\beta }{\gamma 1}+\frac{1}{2}\right]\frac{\mathrm{\Pi }}{\mathrm{\Sigma }}+\frac{1}{2}(v_r^2+v_\phi ^2)+\psi _0(r),$$
(10)
where the first term on the right-hand side is the internal energy of gas ($`\gamma `$ is the adiabatic index and we set $`\gamma =5/3`$ in the present calculation), the second term represents the kinetic energy, and the last term is the potential energy on the equatorial plane, $`\psi _0(r)GM/(rr_\mathrm{g})`$. Radiative cooling flux per unit surface area of optically thick medium is given by

$$F^{}=\frac{8acT_c^4}{3\tau },$$
(11)
where, $`\tau `$ is the optical depth:

$$\tau =\overline{\kappa }\mathrm{\Sigma }=(\kappa _{\mathrm{es}}+\kappa _{\mathrm{ff}})\mathrm{\Sigma },$$
(12)
where $`\overline{\kappa }`$ is average opacity, $`\kappa _{\mathrm{es}}=0.4`$ is a opacity of the electron scattering, $`\kappa _{\mathrm{ff}}=0.64\times 10^{23}\overline{\rho }\overline{T}^{7/2}`$ is the absorption opacity via thermal Bremsstrahlung, and $`\overline{\rho }=16/35\rho `$ and $`\overline{T}=2/3T_c`$ are vertically averaged density and temperature, respectively. We solve equations (5)–(7) and (9) by modified Lax-Wendroff method with artificial viscosity.

To examine the behavior of the innermost part of the accretion flow, the calculations are performed from the outer radius at $`2000r_\mathrm{g}`$ down to the inner radius $`2.2r_\mathrm{g}`$ through the transonic point (Matsumoto et al. 1984). The number of total mesh is $`400500`$ and each mesh point is distributed uniformly on the logarithmic scale at large radii, $`r>3.2r_\mathrm{g}`$, and also uniformly but on the linear scale at smaller radii, where each mesh size is $`0.01r_\mathrm{g}`$. We fixed the black-hole mass to be $`m=10`$ ($`m=M/M_{}`$, where $`M_{}`$ is solar mass) and the viscosity parameter to be $`\alpha =0.1`$. The critical accretion rate, $`\dot{M}_{\mathrm{crit}}`$, is defined by $`L_\mathrm{E}/c^2`$, where $`L_\mathrm{E}`$ is the Eddington luminosity and $`c`$ is the speed of light. We define the dimensionless accretion rate to be $`\dot{m}\dot{M}/\dot{M}_{\mathrm{crit}}=\dot{M}c^2/L_\mathrm{E}`$ throughout the present study.

## 3 Time Evolution of the Inner Disk

### 3.1 Bolometric Light Curves

We first show the calculated light curves for models with $`\mu =0`$ (by the solid curve), 0.1 (by dotted curve) and 0.2 (by the dashed curve) in figure 1. The bolometric luminosity is calculated by integrating $`2\pi F^{}(r)`$ from the inner edge at $`r=2.2r_\mathrm{g}`$ (at $`r=2.4r_\mathrm{g}`$ for $`\mu =0.1`$ case) to the outer edge at $`r=2000r_\mathrm{g}`$. Eddington luminosity for a $`10M_{}`$ black hole is $`1.25\times 10^{39}`$ ($`M/10M_{}`$) erg $`\mathrm{s}^1`$. For simplicity, we assume blackbody radiation field (see section 3.4 for corrections).

In model with $`\mu =0`$, the disk luminosity suddenly rises on timescale of $`5(M/10M_{})`$s, keeps a high value for $`10(M/10M_{})`$s, and then decays. These time scales roughly correspond to the viscous timescale, $`t_{\mathrm{vis}}\alpha ^1\mathrm{\Omega }^1(r/H)^2`$, in the inner region. Namely,

$$t_{\mathrm{vis}}1.6\left(\frac{M}{10M_{}}\right)\left(\frac{r}{5r_\mathrm{g}}\right)^{3/2}\left(\frac{\alpha }{0.1}\right)^1\left(\frac{H/r}{0.1}\right)^2sec.$$
(13)
The amplitude of the burst is about one order of magnitude. Here, we calculated the same model as that in HMK91, i.e., $`\mu =0`$, $`\dot{m}=0.96`$, $`m=10`$, and $`\alpha =0.1`$. We confirmed that the amplitude and the duration of the burst in our calculation agree with those of HMK91 ($`t_{\mathrm{burst}}15`$ second for this parameter set). We continued calculations until $``$ 2000 second (that is, 2$`\times 10^7`$ ($`r_\mathrm{g}/c`$) calculation steps), however, we could not obtain a next burst in the same parameters as those in HMK91. In our calculation, small oscillations (or small bursts) took place near to the critical accretion rate for many times. Moreover, $`\mu =0`$ model takes a lot of recurrence time so that observational feature of GRS1915+105 does not agree with $`\mu =0`$. Hence, we will show $`\mu =0.1`$ model instead of $`\mu =0`$ model in the next section.

In model with $`\mu =0.2`$, in contrast, the amplitude of the burst is less, only about $``$ 2 mag (about a factor of $`6`$ variations). The burst rise time is somewhat longer than in the previous case but not significantly. Both of the bolometric light curves look quite similar to those obtained by HMK91, who also found that the burst amplitude and time scale are sensitive to the parameter, $`\mu `$. The different burst amplitudes can be understood in terms of the different shapes of the S-shaped thermal equilibrium curves, since the ranges of surface density (and thus of effective temperature), in which solutions are unstable, shrinks with increasing $`\mu `$ (HMK91). We will explain detailed evolution in the next two subsections.

### 3.2 State transition: case of $`\mu =0.1`$

As we mentioned in previous subsection, we can not obtain perfectly relaxed solutions in case of $`\mu =0`$ so that we will show $`\mu =0.1`$ case. We apply the data at third burst in figure 1, because the first burst is affected by initial conditions.

Figure 2 illustrates the evolution of the radial distributions of the surface density (top) and effective temperature (bottom). We set mass input rate of $`\dot{m}_{\mathrm{in}}=3.0`$ and $`\mu =0.1`$. The numbers in each panel indicate the order of the evolution sequence from third burst of $`\mu =0.1`$ model in figure 1.

The solid curve indicated by the number ‘1’ represents the state before the burst. In the prior burst state, radiation pressure exceeds the gas pressure from $`5r_\mathrm{g}`$ to $`20r_\mathrm{g}`$.

Initially, there is little mass inside $`3r_\mathrm{g}`$ (this radius corresponds to the radius of the marginally stable circular orbit around a non-rotating black hole). Surface density abruptly drops by 4 orders of magnitude inside $`3r_\mathrm{g}`$ in the initial state, because the radial velocity of the gas prominently increases near the transonic point. After $``$ 5$`(M/10M_{})`$s, a burst starts. A thermal instability is ignited in the innermost region and an upward transition from the radiation pressure-dominated standard-disk branch to the slim-disk one is initiated there. Then, a transition wave propagates outward and the outer zone subsequently undergoes an upward transition. Once the instability sets out, the efficiency of angular momentum extraction becomes largely enhanced, thus, gas in the disk being able to fall into the black-hole rapidly. Thus, surface density in the inner part shows a rapid decline. This enhanced mass flow fills the inner empty zone, and, hence, the surface density inside $`3r_\mathrm{g}`$ suddenly increases by more than one order of magnitude.

The effective temperature distribution also shows interesting features. The most important feature seen around the peak is that because of enhanced mass flow toward the innermost part the temperature inside $`3r_\mathrm{g}`$ keeps a rather high value, $`>10^7`$ K, so that significant amount of radiation is now expected from there. This property can be seen in figure 11 of Abramowicz et al. (1988) and in figure 2 of Watarai et al. (2000) have explicitly demonstrated that this is observable. In the high, slim-disk state, the maximum temperature is slightly higher and the inner disk radius becomes smaller than in the low, standard-disk state. The radial distribution of the effective temperature in time sequence numbers 3 are proportional to $`r^{1/2}`$, which is in good agreement with the analytical prediction (Watarai & Fukue 1999; Wang & Zhou 1999). We plotted $`r^{1/2},r^{3/4}`$ slopes in figure 2b, 3b for reference.

Another notable feature is found in the subsequent evolution. When the instability wave reaches $`100r_\mathrm{g}`$, a downward transition from the slim-disk branch to the standard-disk one begins from the inner region and again propagates outward. In both phases the transition is initiated from the inner region and the instability is propagated outward. Although these features are basically the same as those of the previous works (HMK91), we find interesting features in the behavior of the innermost region, which has been unclear in the previous calculations because of the larger truncation radius. This produces interesting observable effects. It should be noted that effective optical thickness at the inner region is more than unity during the burst phase. Hence, we can expect thermal emission around the disk inner region.

On the other hand, we have confirmed that a relatively large viscosity parameter, $`\alpha 0.1`$, is required to reproduce the observations of GRS1915+105. If we adopt a smaller viscosity parameter, which tends to increase surface density and thus optical depth, the duration becomes longer than the observed value. Therefore, we conclude that $`\alpha `$ is likely to be about 0.1 and some additional sources, such as metals, may contribute and raise opacities.

### 3.3 Different viscosity prescription: case of $`\mu =0.2`$

Honma et al. (1991) used the viscosity prescriptions given by equation (8) and discussed how different prescriptions affect the overall disk evolution. They have demonstrated that in the case with $`\mu =0.25`$ the burst duration is shorter and the amplitude is smaller both by a factor of $`2`$, compared with the case with $`\mu =0`$. That is, modifying viscosity prescription leads to drastic changes in the amplitudes of light variations and in the burst duration. We also test the case with $`\mu =0.2`$ and display the resultant light variations in figure 1 and the structural evolution in figure 3. Radiation pressure dominates over gas pressure from $`4r_\mathrm{g}`$ to $`65r_\mathrm{g}`$ in the initial state ($`\dot{m}_{\mathrm{in}}=3.0`$).

Interestingly, some differences are found. Although the surface density and the temperature of the innermost region are still variable, the way of time changes are qualitatively different from the previous case. The upper panel of figure 3 indicates that the instability evolves less efficiently than in the case with $`\mu =0`$. Especially, surface density remains very low inside $`3r_\mathrm{g}`$ even at the peak and the maximum temperature does not exceed $`10^7`$ K. Moreover, the maximum temperature is reached at $`46r_\mathrm{g}`$ unlike the previous case, in which temperature is maximum at the inner edge.

The reason for this change can be understood as follows: Non-zero $`\mu `$ values lead to a reduction of $`\alpha `$ when the radiation pressure is larger than the gas pressure ($`p_{\mathrm{rad}}>p_{\mathrm{gas}}`$). This reduction tends to increase the surface density of the equilibrium solution at a fixed mass accretion rate (since $`\alpha \mathrm{\Sigma }\dot{M}`$ in the steady condition), and the amount of this increase gets larger, as $`\dot{M}`$ increases. In other words, the S shape of the thermal equilibrium curves get less prominent, as $`\mu `$ increases (see Kato et al. 1998, figure 10.4). As a result, the thermal instability tends to be suppressed.

### 3.4 Model Spectral Fitting

We next perform the spectral fitting to the calculated spectra in the same manner as that adopted by WMM01. First, we calculate the spectra based on our calculation, assuming blackbody radiation at each radius for each time step. Then, we fit the spectra with the MCD model (Mitsuda et al. 1984), from which we can derive the disk maximum temperature, $`T_{\mathrm{in}}^0`$, and the radius, $`R_{\mathrm{in}}^0`$, of the region emitting blackbody radiation with $`T_{\mathrm{in}}^0`$. Before making a direct comparison with the observations we need the following corrections: (i) The derived value of $`R_{\mathrm{in}}^0`$ is influenced by the effects of the inner boundary condition as well as general relativistic effects; that is, the value tends to be overestimated, in general, compared with the real inner-edge radius. We thus need to introduce a correction factor $`\eta =0.41`$ and set $`R_{\mathrm{in}}=\eta R_{\mathrm{in}}^0`$ (Kubota et al. 1998). (ii) Second, the disk spectra are not perfectly blackbody but are modified by the Compton scattering effects, which lead to spectral hardening (Ross, Fabian, & Mineshige 1992; Shimura & Takahara 1995). Zampieri et al. (2001) also estimated the spectral hardening factor through solving the radiation transfer in the vertical direction, and found results, which mostly agree with Shimura & Takahara (1995). The actual value is $`\kappa 1.7`$ for stellar mass black hole and increases with increasing $`\dot{M}`$.

Unfortunately, it is difficult to precisely predict the spectral hardening factor due mainly to the uncertainties in metal opacities (which depend on the metal ionization stages, amount of incident X-ray photons, etc). Poorly understood disk-corona structure should also modify the model spectra. Also it strongly depends on the line of sight of observation, and the radial motion of the flow. Previous works (Shimura & Takahara 1995; Wang et al. 1999; Zampieri et al.2001) solved the radiation transfer in the vertical direction, assuming that, the photons can escape only in the vertical direction. For high accretion rate disk, however, the advected flow can trap significant fraction of photons; the photons propagate not only in the vertical direction but also in the radial direction. According to the preliminary results by Ohsuga, Mineshige, & Watarai (2003), the spectral hardening factor changes from 1.2 to about 2.3. This hardening factor depends on the disk model. Thus, the spectral hardening factor possibly exceeds 2.0. In recent study Shimura & Manmoto (2002) also show that the general relativistic effect is important for solving the radiation transfer equations including the photon trapping. According to their results the hardening factor is more than $`\kappa >100`$ in extreme Kerr hole and high mass accretion rate case ($`\dot{m}`$=100). We have not yet reached the consensus what the most appropriate value is. Therefore, we treated the hardening factor as a parameter.

Here we adopt relatively large value $`\kappa =2.0`$ because of high $`\dot{M}`$ cases being considered here and set the color temperature to be $`T_{\mathrm{in}}=\kappa T_{\mathrm{in}}^0`$. The effects of inclination angle and Doppler beaming are ignored in the present fitting.

Figure 4 displays the time variations of the X-ray luminosity, maximum temperature ($`T_{\mathrm{in}}`$), and the inner edge radius of the disk ($`R_{\mathrm{in}}`$), from the top to bottom, respectively (we also fitted the initial burst for $`\mu =0`$ model).

We calculated the X-ray luminosity by using $`R_{\mathrm{in}}`$ and $`T_{\mathrm{in}}`$, following equation (4) in Makishima et al. (2000). In the rise phase the luminosity increases by about one order of magnitude, $`T_{\mathrm{in}}`$ increases by a factor of $``$ 3, and $`R_{\mathrm{in}}`$ decreases by a factor of $`2`$. At the end of the decay phase (corresponding to the time sequence 4) we see a transient rise in $`R_{\mathrm{in}}`$. This is because in this phase the temperature profile is rather flat up to $`50r_\mathrm{g}`$ so that the fitting is affected by the emission from this extended regions. However, the transient rise may depend on the energy range used for fitting.

We note that to see the significant variation of $`R_{\mathrm{in}}`$, we should set the inner boundary of the disk as small as possible. Interestingly, however, evolutionary paths are different for the rise and decay phases in these diagrams. We will discuss about the comparison between our model and the observed objects in the next section.

## 4 Discussion

### 4.1 Comparison with the observation of Microquasar GRS 1915+105

We have calculated the time evolution of a near-critical accretion disk undergoing the relaxation oscillations between the low mass accretion rate branches and high mass accretion rate branches. We, especially, focus our study on the time changes of the apparent inner-edge radius, $`R_{\mathrm{in}}`$, and the maximum temperature, $`T_{\mathrm{in}}`$, finding peculiar behavior during the course of one burst cycle. In this subsection we discuss to what extent our model calculation can explain the X-ray observations of GRS 1915+105.

Figure 5 shows the results of the fitting to the observational data of GRS 1915+105 made by Yamaoka et al. (2001) based on the MCD model. It is clearly seen in this figure that the temperature increases and the inner radius decreases in the burst phase (see also Taam, Chen, & Swank 1997). This is what we obtained in the calculations with $`\mu =0.1`$, thus we can conclude that qualitative features can be explained by our model.

However, the precise fitting is not easy. According to model calculations, the burst properties can vary by changing the black-hole mass, $`\alpha `$ parameter, or viscosity prescriptions. Although the observed time changes of $`T_{\mathrm{in}}`$ and $`R_{\mathrm{in}}`$ agree well when $`\mu =0.1`$, the burst amplitude is too large to reconcile with the observation. Model with $`\mu =0.2`$, in contrast, can reproduce the observed amplitudes of the luminosity variations during the burst, but the shape of the light curve does not match the observation in the sense that the burst rise is slower than that of the observation. The observed time changes in $`T_{\mathrm{in}}`$ and $`R_{\mathrm{in}}`$ are not reproduced, neither. That is, we cannot perfectly fit the observations by simply changing the $`\mu `$ value. Likewise, changing $`\alpha `$ values does not alter the results qualitatively, but just results in changing the timescales. If we take smaller $`\alpha `$, say $`\alpha =0.01`$, we will obtain longer burst duration, about 19 s, since the burst duration time is $`\alpha ^{0.64}`$ (HMK91).

We should also remark that the absolute value of temperature is not well reproduced, since $`T_{\mathrm{in}}1`$keV in the observations, whereas $`T_{\mathrm{in}}0.6`$ keV in our model, even if we adopt a large spectral hardening factor of $`\kappa 2.0`$. To fit the observations during the low-luminosity state, we need to require that there should be additional mechanisms to further harden the spectral so as to achieve $`\kappa 3`$. As mentioned Introduction, this issue is common to other high-luminosity sources, such as ULXs.

There also exist other uncertainties arising due to jet or wind, corona, Compton scattering effects, and so on. Also, the determination of the inclination angle is very important. Recently some authors argue that dissipation of energy in disk corona or in jets can affect the disk evolution (Nayakshin et al. 2000; Janiuk, Czerny, & Siemignowska 2000). They propose the phenomenological model, in which they assume the total gas energy in both disk and jet components to be conserved and seek for the functional form, which matches the observed burst duration and luminosity. It is then found that that the fraction of the jet energy has to increase during the burst phase to match the observation. However, there is no strong physical basis on this suggestion.

If some dissipation processes like corona, jet or wind work around the peak, the amplitude of the burst get smaller, which might help. In other words, we need to somehow modify the S-shaped thermal equilibrium curves by changing $`\alpha `$ as a function of temperature, by including dissipation processes in jets or corona at the peak, or something else. It might be interesting to note that the disk-instability model for dwarf-nova outbursts has a similar problem; the constant $`\alpha `$ model cannot produce clear-cut outburst-quiescence light curves but instead only exhibits small-amplitude fluctuations (see, e.g., figure 1 of Mineshige & Osaki 1985). To fit the observations, therefore, one needs to increase $`\alpha `$ values in the hot state, compared with that in the cool state, by a factor of $`3`$.

### 4.2 Is there a Kerr black hole?

In the high mass-accretion-rate flow, the disk inner edge can be smaller than $`3r_\mathrm{g}`$ even for a non-rotating black-hole (Watarai et al. 2000). On the other hand, it is well known that the case of a Kerr black-hole can also explain the high disk temperature and small inner-disk radius, as are observed, because of smaller radius of the marginally stable orbit. If the black hole is an extreme Kerr black hole, especially, the disk inner edge should be as small as $`0.5r_\mathrm{g}`$.

Recently, mass estimation was made for the black hole of the microquasar GRS 1915+105, which is $`14\pm 4M_{}`$ (Greiner et al. 2001). The corresponding Schwarzschild radius is 42$`\pm 12`$ km. If we take the face value derived by Yamaoka et al. (2001), the inner radius of the disk should be $`R_{\mathrm{in}}R_{\mathrm{in}}\sqrt{\mathrm{cos}i}30`$ km during the burst, while $`R_{\mathrm{in}}21`$ km for the extreme Kerr hole. This result seems to rule out the extreme Kerr hole interpretation at least for GRS 1915+105, however, we should keep in mind that we are only able to detect photons practically from the region outside $`2r_\mathrm{g}`$ ($``$ 82 km) owing to the significant photon redshifts. It thus follows that it is not easy to distinguish whether there exists an extreme Kerr black hole or not.

At present, therefore, both possibilities (the Kerr hole and the super-critical accretion) are viable, but we can at least rule out the extreme Kerr hole case, since then $`R_{\mathrm{in}}`$ should stay constant at a small value, not depending on the luminosity of the disk, contrary to the observations.

### 4.3 Uniqueness of GRS 1915+105

The greatest mystery involved with GRS 1915+105 resides in why only this source exhibits various types of activities, including quasi-periodic light oscillations. High mass accretion rates may not be the only reason, since other luminous sources, like GS 2023+338, whose luminosities seem to reach the Eddington luminosity (Tanaka 1989) did not record similar behavior. Further, Kubota et al. (2001) analyzed the X-ray spectrum in some Galactic black hole candidates, revealing that XTE J1550–564, one of the most well-observed microquasars with RXTE, does not behave like GRS 1915+105 during the outburst, although the spectrum evolution is somewhat similar to that of GRS 1915+105; i.e., $`T_{\mathrm{in}}`$ increases and $`R_{\mathrm{in}}`$ decreases, as the luminosity increases in their burst phase. They also indicate that the spectrum of XTE J1550–564 can be best fit if we add the anomalous (Comptonized) component in the decay phase, whereas no such component is needed in the rise phase. This means, the disk structure differs in the rise and decay phases. What discriminates these objects and other canonical black hole X-ray transients? Do other bright ultra-luminous X-ray sources show bursting properties similar to GRS 1915+105 but not observed yet? These still remain unsolved problems.

It is interesting to note in this respect that for relatively large viscous parameter, $`\alpha 0.3`$, a transition from the optically thick standard-disk branch to the optically thin ADAF occurs (Takeuchi & Mineshige 1998). The transition timescale is less than a few sec for a stellar mass black-hole ($`M_{\mathrm{BH}}=10M_{}`$). It might be that the solution is not unique and some hidden parameter triggers the different behavior of near-critical accretion flow.

## 5 Conclusions

We calculated the evolution of the inner part of the high-luminosity accretion disk, carefully considering the transonic nature of the flow in the innermost region. We obtain the following results:

1. Our calculations confirm the previous results that the state transitions between the standard-disk and slim-disk branches occur during the burst phase. Around the burst peak, the inner edge of the disk lies at $`<3r_\mathrm{g}`$ and there are significant photons coming from the region inside $`3r_\mathrm{g}`$.

2. To inspect the X-ray observational features in GRS1915+105, we performed spectral fitting to the calculated spectra with the fitting parameters of $`T_{\mathrm{in}}`$ and $`R_{\mathrm{in}}`$. As to model with $`\mu =0.1`$ we find an abrupt transient increase in $`T_{\mathrm{in}}`$ and a temporary decrease in $`R_{\mathrm{in}}`$ during an outburst, both of which were actually observed in microquasar GRS 1915+105. However, the amplitudes of light variations and the absolute value of $`T_{\mathrm{in}}`$ are too large to reconcile with the observational data. When $`\mu =0.2`$, in contrast, the burst amplitude agrees with the observed one, but the behavior of $`T_{\mathrm{in}}`$ and $`R_{\mathrm{in}}`$ differ from the observed one. It is difficult to match the observation by simply varying $`\alpha `$ or $`\mu `$. It is also important to note that the fitting results are sensitive to the photon energy range used for the fitting.

3. To summarize, although the precise fitting is left as future work, we can safely conclude that our “limit-cycle oscillation” model can explain the basic observational features of GRS 1915+105. It is then required that the spectral hardening factor at high luminosities close to the Eddington luminosity should be about 3 instead of less than 2 as is usually assumed. This means, we need additional mechanism that causes a further spectral hardening.

The authors thank the Yukawa Institute for Theoretical Physics at Kyoto University. Discussions during the YITP workshop YITP-W-01-17 on “Black Holes, Gravitational Lens, and Gamma-Ray Bursts” were useful to complete this work. Especially, we would like to thank R. Narayan, C. Done, K. Yamaoka and T. Kawaguchi for useful comments and discussions. Numerical computation in this work was carried out at the Yukawa Institute Computer Facility. This work was supported in part by the Grants-in Aid of the Ministry of Education, Science, Sports, and Culture of Japan (14001680, KW).


text:
# FAMOS: A Dynamically Configurable System for Fast Simulation and Reconstruction for CMS

## I Introduction

The detailed simulation and reconstruction of physics events at the Large Hadron Collider (LHC), $`p\overline{p}`$ collisions at 14 TeV, where a large number of high-energy particles are produced and up to 200 collisions overlay, is extremely time consuming even on todays fastest CPUs. However, to investigate the physics potential of the CMS detector a large number of collisions must be studied. This is especially true for scanning the multi-dimensional parameter space of SUSY models. By replacing the detailed modeling of energy losses and material effects in simulation as well as the reconstruction algorithms by simple parameterizations significant improvements concerning speed can be at the price of precision losses.

## II The Simulation, Reconstruction and Analysis chain

The procedure to analyze physics models can be decomposed into four major steps:

1. Generation of Monte-Carlo events.
The Monte-Carlo events are created using generators like Pythia, Herwig or ISAJET. These generators produce a list of particles - stable and decaying - and their four-vectors. In CMS the information is stored in form of HBOOK-Ntuples. The time required to generate an event is normally less than 100 milliseconds222All times are given as measured on a 1 GHz Pentium III..
2. Simulation of material effects.
This is the most time consuming step required. Currently, the simulation is still done with a Geant3 g3  based FORTRAN program: CMSIM cmsim . The time varies significantly depending on the number of particles and their energies in the Monte-Carlo events. The average is between 100 and 200 seconds per event. The information that is stored as output of this step are called SimHits. They contain the information about the energy stored in different detector elements at different times. A new, object-oriented simulation program OSCAR oscar , based on Geant4 g4 , is currently being validated.
3. Simulation of readout electronics (digitization).
The detector converts the energy deposited by the particles into electronic signals that are converted to digital information by ADCs and TDCs. At the high luminosity at the LHC the detector “sees” the overlap of up to 200 minimum bias events with a single signal event. Since the simulation of material effects requires large amount of CPU time, the minimum bias events are randomly selected from a large pool of simulated events and combined with the simulated signal events. So even being technically part of the simulation, the combination of minimum bias events with a signal event and the simulation of the detector response to the energy deposition are performed by the reconstruction software: ORCA orca . The time for this depends on the simulated luminosity and the event type. It is 1 to 10 seconds per collision. The output created in this step is called DIGIs.
4. Reconstruction of physics/analysis objects.
The reconstruction is performed by ORCA in several sub-steps. First the DIGIs are combined to reconstructed hits, RecHits, which for example combine several strips of the silicon tracking detectors. Similarly, level-1 trigger objects are built. Then RecHits are used to find tracks in the inner tracker and the muon chambers and clusters in the calorimeters. The reconstruction can produce more complicated objects like jets or information about the missing energy and finally physics objects like electrons, photons, muons etc. The time spent on this can vary largely but is typically 10 to 100 seconds per collision.

The total CPU time required before the analysis of a collision can be started is 3-5 minutes.

## III Structure of FAMOS

To give physicists the possibility to study large event sample fast, the FAMOS333FAst MOnte-carlo Simulation, see famos  project has been developed in CMS. The main design concept is to provide fast modules, able to replace each single step in the chain, as well as several steps in one go and this complying to the same interfaces as the full simulation and reconstruction.

To achieve this FAMOS is divided into several subsystems that focus on the different steps. Figure 1 shows the different subsystems and their relations to the other CMS object-oriented software projects. COBRA/CARF carf  is the general framework that controls event handling, DDD ddd  provides services for geometry handling, OSCAR is the simulation software and ORCA the reconstruction software. The subsystems for a particular task can be accessed via several interfaces: One providing the same interfaces as for the corresponding detailed simulation or reconstruction, one for simplified and easy standalone usage. The implementation of OSCARInterface and ORCAInterface allows to dynamically change from detailed simulation and reconstruction to the FAMOS implementations.

The first implementation focuses on the replacement of all steps for simulation and reconstruction and produces physics objects directly from the Monte-Carlo four-vectors. This corresponds to the hashed area in Figure 1. The rest of this article describes the current implementation.

## IV Event handling and simulation modules

When implementing the direct way from Monte-Carlo events to physics analysis objects, two main issues arise. First the Monte-Carlo event is read and second dedicated simulation modules for the different sub-detectors of CMS are executed.

### IV.1 Monte-Carlo event reading

An important first step in FAMOS is the handling of the original Monte-Carlo event. The event is stored in the class RawHepEvent that is an exact C++ implementation of the HEPEVT FORTRAN CommonBlock that is used by most Monte-Carlo generators. An abstract base class BaseHepEventReader exists that allows to fill the RawHepEvent from a multitude of sources. ASCII files, HBOOK Ntuples, particle guns - mostly used for single particle tests, from the databases used as persistency store by OSCAR and ORCA and directly from the Pythia-6 Monte-Carlo generator. Similar reader modules can be created for other Monte-Carlo generators.

The event reader modules are provided by the CMS framework, COBRA, and used from there.

### IV.2 FAMOS event handling

In FAMOS an special event manager class, FamosEventMgr, takes care of

* reading the Monte-Carlo event and
* calling simulation modules that inherit from FamosSimulator

It interfaces directly to the reconstruction on demand mechanism that is a key point of the CMS framework. Figure 2 shows a collaboration diagram between the FamosEventMgr and the framework classes. At startup FamosEventMgr accepts registrations from the selected FamosSimulator modules and reads the geometry. Concrete classes that inherit from FamosSimulator implement the fast simulation algorithms.

When starting to loop over events, the FamosEventMgr first fills RawHepEvent from the selected input. RawHepEvent is also used to keep track of new particles that can be created by the FamosSimulator classes. This modification of the event is done to take care of effects like Bremsstrahlung, pair production or multiple scattering.

The FamosSimulator classes are called in a well defined sequence. It is important to follow the order in which the particles cross the detector - first the Tracker, then the Calorimeters (ECAL and HCAL) and only last the Muon system.

Currently, the simulators are called in the order of their registration - the user is responsible to do this correctly.

### IV.3 Simulation modules and interfaces to the framework

The classes inheriting from FamosSimulator are implementing the concrete fast simulation parameterizations. Currently, implementations exist for

* material effects
* tracking (Tracker and Muon)
* electromagnetic clustering
* muon level-1 trigger

The simulation modules - residing in the GeneratorToAnalysis subsystem have to implement a method
bool reconstruct(RawHepEvent &);
and can add methods to provide the objects they simulate or reconstruct. All modules provide access to objects that are kept as simple as possible and a special, framework specific mechanism is then applied to provide the results in a form compatible with the interfaces used in normal ORCA jobs: RecObj

classes. This is done in the ORCAInterface subsystem. Direct use of the simulators can be achieved using the StandaloneInterface subsystem.

Figure 3 shows how classes from COBRA (LazyObserver, RecUnit), ORCA (TTrack) and FAMOS (generic: FamosSimulator, GeneratorToAnalysis: BaseFastTracksFromGenerator, ORCAInterface: RecTrackFromParticle, FastTrackReconstructor) work together to provide Tracker tracks produced by the fast simulation modules with identical interface to the ORCA analysis program as tracks from the full reconstruction software.

## V Example configuration and results

Since FAMOS provides ORCA compliant interfaces no change is required to the part of code that uses the objects (e.g. TTracks) when changing from full to fast reconstruction. However, instantiation of the FAMOS objects and their registration to the framework must be done. For Tracker tracks to be provided by the FATSIM simulation module this looks like the following:

```
myTrackFinder = new RecTrackFromRawParticle(
                  new FATSIM( &myPtEG, &myAngleEG,
                              &myEff,
                              myImpactParameterEG));
myTrackReconstructor = new FastTrackReconstructor(
                          myTrackFinder, "FATSIM");
[...]
RecCollection<TTrack> MyTracks(ev->recEvent(),
                               "FATSIM");

```

The first lines instantiates an ORCA track-finding object with the FATSIM simulation module as argument. This simulation module uses internal classes that provide the parameterizations of transverse momentum, angular distributions and the efficiencies. Then the ORCA track reconstructor is created an a string - “FATSIM” - is used to identify it. Creating the collections of TTrack specifying the same string will issue the reconstruction on demand, i.e. when iterated over the collection.

This can be compared to the registration of the regular track finding algorithm:

```
myTrackFinder = new TrackReconstructor(
                 new CombinatorialTrackFinder,
                 "FkFTracks");
RecCollection<TTrack> MyTracks(ev->recEvent(),
                               "FkFTracks");

```

Again a track-finder is instantiated and then a collection. For the collection the ONLY difference is the string to specify which reconstruction algorithm is to be used when the collection is accessed.

Similarly, when linking the executable different sets of libraries must be selected when using full or fast reconstruction. The main use-case is to run the fast simulation and reconstruction directly with a Monte-Carlo generator as input. However, it is also possible to use the FAMOS in parallel to a regular full reconstruction from a prepared database. This is especially useful for direct comparison between full and fast simulation and reconstruction. The original Monte-Carlo event that was simulated in a time-consuming process is stored in the database, read from there by the FamosEventMgr and then handed to the fast simulators. That allows for example to analyze two RecCollections of TTracks in the same program since the different collections are mapped via the string identifier to different track-finding algorithms: detailed reconstruction with RecHits as input and FAMOS with the Monte-Carlo event as input.

## VI Timing

The main aim of the current first implementation was to achieve the highest possible speed with reasonable agreement. The measured results for high-multiplicity events without minimum bias pileup is shown in Table 1.

The FAMOS time is dominated by the Monte-Carlo generator (Pythia) and the Calorimeter simulation module. All other modules - in particular the FAMOS framework (FamosEventMgr) contribute only insignificantly. The performance achieved is more than 3000 times faster than the detailed simulation and reconstruction.

It is possible to simulate the CMS detector and in particular the material effects more precisely. A first prototype has been developed that takes pair production and Bremsstrahlung into account. Without optimization for performance this reduces the speed of FAMOS to be only factor 100 faster than CMSIM+ORCA. The modularity - this is just one additional simulation module to register to FamosEventMgr \- give the user full flexibility to adapt the simulation to the precision his concrete physics analysis requires.

## VII Summary

FAMOS is a high-performing, flexible and dynamically configurable mini-framework for fast simulation and reconstruction. It is fully integrated in the general CMS framework but the components can be used independently. It is possible to mix full and fast simulation with minimal changes to the user code. The first modules implemented provide high-level reconstructed objects (Tracks, Muon, Muon Lvl-1 trigger, Calorimeter clusters) at reasonable agreement in precision about 3000 times faster than full simulation and reconstruction.

In future fast modules for individual steps in the analysis chain will be developed and geometry reading will be using the same XML description as used by the full simulation to ensure consistency.

###### Acknowledgements.

Special thanks to Marco Battaglia, Filip Moortgat and Artur Kalinowsky for their implementations and the plots showing example results.


text:
# Reconstruction of electrons with the Gaussian-sum filter in the CMS tracker at LHC

## I Introduction

Modern track detectors based on semiconductor technologies contain larger amounts of material than gaseous detector types, partially due to the detector elements themselves and partially due to additional material required for on-sensor electronics, power, cooling, and mechanical support. A precise modelling of material effects in track reconstruction is therefore necessary to obtain the best estimates of the track parameters. Such material effects are particularly relevant for the reconstruction of electrons which, in addition to ionization energy loss and multiple Coulomb scattering, suffer from large energy losses due to bremsstrahlung.

A well-known model of the bremsstrahlung energy loss is due to Bethe and Heitler BetheHeitler34 . In this model, the probability density function (PDF), $`f(z)`$, of the energy loss of an electron is

$$f(z)=\frac{\left[\mathrm{ln}z\right]^{c1}}{\mathrm{\Gamma }(c)},$$
(1)
where $`c=t/\mathrm{ln}2`$, $`t`$ is the thickness of material traversed by the electron (in units of radiation length), and $`z`$ is the fraction of energy remaining after the material layer is traversed. The probability of a given fractional energy loss is assumed to be independent of the energy of the incoming particle. This PDF is shown in Fig. 1 for different thickness values.

The baseline for track reconstruction in the CMS tracker is the Kalman filter Rudi87 . Throughout the filter tracks are described by a five-dimensional state vector, containing the information about the momentum, the direction and the position at some reference surface. The material effects are currently assumed to be concentrated in the active elements of the detector layers. In this context the optimal treatment of radiative energy loss is to correct the momentum with the mean value of energy loss and to increase the variance of the momentum by adding the variance of the energy loss distribution. This procedure should ensure unbiased estimates of the track parameters and of the associated uncertainties Stampfer94 . The Kalman filter is a linear least-squares estimator, and is proved to be optimal only when all probability densities encountered during the track reconstruction procedure are Gaussian. The implicit assumption of approximating the Bethe-Heitler distribution with a single Gaussian is quite crude. It is therefore plausible that a non-linear estimator which takes the actual shape of the distribution into account can do better.

A non-linear generalization of the Kalman filter (KF), the Gaussian-sum filter (GSF) Rudi97 ; Rudi98 , has therefore been implemented in the reconstruction software of the CMS tracker ORCA . In the GSF the distributions of all state vectors are Gaussian mixtures, i.e. weighted sums of Gaussians instead of single Gaussians. The algorithm is therefore appropriate if the probability densities involved in track reconstruction can be adequately described by Gaussian mixtures. The basic idea of the present work is to approximate the Bethe-Heitler distribution as a Gaussian mixture rather than a single Gaussian, in which the different components of the mixture model different degrees of hardness of the bremsstrahlung in the layer under consideration. The resulting estimator resembles a set of Kalman filters running in parallel, where each Kalman filter corresponds to one of the components of the mixture describing the distribution of the state vector.

## II Approximating the fractional energy loss distribution

An important issue with the GSF reconstruction of electrons is to obtain a good Gaussian-mixture approximation of the Bethe-Heitler distribution. The parameters to be obtained are the weights, the mean values and the variances of each of the components in the approximating mixture. The parameters are determined by minimizing the following two distances:

$`D_{\text{CDF}}`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\left|F(z)G(z)\right|𝑑z,`$ (2)
$`D_{\text{KL}}`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\mathrm{ln}\left[f(z)/g(z)\right]f(z)𝑑z,`$ (3)
where $`f(z)`$ and $`F(z)`$ are the PDF and cumulative distribution function (CDF) of the model distribution and $`g(z)`$ and $`G(z)`$ are the PDF and CDF of the Gaussian mixture, respectively. The distance $`D_{\text{KL}}`$ is the so-called Kullback-Leibler distance between the model distribution and the mixture. Hereafter, the mixtures obtained by minimizing $`D_{\text{CDF}}`$ are called CDF-mixtures, whereas the mixtures obtained by minimizing $`D_{\text{KL}}`$ are called KL-mixtures. The minimizations have been done independently on a set of discrete values of $`t`$, ranging from 0.02 to 0.20. Figures 2 and 3 show the resulting distances as a function of thickness for a varying number of components in the approximating mixture.

In order to obtain mixtures for arbitrary values of the thickness, fifth-degree polynomials have been fitted to the parameters as a function of $`t`$. Due to the fast access to the parameters from the polynomials, the calculation of the mixture is done on the fly during reconstruction, using the effective thickness of a detector layer from the knowledge of the incident angle of inclination.

## III Reducing the number of components

The approximation of energy loss by a Gaussian mixture amounts to a convolution of this mixture with the current state, which in general is also composed of several Gaussian components. The strict application of the GSF algorithm therefore quickly leads to a prohibitively large number of components due to the combinatorics involved each time a layer of material is traversed.

In a realistic implementation of the GSF the number of components must repeatedly be reduced to a predefined maximum. As little information as possible should be lost in this procedure. Two strategies have been tested:

1. Only the N components with the largest weights are kept;
2. Components are merged into clusters, according to a given metric.

The first option has the advantage of being computationally light, but it turns out to be inferior. Even the first two moments of the estimated parameters are not described correctly.

In the second approach, the component with the largest weight is merged with the one closest to it, and this procedure is repeated until the required number of components is reached. The results below have been obtained by using the Kullback-Leibler distance – defined in Equation (3) – as a measure of distance.

## IV Results from simulated tracks in the CMS tracker

First, results from the reconstruction of data originating from a simplified simulation are shown. In this simulation multiple scattering and ionization energy loss are turned off, all the material is concentrated on the detector units, and the exact amount of material used in the simulation is known by the reconstruction program. Single electron tracks with $`p_T=10`$ GeV/$`c`$ have been simulated for absolute values of $`\eta `$ less than 1.0 . Reconstructed hits have been collected using the knowledge of the associated simulated hits, so no pattern recognition has been involved. The following results all refer to the quantity $`q/p`$ (charge over absolute value of the momentum) recorded at the point of closest approach to the vertex in the transverse plane – the transverse impact point (TIP) – after a fit going from the outside towards the inside of the tracker. Figure 4 shows an example of the estimated $`q/p`$ for one single track, both for the KF and for the GSF.

Figures 5 and 6 show probability distributions for the estimated $`q/p`$ of the KF and the GSF with a varying maximum number of components kept during the reconstruction. Given the estimated PDF (a single Gaussian for the KF, a Gaussian mixture for the GSF), each entry in the histogram amounts to the integral from $`\mathrm{}`$ to the true value of $`q/p`$. If the estimated PDF is a correct description of the real distribution of the parameter, the corresponding histogram should be flat.

The deviation from flatness can be quantified by the $`\chi ^2`$ of the difference between the probability distributions of $`q/p`$ and the flat distribution. This $`\chi ^2`$ per bin is shown in Fig. 7 for a set of different mixtures as a function of the maximum number of components kept. The CDF-mixtures are superior to the KL-mixtures concerning the quality of the estimated $`q/p`$. The main trend seems to be related to the maximum number of components kept rather than the number of components in the mixture describing the energy loss, even though the mixtures with five and six components are best in the limit of keeping a large number of components.

Figure 8 shows the residuals of the estimated $`q/p`$ of the GSF and the KF with respect to the true value of the parameter. The estimated $`q/p`$ for the GSF is the mean value of the state vector mixture, and the mixture used for this specific plot is a CDF-mixture with six components. In order to quantify the difference between the GSF and the KF residuals, the full-width at half-maximum (FWHM) and the half-width of intervals covering 50% and 90% of the distribution have been considered. The covering intervals have been chosen to be symmetric about zero. The FWHM and the half-widths of the covering intervals are shown in Figs. 910 and 11. The different flavours of the GSF in these figures are the same as those described in the caption of Fig. 7.

The GSF and the KF have also been run on tracks from a full simulation using the official CMS simulation program CMSIM . The $`p_T`$ and the $`\eta `$ range are the same as in the simplified simulation, but the amount and spatial distribution of the material are different. Probability distributions of the estimated $`q/p`$ for the GSF and the KF are shown in Fig. 12.

The probability distribution of the GSF exhibits no large deviation from flatness, indicating that the estimated PDF of $`q/p`$ describes reasonably well the actual PDF of $`q/p`$. This observation is all the more remarkable since, with the full simulation, the energy loss is not generated by the simple Bethe-Heitler model, and neither the exact amount nor the exact location of the material are known to the GSF.

The corresponding residuals of the estimated $`q/p`$ with respect to the true value are shown in Figs. 13 and 14. The residuals shown in Fig. 14 have been obtained by including a vertex constraint in the fit. Such a constraint allows the momentum to be measured in the innermost part of the track and thus gives a handle on possible radiation in the first two layers. The result of including this constraint is a less skew distribution with the mode being moved closer towards zero, and the amount of tracks in the tails is also reduced. Even though the results from the full simulation qualitatively seem to confirm those from the simplified simulation, more studies are needed to understand the differences in detail.

## V Conclusion

The Gaussian-sum filter has been implemented in the CMS reconstruction program. It has been validated with electron tracks with a simplified simulation in which the energy loss distribution (Bethe-Heitler model), the exact amount of material and its exact location are known to the reconstruction program. It has been shown that the quality of the momentum estimate depends mainly on the number of mixture components kept during reconstruction, and to some extent also on the number of components in the mixture approximation to the energy loss distribution. A comparison with the best linear unbiased estimator, the Kalman filter, shows a clear improvement of the momentum resolution. Remarkably, a similar improvement can be seen with electron tracks from the full simulation, although in this case neither the exact energy loss distribution nor the precise amount and location of material are known to the reconstruction program. More systematic studies with electrons from the full simulation are clearly needed, but it seems safe to conclude that in electron reconstruction the Gaussian-sum filter yields a substantial gain in precision as compared to the Kalman filter.