text:
# Role of finite layer thickness in spin-polarization of GaAs 2D electrons in strong parallel magnetic fields

## Abstract

We report measurements and calculations of the spin-polarization, induced by a parallel magnetic field, of interacting, dilute, two-dimensional electron systems confined to GaAs/AlGaAs heterostructures. The results reveal the crucial role the non-zero electron layer thickness plays: it causes a deformation of the energy surface in the presence of a parallel field, leading to enhanced values for the effective mass and g-factor and a non-linear spin-polarization with field.

The spin-polarization of an interacting, dilute two-dimensional (2D) carrier system has been of interest for decades. It has long been expected that because of Coulomb interaction the product $`g^{}m^{}`$, which determines the spin susceptibility of the 2D system, increases as the 2D density ($`n`$) is lowered and eventually diverges as the system makes a transition to a ferromagnetic state at sufficiently low $`n`$ tanatar ; attaccalite  ($`g^{}`$ and $`m^{}`$ are the carrier Landé g-factor and effective mass, respectively). Recently, there has been much renewed interest in this problem, thanks to the availability of high-quality dilute 2D systems, and the belief that it may shed light on the controversial issue of a metal-insulator transition in 2D abrahams . A technique commonly used to study the spin-polarization is to measure the response of the 2D system to a tilted or parallel magnetic field okamoto ; tutuc ; vitkalov ; pudalov ; tutuc02 ; proskuryakov ; noh ; etienne . The results of some of these measurements tutuc02 ; proskuryakov ; noh , however, appear to be at odds with what is theoretically expected tanatar ; attaccalite  for a dilute, interacting 2D system that is otherwise ideal, i.e., has zero layer thickness and is disorder-free. In particular, when $`g^{}m^{}`$ is deduced from parallel magnetic field at which the 2D system becomes fully spin-polarized, then the experimental results for GaAs 2D electrons tutuc02  and holes proskuryakov ; noh  suggest a decreasing value of $`g^{}m^{}`$ with decreasing $`n`$, opposite to the theoretical predictions.

Here we report a combination of measurements and calculations for the parallel magnetic field-induced spin-polarization of 2D electrons at the GaAs/AlGaAs heterojunction. The results highlight the importance of the finite thickness of the electron layer and the resulting deformation of the energy surface $`E(𝐤_{})`$, where $`𝐤_{}`$ is the in-plane wave vector, that occurs in the presence of a strong parallel field. This deformation induces an enhancement of both $`m^{}`$ and $`g^{}`$ and leads to a non-linear spin-polarization in a parallel field. We find that, once the effect of the finite layer thickness and interaction is taken into account, there is reasonable agreement between the experimental data and calculations.

We used five samples from three different wafers (A, B, and C). The samples were all modulation-doped GaAs/AlGaAs heterojunctions with $`n`$ in the range $`0.8`$ to $`6.5\times 10^{10}`$ cm-2. Their low-temperature mobility varied depending on the sample and $`n`$; at $`n=2\times 10`$ cm-2, it ranged from about $`2\times 10^5`$ to $`2\times 10^6`$ cm2/Vs. Samples were patterned in either van der Pauw or Hall bar shapes, and were fitted with back- or front-gates. To tune $`n`$ in samples A, B1, and B2, following illumination with a red LED, we used front-gate bias; for B3 (in the range $`n<4\times 10^{10}`$ cm-2) and C we used back-gate bias and no illumination. For B3, the highest density ($`n=6.5\times 10^{10}`$ cm-2) was obtained after illumination, followed by back-gating to reduce $`n`$ to $`4.5\times 10^{10}`$ cm-2. Measurements were done down to a temperature of 30 mK, and a rotating platform was used to tune the angle between the applied magnetic field and the sample plane.

Figure 1 summarizes our data taken on the different samples. Plotted are the values of $`g^{}m^{}/g_bm_b`$, determined from the parallel magnetic field, $`B_P`$, at which the 2D system becomes fully spin polarized ($`m_b=0.067m_0`$ and $`g_b=0.44`$ are the band effective mass and Landé g-factor for GaAs electrons; $`m_0`$ is the free electron mass). The parallel magnetic field ($`B_{}`$) leads to the formation of two energy subbands, one for each spin, and separated by the Zeeman energy, $`E_Z=|g^{}|\mu _BB_{}`$, where $`\mu _B`$ is the Bohr magneton. The 2D system becomes fully spin-polarized above a field $`B_P`$ at which $`E_Z`$ equals the Fermi energy. The equality leads to an expression for $`B_P`$: $`B_P=(h^2/2\pi \mu _B)(n/|g^{}|m^{})`$, from which we determine $`g^{}m^{}`$ that are plotted in Fig. 1.

The procedures we have used to experimentally determine $`B_P`$ have been described elsewhere tutuc ; tutuc02 ; here we give a brief summary. We determine $`B_P`$ from two independent sets of experiments: Shubnikov - de Haas (SdH) measurements in a nearly parallel magnetic field, and magnetoresistance measurements in a strictly parallel field. In the first type of experiment, we apply a constant magnetic field ($`B_{tot}`$) whose initial direction is parallel to the 2D electron plane, and then slowly rotate the sample while recording the sample resistance as a function of the angle between the 2D plane and the field direction. If we limit ourselves to small angles, the field’s parallel component ($`B_{}`$) remains almost constant (equal to $`B_{tot}`$ typically to better than 1$`\%`$) during the rotation, while its perpendicular component ($`B_{}`$) changes sufficiently to probe the SdH oscillations. We then Fourier analyze the SdH oscillations to obtain the populations of the two spin subbands. These densities provide a direct measure of the spin-polarization of 2D electron system and allow us to determine the field $`B_P`$ above which the system becomes fully spin-polarized. In the second type of experiment we measure the sample resistance as a function of a magnetic field applied strictly in the 2D plane. As shown elsewhere tutuc ; tutuc02 , the in-plane magnetoresistance shows a marked change in its functional form at the field $`B_P`$. For a given sample, the value of $`B_P`$ obtained from the two types of experiments (open and closed symbols in Fig. 1) are in agreement.

Data of Fig. 1 illustrate that the product $`g^{}m^{}`$, deduced from $`B_P`$ as described in the last two paragraphs, deviates substantially from what is expected for an ideal, interacting 2D electron system, namely a monotonically increasing $`g^{}m^{}`$ as $`n`$ is lowered tanatar ; attaccalite . The data also reveal that the measured $`g^{}m^{}`$ is sample dependent and not a unique function of $`n`$. A possible reason for this non-uniqueness may be the sample disorder that indeed varies between different samples. An examination of the data, however, argues against this hypothesis: considering the data at a given density, it is clear that there is no simple trend linking the sample disorder, as deduced from the low-temperature mobility, to the measured $`g^{}m^{}`$. As we demonstrate below, another factor that renders the experimental 2D electrons non-ideal, namely their finite layer thickness, appears to be responsible for the sample-dependent $`g^{}m^{}`$ and the difference between the observed and expected density dependence of $`g^{}m^{}`$.

In the presence of a large $`B_{}`$, when the magnetic length ($`=\sqrt{\mathrm{}/eB_{}}`$) becomes comparable to or smaller than the thickness of the electron layer, the energy surface $`E(𝐤_{})`$ of the electrons gets deformed in the in-plane direction perpendicular to $`B_{}`$. The deformation leads to an increase of the in-plane effective mass, $`m^{}`$, which, in second order perturbation theory, is given by stern

$$m^{}(B_{})=m_b/\sqrt{1\frac{2e^2B_{}^2}{m_b}\underset{j0}{}\frac{|z_{0j}|}{E_jE_0}}$$
(1)
where the sum runs over all excited subbands $`j`$ and $`z`$ is the quantization axis das sarma . Data of Fig. 2 provide an experimental demonstration of this effect in our samples. In Fig. 2(a) we show, as a function of $`B_{}`$, the measured magnetoresistance of samples A and B1, both at a density of $`2.7\times 10^{10}`$ cm-2. The magnetoresistance for each sample shows a clear change in its dependence on $`B_{}`$ at a field marked by a vertical arrow as $`B_P`$. As demonstrated previously tutuc ; tutuc02 , the field $`B_P`$ marks the onset of full spin-polarization. Note in Fig. 2(a) that $`B_P`$ is larger for sample B1 than A even though they have the same density. We also measured $`m^{}`$ in the two samples as a function of $`B_{}`$, as shown in Fig. 2(b) footnote1 . We determined $`m^{}`$ from the temperature dependence of the amplitude of the SdH oscillations, measured as the sample was slowly rotated in an almost parallel field. We performed a standard analysis, fitting the amplitude of the SdH oscillations ($`\mathrm{\Delta }R`$) to the Dingle formula, $`\mathrm{\Delta }R\xi /\mathrm{sinh}\xi `$, where $`\xi 2\pi ^2k_BT/\mathrm{}\omega _c`$ and $`\omega _c=eB/m^{}`$. It is clear in Fig. 2(b) that $`m^{}`$ for both samples exhibits a strong enhancement with increasing $`B_{}`$, consistent with Eq. (1). Moreover, $`m^{}`$ has a larger enhancement for sample A than for B1. This stronger enhancement correlates with the smaller $`B_P`$ measured for sample A. We believe that the main difference between the two samples in Figs. 2(a) and (b) is that sample A has a larger layer thickness: because of a decrease in subband separation, the larger thickness leads to the larger enhancement of $`m^{}`$, consistent with the smaller $`B_P`$ for sample A compared to B1.

To further substantiate the connection between layer thickness and $`m^{}`$ enhancement, in Fig. 2(c) we show data for sample A at a higher carrier density of $`4.2\times 10^{10}`$ cm-2. The measured $`m^{}`$ enhancement is smaller for the higher $`n`$ state. This is consistent with layer thickness being responsible for the $`m^{}`$ enhancement: as we use a more positive front-gate bias to increase $`n`$ in this sample, the electron wavefunction is squeezed more towards the interface so that the layer thickness is reduced \[see the insets to Fig. 2(c)\], consistent with the smaller measured $`m^{}`$ enhancement.

In order to quantitatively understand the experimental data, we have done self-consistent density-functional calculations of the subband structure in the presence of $`B_{}`$. We used the recent parameterization of the exchange-correlation energy by Attaccalite et al. attaccalite . The energy and length scales for electrons in a semiconductor are characterized by the effective Rydberg and the Bohr radius according to the effective mass $`m^{}`$ and the dielectric constant of the material. In our calculations it was crucial that $`m^{}`$ was determined as a function of $`B_{}`$ from the self-consistently calculated subband dispersion $`E(𝐤_{},B_{})`$. These calculations confirm the qualitative trends expected from Eq. (1). Band structure effects beyond the effective-mass approximation are not consisdered here. We have checked that these effects are of minor importance winkler93 . In the calculations, the field $`B_P`$ is defined as the smallest value of $`B_{}`$ for which the fully spin-polarized configuration has the lowest total energy.

Figure 3 provides an example of the results of the calculations for $`n=4.2\times 10^{10}`$ cm-2; shown are (a) the spin-subband densities $`n_{}`$ and $`n_{}`$, (b) $`m^{}`$, and (c) $`g^{}`$, as a function of $`B_{}`$. The calculations were done using the parameters of sample A (spacer thickness and barrier height), assuming a p-type background doping of $`2.7\times 10^{13}`$ cm-3, and a binding energy of 90 meV for the dopants (Si) in the barrier. These are reasonable values, consistent with our estimate of the unintentional (residual) doping in our molecular beam epitaxy system and the binding energies quoted in the literature schubert . The calculations predict a non-linear but smooth increase of the spin-polarization as a function of $`B_{}`$.

In Figs. 3 (a) and (b) we have also included our measured spin-subband densities and $`m^{}`$, determined from SdH oscillations in a nearly parallel magnetic field. Overall, there is good agreement between the experimental data and calculations for both spin-polarization and $`m^{}`$ as a function of $`B_{}`$. The calculated $`B_P=14.5`$ T agrees well with the measured $`B_P16`$ T, and is much smaller than $`B_P48`$ T expected for an interacting GaAs 2D electron system with zero layer thickness at $`n=4.2\times 10^{10}`$ cm-2 (see dashed curve in Fig. 4).

To understand the density dependence of $`B_P`$ for a given sample, we also calculated $`B_P`$ as a function of $`n`$, which is tuned either with a front- or back-gate bias. For these calculations, we kept the sample parameters fixed, and only changed the boundary conditions for the Hartree potential thus simulating the effect of the gate bias. The results for the case where the front-gate is used are shown in Fig. 4 (solid curve), along with the experimental data for samples A and B1. There is reasonable agreement between the calculations and the experimental data footnote2 .

In Fig. 4 we also plot $`B_P`$ vs. $`n`$ (dashed curve) calculated based on Ref. attaccalite  where an interacting 2D electron system with zero layer thickness is assumed, as well as $`\mathrm{"}B_P\mathrm{"}`$ (dotted curve) determined from the calculated spin susceptibility in the limit $`B_{}=0`$ attaccalite  and assuming a linear spin-polarization as a function $`B_{}`$. Three noteworthy trends are observed in Fig. 4. First, for $`n>0.5\times 10^{10}`$ cm-2 the dashed and dotted curves are very close to each other, meaning that the spin polarization of an ideal, zero-thickness 2D system is approximately (within 5%) linear with $`B_{}`$ in this density range or, equivalently, the Zeeman splitting can be expressed well in terms of an effective g-factor independent of $`B_{}`$. This implies that the finite layer thickness is the key factor that leads to the observed non-linearity in spin-polarization with $`B_{}`$ and the resulting reduction of $`B_P`$. The mechanism responsible for this non-linearity can be summarized as follows: $`B_{}`$ induces an increase of $`m^{}`$ due to the finite layer thickness. This has a twofold effect. It directly reduces $`B_P`$ because $`B_P1/m^{}`$. Furthermore, the increase of $`m^{}`$ reduces the effective Bohr radius and thus increases $`r_s`$, the average electron spacing measured in units of effective Bohr radius. The increase in $`r_s`$ in turn yields an increase of $`g^{}`$ due to the Coulomb interaction. Second, the solid and dashed curves merge as the density is lowered, consistent with the expectation that, because of the smaller $`B_P`$, the finite layer thickness induced $`g^{}m^{}`$ enhancement becomes less important. Third, at ultra-low densities the dashed and dotted curves start to diverge (see the lower inset), implying that the interaction alone can induce a significant non-linearity of spin-polarization with $`B_{}`$ in a very dilute 2D system attaccalite .

Finally, if we use the expression $`B_P=(h^2/2\pi \mu _B)(n/|g^{}|m^{})`$ to determine $`g^{}m^{}`$ from $`B_P`$ as frequently done in the literature, we obtain the curves shown in Fig. 4 upper inset. These plots emphasize that ”$`g^{}m^{}`$” deduced from $`B_P`$ for a 2D system with finite layer thickness (solid curve) is significantly enhanced with respect to the ideal 2D system and can show a non-monotonic dependence on $`n`$. The results therefore caution against extracting values for $`g^{}m^{}`$ in the limit of zero magnetic field from measurements of $`B_P`$ at large parallel fields zhu . Moreover, $`B_P`$ and $`g^{}m^{}`$ are not unique functions of $`n`$; they depend on the electron layer thickness which in turn depends on sample parameters and experimental conditions.

We thank NSF and DOE for support, D.M. Ceperley and S. Moroni for helpful discussions. Part of this work was done at NHMFL; we thank T. Murphy and E. Palm for support.


format: text

text:
# Distribution of infected mass in disease spreading in scale-free networks

## Abstract

We use scale-free networks to study properties of the infected mass $`M`$ of the network during a spreading process as a function of the infection probability $`q`$ and the structural scaling exponent $`\gamma `$. We use the standard SIR model and investigate in detail the distribution of $`M`$, We find that for dense networks this function is bimodal, while for sparse networks it is a smoothly decreasing function, with the distinction between the two being a function of $`q`$. We thus recover the full crossover transition from one case to the other. This has a result that on the same network a disease may die out immediately or persist for a considerable time, depending on the initial point where it was originated. Thus, we show that the disease evolution is significantly influenced by the structure of the underlying population.

There has been a growing interest recently in the network structure  and dynamics  of real-life organized systems. Many such systems, covering an extremely wide range of applications, have been recently shown  to exhibit scale-free character in their connectivity distribution, meaning that they obey a simple power law. Thus, the distribution of the connectivity of nodes, follows a law of the form

$$P(k)k^\gamma ,$$
(1)
where $`k`$ is the number of connections that a node has and $`\gamma `$ is a network parameter which determines the degree of its connectivity. These networks have some unusual properties, thus justifying the heavy interest in recent years. For example, spreading processes in scale-free networks show a dynamic behavior which is different than in other classes of networks. These processes are plausible models for the spreading of diseases, epidemics, etc. Several models for spreading exist in the literature built on different algorithms, such as the SIR model , the SIS model , the SIRS model , etc.

In a recent paper Newman  studied analytically the SIR model in scale-free networks. The connectivity distribution had an exponential cutoff of the form

$$P(k)k^\gamma e^{k/\kappa },$$
(2)
where $`\kappa `$ is an arbitrary cutoff value for $`k`$. This work gave a closed form solution for the epidemic size and the average outbreak size as a function of the infection probability. It showed that there is a critical infection threshold ($`q_c`$) only for small finite values of $`\kappa `$, but as $`\kappa `$ increases $`q_c`$ decreases, apparently resulting in $`q_c0`$ (no critical threshold at all) in the limit of $`\kappa \mathrm{}`$. The same result was also shown for the SIS model by Pastor-Satorras and Vespignani .

The absence of a critical threshold is not a universal network property. Actually, in the well studied cases of lattice networks and small-world networks the opposite is true . Such a threshold  is always present, which separates the infected from the uninfected regions. This has as a result that these models do not offer a very realistic picture. Recently, Warren et al.  have used a heterogeneous distribution for the infection probabilities both for lattices and small-world networks. They model the variability in a population which results in a broadening of the transition regime; however, a threshold still exists and the behavior of the transition is qualitatively similar to the case of the simpler SIR model on a lattice. Because of this, scale-free networks are distinctly different regarding the predictions on the rate and efficiency of spreading. This is clearly much closer to what it is intuitively expected, and can provide useful estimates for the properties of epidemics of any kind. The same type of model could also describe a diverse set of networks, such as social networks, virus spreading on the Web, rumor spreading, signal transmission etc.

In the present study we calculate the distribution of the epidemic size, i.e. the distribution of the infected mass, for several different $`\gamma `$ values. This property helps us to better understand the importance of the starting point (origin) of the disease. It turns out that this distribution is not a smooth function for all networks, but depends strongly on the network density, i.e. the value of $`\gamma `$.

We use a simulation algorithm to construct a scale-free network comprised of $`N`$ nodes. We follow a network generation method which enables us to freely vary the connectivity distribution of the network. We assign a number of edges $`k`$ for each node by using a power-law distribution $`P(k)k^\gamma `$. Starting from the highest connectivity nodes we create links by randomly choosing $`k`$ other nodes. Care is taken that no duplicate links are established between the same two nodes and once a node has reached the number of edges initially assigned to it, it no longer accepts any new connections. The cutoff value for the maximum possible connectivity of a node was fixed to $`N/2`$.

The spreading of a disease follows the standard SIR (Susceptible, Infective, Recovered) model. Initially, all nodes are in the susceptible (S) status, and a random node is infected (I). During the first time step it tries to infect with probability $`q`$ the nodes linked to it, and when the attempt is successful the status of the linked node switches from S to I. The process is repeated with all infected nodes trying to influence their susceptible (S) neighbors during each time step. After trying to infect its neighbors the status of an infected site changes to recovered (R) and can no longer be infected. The simulation stops when there are no infected nodes in the system or when all nodes have been infected. Small-world networks were constructed as described in reference .

We monitor the percentage $`M`$ of nodes infected and the duration of a disease (i.e. the time needed for the disease to either disappear or cover the entire network). In figure 1 we show for comparison purposes the percentage $`M`$ of infected sites as a function of $`q`$, for three different network types. The form of the curves for the lattice and the small world networks is similar, i.e. in all cases there is a sharp transition and a critical point. However, this behavior is unrealistic, as it does not follow the majority of actual situations, such as disease spreading in real-life networks . If this were to happen, e.g. the Web would be in a state of either no virus present or the entire Web (all computers in the world) would be infected, with a very small probability of having an intermediate situation with only a certain fraction of computers infected, which is the realistic picture. Scale-free networks with a high degree of connectivity ($`\gamma =2.0`$) follow a much smoother spreading evolution, as the mass of the infected population increases almost linearly with the infection probability q and there is no transition regime. This is in agreement with the recent formalism of Newman . The linear behavior can be understood as follows: On a scale-free network there exist nodes with a wide range of connectivity. For fixed infection probability $`q`$, the average probability for an infected node with $`k`$ links to spread the disease is $`kq`$. If this number is greater than 1, it is statistically certain that a neighbor node will be infected. If it is significantly less than 1 the disease will die out. On the other hand, for a small-world network and for a lattice network there is a characteristic mean number of links $`k`$ assigned to each node.

Here, for a scale-free network, we look at the distribution of the quantity $`M`$. This is shown in figure 2, where we see that for different $`\gamma `$ values we have two opposite situations. First, for $`\gamma =2.0`$ we see that for fixed $`q`$ we can have both cases coexisting, depending on the connectivity of the initially infected node. For this case the distribution $`\mathrm{\Phi }(M)`$ of the infected population for fixed $`q`$ is bimodal and comprises of two distinct parts, a strong peak at $`M=1/N`$ (only one node is infected) attributed to initially infected nodes with $`k`$=1, or more accurately $`kq1`$, and a Gaussian-shaped part at some higher value of $`M`$. The almost linear increase of $`M`$ with $`q`$ in figure 1 is due to the fact that for higher $`q`$ the peak of $`1/N`$ decreases ($`kq`$ increases), and the average value of the Gaussian distribution increases accordingly. We also notice in figure 2 that the gaussian part of the distribution has an average value which is larger than $`q`$. For example, for $`q=0.1`$ the peak of the distribution is close to 0.25, and when $`q=0.8`$ it is closer to 0.9. This means that if the peak in $`1/N`$ did not exist, the behavior of the curve in figure 1 would be superlinear, i.e. $`M`$ would always be larger than $`q`$ as a result of the complex connectivity of the network, but it would not reach the value of $`M=1`$ for infection probabilities less than 1. The interplay between the peak at $`1/N`$ and the gaussian part yields the final curve which follows roughly a linear increase.

This result is a superposition of the two possible states present in the case of a simple lattice, where below the threshold the distribution is a simple peak at $`M=1/N`$, while above the threshold the distribution peaks around $`M=1`$. On a scale-free network there exists a finite probability for the disease to be either eliminated immediately or cover a considerable portion of the network. Thus, for small $`\gamma `$ values, the prediction on the future of a disease largely depends on the place where it originates, because of the existence of the gap. This bimodality implies that every disease which survives the initial step(s) spreads over a non-zero portion of the network population. This is not the case for large $`\gamma `$ values, where recently , results for the SIR model on a scale-free network (with $`\gamma =3.0`$) were presented. The distribution of $`M`$ did not show any bimodality, similarly to our solid curve on the left of figure 2, which corresponds to $`\gamma =2.9`$ and $`q=0.5`$. We see that the bimodality is now lost, and it is replaced by a smoothly decreasing function with increasing $`M`$. This result is in agreement with the work of Moreno et al. . Thus, the distribution changes from bimodal to unimodal as we go from a dense network to a sparse network. The crossover from one to the other takes place as it is given in the inset in figure 2, where we plot $`\gamma `$ vs $`q`$, i.e. each ($`\gamma `$, $`q`$) pair is exactly at the corresponding transition point.

These observations show that in scale-free networks nodes of high connectivity act as “boosters” to the disease spreading; even if very few nodes remain infected, by the time a high connectivity node is infected it spreads the disease over a significant number of its neighbors, even for small $`q`$. This fact stresses the importance of the ‘hubs’, as it has also been observed in the past in studies of the static properties for such networks . Similarly, the low-connectivity nodes may serve to isolate large clusters of the network. These clusters are effectively screened by the disease via the presence of the low-connectivity node, especially in the case of high $`\gamma `$, where the network is loosely connected.

Similar conclusions can be drawn for the duration of a disease. For small-world networks and regular lattices the duration of a disease when epidemics takes place is practically constant for infection probabilities greater than the threshold (figure 3). For scale-free networks, on the contrary, there is a slow increase of the duration as a function of $`q`$. However, the duration is much smaller now, which implies that even for considerable infection probabilities a disease cannot last for a long time, and a considerable portion of the network can be infected in a practically small and constant time. The duration of the disease at $`q=1`$ is also a measure of the network ‘diameter’, since it represents the average number of links needed to cross, before reaching all the system nodes.

Upon monitoring the distribution of uninfected sites as a function of time we observed that it followed a power law at all times. The exponent of this power law was always the same as the one used for the initial connectivity distribution, with the curve scaled down by a constant factor. Thus, sites of different connectivity are infected with the same relative rate.

Summarizing, we have investigated spreading properties on scale-free networks. For the SIR model we studied we find that the starting point of the disease is very important, because it can either stop the disease or facilitate its spreading. This phenomenon for dense networks yields a bimodal distribution for the infected mass, with a peak close to 0 and a gaussian part around a finite value of $`M`$. Despite the smooth increase of the infected mass with $`q`$ the disease spreads rapidly on the network in a practically constant time, almost independently of $`q`$. This rapid spreading manifests the compactness of the network (as compared to lattice and small-world networks) and its small diameter which is related to the short path length from any site of the network to another. For sparse networks this is not the case, as it also has been previously observed.


text:
# Large Transport Properties of Bi2Sr2-xLaxCaCu2O8+δ Single Crystals Grown by a Floating-Zone Method

## I Introduction

It is well known that the high $`T_c`$ superconductivity occurs with doping carriers into the antiferromagnetic (AF) insulator. The AF insulator is quite different from a conventional insulator, and shows various peculiar features. One of the most peculiar features is that in the lightly doped La2-xSrxCuO4 (LSCO) metallic behavior of the in-plane resistivity $`\rho _{ab}`$ ($`d\rho _{ab}/dT>`$ 0) near room temperature was observed, even though the magnitude of the in-plane resistivity is far above the Mott limit ando . This metallic behavior is considered to be responsible for the appearance of band dispersion within the Mott-Hubbard gap seen in the angle-resolved photoemission spectroscopy (ARPES) measurement teppei . The electric states of the doped parent insulator, on the other hand, seem to be quite different between LSCO and the Bi-based cuprates. The charge-stripe order is stabilized in LSCO, whereas it is much weaker in the Bi-based cuprates. LSCO has been most frequently used for systematic studies of the heavily underdoped region of the cuprates, because a heavily underdoped sample is easy to synthesize. By contrast, there are only a few reports on the lightly doped Bi2Sr2CaCu2O8+δ terra ; kend , because it is difficult to grow high-quality single crystals of lightly doped region. We have previously grown the parent insulator of Bi2Sr2CaCu2O8+δ by flux method, where the trivalent lanthanide, such as Dy and Er, is substituted for divalent Ca terra . Here, we have successfully grown the parent insulator of Bi2Sr2-xLaxCaCu2O8+δ by a traveling solvent floating zone (TSFZ) method.

## II Experimental

Two crystals of 40% La doped Bi2Sr2CaCu2O8+δ were used in this study. One is a crystal in which all La atoms were substituted for Sr site (Bi2(Sr,La)2CaCu2O8+δ), and the other is a crystal in which La atoms were substituted for Sr and Ca sites (Bi2(Sr,Ca,La)3Cu2O8+δ). They were grown by TSFZ method with a growth velocity of 0.5mm/h in a mixed gas flow of O2 (20%) and Ar (80%). The starting composition of the feed rods were Bi2.1Sr1.1La0.8Ca1.0Cu2.0O8+δ for Bi2(Sr,La)2CaCu2O8+δ and Bi2.1Sr1.5La0.8Ca1.0Cu2.0O8+δ for Bi2(Sr,Ca,La)3Cu2O8+δ. Very large single crystals with a size of 2 $`\times `$ 4 $`\times `$ 0.1 mm3 were successfully grown. By energy dispersive x-ray analysis (EDX), the actual compositions were determined to be Bi2.01Sr1.11La0.84Ca1.00Cu2.00O8+δ for Bi2(Sr,La)2CaCu2O8+δ, and Bi2.04Sr1.31La0.86Ca0.76Cu2.00O8+δ for Bi2(Sr,Ca,La)3Cu2O8+δ. In the case of Bi2(Sr,Ca,La)3Cu2O8+δ, excess La replaces the Ca site. Figure 1 shows the XRD pattern for the as-grown Bi2(Sr,La)2CaCu2O8+δ single crystal. All the peaks were indexed as (0 0 $`l`$) peaks of the Bi-2212 phase with no trace of impurities. The $`c`$-axis length was evaluated to be 30.18 Å for Bi2(Sr,La)2CaCu2O8+δ, and 30.38 Å for Bi2(Sr,Ca,La)3Cu2O8+δ, which is smaller than that of Bi2Sr2CaCu2O8+δ (30.86 Å). This is due to the larger ion radius of the Sr atom (1.13 Å) than that of the La (1.05 Å) and Ca (1.00 Å) atoms. The in-plane resistivity of Bi2(Sr,La)2CaCu2O8+δ was measured using standard four-probe method, whereas the in-plane resistivity of the highly insulating sample Bi2(Sr,Ca,La)3Cu2O8+δ was measured using a two-probe method with constant voltage (where the sample resistance is much higher than the contact resistance). The voltage for both samples was confirmed to vary linearly with the current in the measured range from 4.2 to 300 K. The thermopower was measured using a steady-state technique from 4.2 to 300 K. The Hall coefficient measurement was done by sweeping the magnetic field which was applied along the c axis (c // H) from -7 to 7 T at fixed temperatures. The Hall resistance was confirmed to be linear in magnetic field. The crystals were annealed in the O2 or Ar atmosphere at 600 C for 10 h to change the oxygen contents as is similar to Bi2Sr2CaCu2O8+δ wata .

## III Results and discussion

Figure 2 shows the temperature dependence of the in-plane resistivities for Bi2(Sr,La)2CaCu2O8+δ (closed circles) and Bi2(Sr,Ca,La)3Cu2O8+δ (open circles). Bi2(Sr,Ca,La)3Cu2O8+δ shows insulating behavior ($`d\rho _{ab}/dT`$ $`<`$ 0) in the measured temperature range. The magnitude of the in-plane resistivity at room temperature is about 30 m$`\mathrm{\Omega }`$cm for an as-grown sample. It increases to 120 m$`\mathrm{\Omega }`$cm by annealing in Ar, which suggests that the carrier concentration decreases with decreasing the oxygen content. Bi2(Sr,La)2CaCu2O8+δ shows metallic behavior ($`d\rho _{ab}/dT`$ $`>`$ 0) above 150 K, and the magnitude of the in-plane resistivity is about 5 m$`\mathrm{\Omega }`$cm at room temperature. In the inset of Fig. 2, the in-plane resistivity of Bi2(Sr,La)2CaCu2O8+δ is plotted in linear scale to see the metallic behavior clearly. Although observed metallic behavior seems to suggest the Boltzmann transport, the magnitude of the in-plane resistivity is far above the Mott limit for metallic conduction, which is quite different from a conventional insulator (usually called as a “bad metal” emery ).

Figure 3 shows the temperature dependence of the in-plane thermopower for Bi2(Sr,La)2CaCu2O8+δ (closed circles) and Bi2(Sr,Ca,La)3Cu2O8+δ (open circles). For both samples, the thermopower decreases with decreasing temperature, which is consistent with our previous results of Bi2Sr2RCu2O8+δ (R=Dy, Er) terra2 . The magnitude of the thermopower increases by Ar annealing and decreases by O2 annealing, which means that the thermopower decreases with increasing oxygen contents. Since the room temperature thermopower is considered to be a universal measure of the doping level obert , this result indicates that the carrier doping is caused by the excess oxygen. Note that the room-temperature thermopower and the magnitude of the in-plane resistivity of Bi2(Sr,La)2CaCu2O8+δ are smaller than those of Bi2(Sr,Ca,La)3Cu2O8+δ. Although these two compounds contain approximately the same amount of La ($`x`$ $``$ 0.8), the carrier concentration is different. For example, the carrier concentration per Cu $`p`$ are estimated to be 0.02 for as-grown Bi2(Sr,Ca,La)3Cu2O8+δ, and 0.04 for as-grown Bi2(Sr,La)2CaCu2O8+δ. One explanation for this is that the substitution of La for Ca site causes a disorder into the CuO2 plane to make the doped carriers localized. Another explanation is that the amount of the excess oxygen is different.

Figure 4 shows the temperature dependence of the Hall coefficient RH for as-grown Bi2(Sr,La)2CaCu2O8+δ (closed circles) and as-grown Bi2(Sr,Ca,La)3Cu2O8+δ (open circles). The Hall coefficients monotonically increase with decreasing temperature for both samples. The room temperature carrier concentration is estimated from the relation $`n=(eR_H)^1`$ to be 4 $`\times `$ 1020 cm-3 ($`p`$ 0.04) for Bi2(Sr,La)2CaCu2O8+δ, and 1 $`\times `$ 1020 cm-3 ($`p`$ 0.01) for Bi2(Sr,Ca,La)3Cu2O8+δ, which agree well with the carrier concentration estimated from the room temperature thermopower. The Hall mobility calculated using the Hall coefficient and the resistivity as $`\mu R_H/\rho _{ab}`$ is shown in the inset of Fig. 4. It should be emphasized that the mobility is relatively large; the magnitude of which (2.2 cm2/Vs for Bi2(Sr,La)2CaCu2O8+δ and 1.8 cm2/Vs for Bi2(Sr,Ca,La)3Cu2O8+δ) is only 5 times smaller than that of Bi2Sr2CaCu2O8+δ (10 cm2/Vs) maeda , although the magnitude of the in-plane resistivity of Bi2(Sr,La,Ca)3Cu2O8+δ is 100 times larger than that of Bi2Sr2CaCu2O8+δ. In our previous work of the parent insulator Bi2Sr2RCu2O8+δ (R = Dy, Er), the mobility was too small to measure the Hall coefficient precisely syuuron . A similar large mobility at 300 K was reported in the lightly doped LSCO, where only 1% of hole doping causes the metallic behavior ando .

In the case of a conventional semiconductor, the mobility does not decrease with decreasing carrier concentration. It is rather possible that the mobility becomes large with decreasing carrier concentration because of the suppression of the electron-electron scattering. While, in the case of the doped antiferromagnetic insulator, such as high-$`T_c`$ superconductor, the hole is replaced by the spin in the next site whenever it moves. Then, if the hole moves independently, there is the energy loss of the order of the antiferromagnetic correlation energy and the mean free pass becomes as small as the lattice constant. The fact that the mobility of the parent insulator does not so much differ from that of an optimally doped sample does not agree with above simple consideration. It strongly suggests that the doped holes form a charge order, such as charge stripe or phase separation, and move in a group. Thus, the doped hole in the parent insulator is essentially itinerant and shows metallic behavior.

## IV Conclusion

In conclusion, we have successfully grown large single crystals of the parent insulator of Bi2Sr2CaCu2O8+δ by substituting La for Sr. With using these crystals, we have measured the in-plane resistivity , thermopower and Hall coefficient along the CuO2 plane. the in-plane resistivity of the most conducting sample shows metallic behavior above 150 K. Although the magnitude of the in-plane resistivity of the parent insulator is about 100 times larger than that of optimum doped Bi2Sr2CaCu2O8+δ, the mobility is not so much different from that of the superconducting compounds. These results suggest that the system is not a conventional insulator with a well-defined gap, and the hole in the parent insulator is essentially itinerant, as soon as it is doped. Our single crystal shows relatively large mobility compared to Ca-site-substituted single crystals, which indicates a smaller disorder in CuO2 plane.


text:
# The Magnetic Field Geometry in DR21

## 1 Introduction

The driving mechanism for the bipolar outflows so often seen associated with Young Stellar Objects (YSOs) is still a subject of considerable debate. Königl (1999) reviewed efforts to extend the commonly invoked centrifugally driven wind models of bipolar outflows in low-mass stars to high mass YSO’s. In these models the magnetic field is twisted into a spiral pattern rising vertically off of the star’s dense, rotating circumstellar disk. Material entrained in this field is forced up the vertical axis and produces the observed collimated outflow (Pudritz & Ouyed, 1998; Kudoh & Shibata, 1997).

Determining the magnetic field geometry in bipolar YSO outflows is a difficult observational task. These outflows are most easily studied when the outflow axis lies in the plane of the sky and the star’s equatorial disk is viewed nearly edge-on. This projection produces the classic bow-tie morphology seen in so many sources and allows both poles of the outflow to be easily separated spatially on the sky. If dust grains in the bipolar lobes are aligned by the magnetic field in the outflow, then polarimetry of starlight shining through the outflow can measure the projected magnetic field geometry. If most of the material in the flow contains a spiral magnetic field with a high pitch angle, then we may be able to observationally distinguish the magnetically driven models from a model that invokes some other flow mechanism in which the field is simply dragged along radially outward.

Since these outflows can be dusty and usually occur in star forming regions, infrared polarimetry of background starlight shining through the outflow can more easily reveal the composite magnetic field geometry along that line of sight than optical polarimetry. Traditional aperture photopolarimetry of field stars can be compromised by contamination from scattered light (reflection nebulosity) in the beam. By using imaging polarimetry, a diffuse component in the images can be effectively removed from ’under’ the star’s image. Imaging polarimetry also has the advantage that an image will usually contain several stars that can be measured. Traditional aperture polarimetry in the red (e.g. R and I bands) has been done for several sources and generally shows alignment of the local magnetic field with the outflow direction (Hodapp, 1990, cf.).

A related technique has been employed by Itoh et al. (1999, hereafter ICBHT). They used narrow band imaging polarimetry of knots of molecular hydrogen emission in the massive outflow associated with DR21 to probe the projected magnetic field geometry there. DR21 is one of the most massive and powerful bipolar outflows from young stars known. As such, it is an ideal laboratory for investigating the origin of the outflow. Their technique has the advantage of using emission that is from within the outflow itself, uncontaminated by transmission through background dust. However, it may suffer from contamination by underlying reflection nebulosity that could be hard to remove.

In this paper we report imaging polarimetry of DR21 in the broadband K (2.2$`\mu `$m) filter. We observe both background field stars shining through the outflow and reflection nebulosity in the flows. Note that background starlight samples all of the dust along the line of sight through the outflows, so we can not distinguish between material entrained by the outflow in a surrounding sheath and material interior to the flow itself.

## 2 Observations

The observations were made using NSFCAM on the IRTF in polarimetry mode. The basic technique is described in some detail in Jones (1997). Images are taken of the source and an off-source sky frame at four positions of a half waveplate that can be rotated in the beam. This allows us to form sky subtracted images at polarization position angles of 0, 45, 90 and 135. From these we computed Stokes images Q ($`I_0I_{90}`$) and U ($`I_{45}I_{135}`$) and the total intensity image $`\frac{1}{2}(I_0+I_{45}+I_{90}+I_{135})`$. Details of the flat-fielding and the image sequencing are given in Jones (1997).

Since NSFCAM is not a true imaging polarimeter, but rather a camera retrofitted with some polarization optics, we were limited in the precision of our polarimetry to about $`\pm `$0.3% during our observations of DR21. This limit was due primarily to fluctuations in transmission and sky background. We know this to be the case because the minimum error can be improved by shortening the time between waveplate positions. Unfortunately to do so dramatically reduces our duty cycle and integration time on source. Since the nebulosity and stars in DR21 are faint, we kept the duty cycle at 80%, which corresponds to about 12 seconds between waveplate positions.

Calibration of the efficiency and the position angle was done using observations of S1 in $`\rho `$ Oph with the same observing procedure as for DR21. S1 $`\rho `$ Oph was assumed to have a polarization of P = 1.95% and $`\theta `$ = 28 at K (2.2$`\mu `$m). Instrumental polarization was checked by observing unpolarized standard stars from the UKIRT faint standards list and was found to be too small to be measured. These measurements also established our limiting precision for that night.

Seeing was about 1.1” during our DR21 observations. Images were centered at 3 different locations in DR21. The Central image was centered on the bright stellar cluster at the heart of the DR21 complex. The approximate positions of the Eastern (E), Western (W) and Central images are given in Table 1. The sky background images for each of these positions were directly 2 arc minutes South. Since DR21 lies in a crowded region, the sky background images contained a significant number of field stars (but no nebulosity). These will show up in our images as black spots, since they produce negative values in the sky subtracted image. Although these stars could have been first ’removed’ from the sky frames for cosmetic purposes, to do so would have caused confusion in our polarimetric data analysis. This is because it is unlikely we could have removed the effects of stars in the sky frames to a level of $`\pm `$0.3% in polarimetric precision. By keeping the negative star images in the reduced data, we could more easily avoid measuring polarization in areas close to a subtracted field star from the sky image.

We measured the polarization of both the nebulosity associated with DR21 and a number of stars in the images. Finding charts for the stars and locations of the measured nebulosity for the Eastern and Western lobes are shown in Figure 1. The measured polarizations of these objects are given in Tables 2 and 3. Results for the central image are given in Tables 4 and 5 where the positions are the offsets in arc seconds relative to the brightest source, IRS1, near the bottom center of the core cluster.

The stars in all images were measured by using a set of concentric circular synthetic apertures with diameters of 2.4, 3.0 and 4.8”. The area between the outer two apertures was used to determine the nebular background which could then be subtracted from ’under’ the image of the star in the inner aperture. This allows us to remove any polarized nebulosity from the stellar measurements, and is a major advantage of imaging polarimetry over aperture photopolarimetry. Results for 8 stars in the lobes are given in Table 2 and illustrated in Figure 2. Results for several stars and point-like objects in the central cluster are given in Table 4 and illustrated in Figure 4.

Polarimetry of the nebulosity in the lobes was done using a 10” square synthetic aperture. The results are given in Table 3 and shown in Figure 3. Note that we only have upper limits for locations in the Western lobe, which clearly has less polarized nebulosity than the Eastern lobe. Polarimetry of the nebulosity in the Central image was done using a 2.7” square synthetic aperture. The results are given in Table 5 and illustrated in Figure 5. There were other locations in the images where we attempted to measure polarization but the nebulosity was too faint to yield any detections or any meaningful (5%)upper limits.

## 3 Discussion

### 3.1 Interstellar Polarization

ICBHT thoroughly discuss the issue of foreground interstellar polarization in the direction of DR21. Foreground polarization, if large, can contaminate the polarization due to the dust associated only with DR21. To first order, this foreground polarization can be removed if it is known, but this is rarely the case. Despite a substantial effort, ICBHT were unable to satisfactorily determine the foreground polarization from their analysis of surrounding field star polarimetry from the literature. The surrounding field star polarimetry shows a rather chaotic distribution and is not very useful in determining the foreground polarization.

The strength of the foreground polarization is probably low, however, based on the magnitude of the polarization of nearby stars. ICBHT estimate a magnitude of 0.8% at K attributable to foreground polarization. Unfortunately these optically visible field stars are all quite distant from DR21 on the sky, typically a degree or more. There are a few stars near the central cluster that Davis & Smith (1996) identify as optically visible stars. These stars are probably the best candidates for foreground stars with only foreground polarization that are accessible in our images. Two of these are designated I4 and I5 by Davis & Smith (1996) and are well enough separated from other sources for us to determine their polarization (Table 4). For both of these stars we have only an upper limit to the polarization of 0.9% (3$`\sigma `$). These measurements support the conclusion by ICBHT that the foreground interstellar polarization is rather low.

The position angle of this weak (if any) foreground polarization is still undetermined. ICBHT investigated two possible cases for subtracting an interstellar polarization component, one parallel to the outflow and one perpendicular to the outflow. Neither case altered their conclusions. Since the position angle of the foreground polarization can not be determined by our observations nor by ICBHT, and the magnitude of the polarization is certainly less than 0.9%, we will ignore it in this paper. We present other evidence later that further indicates foreground interstellar polarization is insignificant.

### 3.2 Stars in the Central Cluster

The polarization vectors for most stars in the central cluster (Figure 4) are in fair agreement with the bright central source, IRS1. Two stars, Nos. 8 and 9, have position angles more consistent with scattering than extinction and may be knots of reflection nebulosity. These sources are prominent in the broadband image taken by Davis & Smith (1996), but not their continuum subtracted H2 image. This means that these two sources are not knots of H2 emission but must be either stars or reflection nebulosity.

The position angle of the polarization of IRS1 (see Figure 2) is very nicely perpendicular to the 1.3mm polarization measured in a 19” beam by Glenn, Walker & Young (1999). The millimeter polarization is due to emission of polarized light from aligned dust grains. If the same aligned grains emitting the mm flux are also causing the interstellar polarization in transmission at K, then the position angles would naturally be perpendicular to one another (Hildebrand, 1983). Since the stars in the DR21 cluster are the source of energy for heating these grains, the mm emission is coming from dust local to DR21. Cold dust along the line of sight behind and in front of DR21 does not contribute to the mm polarization. Since our Near Infrared polarization of IRS1 has a position angle so close to perpendicular to the mm polarization, the extinction to DR21 must be dominated by grains local to the cluster, not in the foreground (unless the foreground interstellar polarization is at exactly the same position angle, of course).

This strongly suggests that at least for the central cluster, the grains causing the polarization we measure must be associated with DR21 itself. We can not make the same statement for the polarization in the lobes. The mm flux drops rapidly moving off from the central cluster and there are no published polarization measurements of thermally emitting dust in the lobes. This is unfortunate, because our probe of the magnetic field geometry in the lobes hinges on having the major fraction of the extinction to background stars associated with the lobes themselves.

Based on their narrow band imaging, Davis & Smith (1996) suggest there is evidence for a second bipolar outflow from the IRS1 region with a North-South axis. We note that star No. 7 has a position angle that lies along the ridge of reflection nebulosity to the North of IRS1 and is suggestive of having a N-S component to the polarization in addition to the polarization seen towards IRS1. Whether or not this is evidence for a N-S flow with a magnetic field aligned in that direction is only speculation at this point. There are a number of stars further North than our observations could reach that will be the subject of future observations.

### 3.3 Nebulosity in the Lobes

The nebulosity in our bandpass consists primarily of scattered continuum light (presumably from the central cluster) and local H2 emission. We were able to measure polarization of the nebulosity at 5 locations in the Eastern lobe (Figure 3). All five polarization vectors are perpendicular to the line of sight to the central cluster, entirely consistent with simple reflection nebulosity if the illuminating source is IRS1 or some combination of stars in the central cluster. This implies that there is a relatively clear, extinction-free path between the central cluster and the dust in the Eastern lobe.

Perhaps surprisingly, we were unable to detect polarization in the Western lobe’s nebulosity at the same level as seen in the Eastern lobe. Although our upper limits are only 4.2 and 6% (3$`\sigma `$) in the Western lobe, polarizations as high as 21% are present in the Eastern lobe. This suggests that the line of sight from the central cluster to the Western lobe may not be as clear and extinction free as in the East, and that most of the emission is from molecular hydrogen, relatively uncontaminated by reflection nebulosity.

### 3.4 Stars in the Lobes

There are 8 stars in the lobes well isolated and bright enough for us to measure their polarization. All 8 stars have similar polarization position angles, basically aligned with the axis of the outflow (Figure 2). There is no indication of the presence of any magnetic field geometry along the line of sight to these stars other than aligned with the outflow. Note that the polarization is due to all of the dust along the line of sight, and our measurements represent a complicated weighted average of the effects of aligned grains along the line of sight (Jones, Klebe & Dickey, 1991, cf.).

It is hard to see how this alignment of the polarization vectors with the outflow axis could be by chance. Even though we can not prove these stars are actually behind the outflow or that dust causing the polarization is primarily associated with the outflow itself, the data strongly suggest that it is so. There are three lines of circumstantial evidence that lead us to conclude that the polarization of the stars in the lobes is due to the local projected magnetic field geometry in the DR21 outflow. These arguments are:

1. The polarization of stars in the central cluster is clearly from dust local to the DR21, not foreground dust (§3.2).
2. The nebulosity in the Eastern lobe shows no significant deviation from a simple scattering geometry even when the polarization strength is as low as 4.8% (EN3), further indicating little contamination from foreground interstellar polarization.
3. Interstellar polarization caused by dust not associated with DR21 (foreground and/or background) would have to be aligned with the outflow axis by chance. There is no indication in the surrounding field stars for such an alignment (ICBHT, §3.1).

ICBHT avoid some of this confusion by measuring the polarization of the H2 emission in the outflow through a narrow band filter. This type of emission is local to the lobe, likely not scattered light from the central cluster, and therefore represents an intrinsically unpolarized source in the outflow itself. Consequently, there can be no contamination by dust from behind the lobes, only foreground interstellar polarization would be a concern (which we have argued is not a significant contributor). Our results for stars in the field are in excellent agreement with ICBHT in the Eastern lobe, where they also find polarization vectors aligned with the outflow axis. Clearly, the magnetic field geometry in the Eastern lobe is well aligned with the flow.

In the Western lobe, or results disagree with ICBHT, who find a polarization position angle perpendicular to the flow axis. Since we are observing stars, we can remove contamination from nebulosity in the star image (§2). This makes our technique relatively less sensitive to polarized reflection nebulosity. ICBHT, however, were measuring nebular emission, and it is possible their results in the Western lobe were contaminated by reflection nebulosity in their beam. Scattering will, of course, produce a polarization perpendicular to the line-of-sight to the illuminating source, which is what they found. However, our inability to find highly polarized nebulosity through a broad band filter in the Western lobe would argue against this explanation. Either the ICBHT results are contaminated by reflection nebulosity in their beam, or there is a real difference between the magnetic field geometry in front of the H2 knots they measured and in the dust along the lines of sight to the 3 stars we measured in the Western lobe.

We find no evidence for a tightly wound spiral magnetic field in either lobe of the DR21 outflow. This does not necessarily mean that no such field geometry exists, it is just that the majority of volume associated with the extinction (dust) in the outflow does not contain such a geometry. It is possible that a spiral magnetic field is still responsible for the outflow and is present in the inner portions of the flow, but that there is relatively little dust optical depth associated with the volume occupied by the spiral field. Large amounts of dusty material entrained by the flow on its outside boundary could easily produce a cocoon around the flow that exhibits a magnetic field aligned with the outflow. More theoretical work needs to be done to determine if this could be the case.

## 4 Conclusions

We have measured the polarization of stars and nebulosity in the field of the DR21 stellar cluster and in the Eastern and Western lobes of its massive bipolar outflow through a broadband K filter. Based on our observations we can reach the following conclusions:

1. If the stars are shining through the outflow from behind and there is minimal foreground and background interstellar polarization, then the magnetic field geometry sampled by the dust in the outflow is well aligned with the outflow axis. There is no evidence for a tightly wound spiral magnetic field nor significant fluctuations or turbulence in the field geometry.
2. The polarization vector of the bright central source IRS1 is nicely perpendicular to the polarized millimeter emission observed towards the same source. This strongly implies that the dust causing the K band polarization for IRS1 is local to DR21, not in the foreground.
3. Through a broadband filter, the nebulosity in the Eastern lobe has significant polarization indicative of simple reflection nebulosity illuminated by the core cluster. Oddly, we were unable to measure polarization this strong in the nebulosity in the Western lobe. These results suggest that a clearer line of sight exists between the nebulosity in the Eastern lobe and the central illuminating cluster than is the case in the Western lobe.
4. Itoh et al. (1999) have made observations complimentary to ours. Our results agree with theirs in the Eastern lobe where we both find position angles aligned with the outflow axis. Our observations disagree with Itoh et al. (1999) in the Western lobe where they find a position angle perpendicular to the outflow axis. Either the Itoh et al. (1999) results are contaminated by reflection nebulosity in their beam, or there is a real difference between the magnetic field geometry in the H2 knots they measured and the dust along the lines of sight to the stars we measured in the Western lobe.

## 5 acknowledgments

We would like to thank the support staff of the IRTF for their assistance helping make imaging polarimetry a viable option on the IRTF.


text:
# Eccentricity Evolution in Simulated Galaxy Clusters

## 1 Introduction

Galaxy clusters are the largest bound objects in the Universe. They provide information on gravitational instability, and understanding their formation is helpful in describing the transition to nonlinearity in structure formation.

In studies of a variety of optical and X-ray samples, Melott, Chambers, and Miller (2001, hereafter MCM) and Plionis (2002) found significant evolution in the gross morphology of galaxy clusters over the redshift range 0 to 0.13. When projected galaxy clusters were fit to ellipses by a variety of procedures, their eccentricities were found to evolve on a timescale comparable to one or two crossing times for the clusters. The inclusion of X-ray data is important, as it is much less subject to projection effects than optical data. In all cases, with varying significance, the clusters were found to have become less eccentric at lower $`z`$. It was argued that this could be understood as a relaxation process in a low-density cosmological background, made possible by a greatly reduced recent merger rate. Mergers tend to occur by infall along supercluster filaments (Shandarin and Klypin 1984; Colberg et. al 1999) which can lead to alignment of the cluster axis (Bingelli 1982; Chambers, Melott, & Miller 2002, and references therein) as well as aligned bulk flows within the cluster (Novikov et al. 1999).

This poses new questions of structure formation for cosmological simulations. Will cluster eccentricity evolution appear there? Does it depend on the initial conditions, the background cosmology, and/or the physics included within given simulations?

Richstone, Loeb, & Turner (1992) showed that cluster morphology could be used as a diagnostic of mass density in the Universe. While a cosmological constant, $`\lambda `$, was argued to have little effect, a large mass density $`\mathrm{\Omega }_m`$ would produce more irregular clusters due to recent formation; on this basis they argued $`\mathrm{\Omega }_m0.5`$. Mohr et al. (1995) looked at a variety of simulations, finding their $`\mathrm{\Omega }_m=1`$ simulations closest to the distribution of axial ratios for observed clusters, since low-$`\mathrm{\Omega }_m`$ simulations produced too-spherical clusters. Jing et al. (1995) disagreed, noting that open or $`\lambda `$-dominated flat cosmologies also produced substantial irregularity in clusters. Using the terminology OCDM and $`\lambda `$CDM for these, and $`\tau `$CDM for the $`\mathrm{\Omega }_m=1`$ case (with the same initial conditions), they found increasing irregularity and eccentricity in the order OCDM, $`\lambda `$CDM, $`\tau `$CDM. They found the effect of $`\lambda `$ was smaller for increasing $`\mathrm{\Omega }_m`$, as one would expect. They did not explore evolution. Thomas et al. (1998) reinforced these conclusions. Suwa et al. (2003) presented results on ellipticity evolution, comparing clusters only at $`z=0`$ with those at $`z=0.5`$ (but not intermediate redshifts). These results suggest that OCDM clusters are slightly rounder than $`\lambda `$CDM, although again the difference is small and not significant in their study. They found a small amount of ellipticity evolution between the two redshifts.

We are interested in exploring the evolution of eccentricity in the recent past, especially for redshifts $`0z0.1`$, where there is a variety evidence for evolution in the data. For this reason we have emphasized the outer regions of clusters, and used methods similar to those applied to many observed clusters.

## 2 Simulations

We examined projected cluster properties in a variety of simulations based on Cold Dark Matter and Gaussian perturbations (hereafter CDM). One set has been previously described and extensively used elsewhere (Loken et al. 2002, and references therein). We refer to this as the $`\lambda `$CDMH simulation, to distinguish its inclusion of hydrodynamics from other, purely N-body simulations we studied. This simulation focused on the most massive objects within a region 256 $`h^1Mpc`$ on a side in a flat, cosmological-constant dominated universe ($`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Omega }_\lambda =0.7`$) with a baryon fraction of $`\mathrm{\Omega }_b=0.026`$ and a Hubble constant of 70 km/s/Mpc. The initial conditions used the Eisenstein & Hu (1999) parameterization of the CDM power spectrum, normalized to an amplitude corresponding to $`\sigma _8`$=0.93, where $`\sigma _8`$ is the rms density fluctuation inside spheres of radius 8 $`h^1Mpc`$. The dark matter particle mass was $`1.3\times 10^{10}M_{}`$ and the evolution of the gas was followed with an Adaptive Mesh Refinement method (see Bryan 1999) that solves the equations of hydrodynamics (not including radiative cooling, thermal conduction, or star formation) on a grid. The method starts off with a relatively coarse mesh but decreases the mesh spacing in dense region so that the smallest cell size (best resolution) achieved was 16 $`h^1`$ kpc. We used X-ray images assuming a metallicity of 0.3 of solar in the 0.5-2.0 keV bandpass. We used previously saved data on 11 clusters taken at $`z`$ = 0.25, 0.1, and 0, identified using the “HOP” algorithm (Eisenstein and Hut 1998, hereafter EH). This algorithm finds regions above some threshold density, and merges them, associating them with regions of some maximal density $`\delta _{peak}`$. These regions are merged into a single region (“cluster”, for our purposes) if they lie within a contour of density $`\delta _{saddle}`$. Any associated particles are also considered part of the cluster if they lie within a connected contour above the density $`\delta _{outer}`$. The parameters used to identify these clusters were were 480, 400, and 160 for $`\delta _{peak}`$, $`\delta _{saddle}`$, and $`\delta _{outer}`$, respectively. For more details on this method see EH. In the $`\lambda `$CDMH, the mesh-refinement method requires following the same volume element, chosen by the presence of a cluster at $`z=0`$, through its evolution. Clusters at higher redshift are taken from this volume; it is possible that in this case they may not be the most massive clusters in the large simulation volume at earlier times. Another substantial difference with the larger cluster sample which will be analyzed next is that these very massive clusters correspond to a Richness in excess of 2. In fact they are even more massive than typical $`R=2`$ clusters, but less than $`R=3`$ clusters. It is not possible to change this given the AMR method, although checking the conclusions against other (N-body) simulations constitutes a partial cross-check.

We conducted another, larger, ensemble of pure Nbody simulations. All used the same initial power spectrum, corresponding to CDM (Bardeen et al. 1986) with $`\mathrm{\Gamma }`$=$`\mathrm{\Omega }`$h=0.225. This was evolved in both an open ($`\mathrm{\Omega }_{m0}`$=0.34, OCDM1), a flat matter-dominated ($`\mathrm{\Omega }_{m0}`$=1, $`\tau `$CDM), and a flat cosmological-constant model ($`\mathrm{\Omega }_{m0}`$=0.34, $`\mathrm{\Omega }_{\lambda 0}`$=0.66, $`\lambda `$CDM1) Inclusion of the $`\tau `$CDM model was motivated by the wish to see whether major differences in the expansion rate affected the eccentricity. All three of these ensembles were run as before, to an amplitude corresponding to $`\sigma _8`$=0.93. A Hubble Constant $`h=2/3`$ was used in the simulation analysis. These N-body runs had three realizations each of $`256^3`$ particles in boxes of 128 $`Mpc`$ diameter. Given the volume and target space density of $`R1`$ clusters, 6 clumps were extracted from each of these three simulation ”boxes”. For consistency with our previous observational analyses, we simulated clusters of Richness $`R1`$. Volume-limited samples exist for such clusters, but not for Richness 0. We used parameters 400, 200, 80 for HOP in this case, in analogy with densities of observed clusters and virialization density thresholds. Having identified clumps with HOP, clusters were identified with the $`N`$ most massive clumps in the simulation corresponding to a mean separation of $`50h^1`$ Mpc. We therefore have a total of 18 clusters from our three realizations. Given two independent axial ratios out of three projections, we have 54 views at each redshift which have the statistical weight of 36 views. In other words, this simulation ensemble produces the equivalent of 36 independent observed clusters for each assumed background cosmology. (We have the equivalent weight of 22 views in the $`\lambda `$CDMH runs.) The dark matter particle mass was $`5.2\times 10^9M_{}`$ in the low $`\mathrm{\Omega }_m`$ model and proportionally larger in the high $`\mathrm{\Omega }_m`$ model. Data was taken at moments corresponding to $`z`$= 0.2, 0.1 and 0 in the OCDM1 and $`\tau `$CDM simulations.

In the OCDM1 boxes, our selected objects ranged from about 8 $`\times `$ $`10^{13}`$ to 1 $`\times `$ $`10^{15}`$ M, for the mass within an Abell radius. This agreement with the mass range of observed clusters is a consistency check on our procedure. In the $`\tau `$CDM box, they were about three times more massive, as expected for consistency with the assumed background cosmology. In the $`\lambda `$CDMH box, the range was from about 5 $`\times `$ $`10^{14}`$ to 2 $`\times `$ $`10^{15}`$ M, appropriate for $`R>2`$ clusters.

Having chosen these clusters, three Cartesian projections were constructed. There are only two independent axial ratios, $`a/b`$, and $`b/c`$, therefore only two independent ratios in projection, whether we project along cluster principal axes or not. We project along the simulation coordinate axes, since we do not want to prefer principal axes; we use all 3 axes for balance. In order to take account of this partial lack of independence, we recognize that the standard deviation of our results is smaller by a factor of $`\sqrt{\frac{2}{3}}`$ than if they were fully independent. This is taken into account in all statistical results quoted here.

After initial analysis of OCDM1, we decided to examine more closely the possible effect of amplitude normalization. We therefore reran the entire ensemble with lower values of $`\sigma _8`$=0.78 and 0.58, called OCDM2 and OCDM3 respectively. Naturally, these clusters are somewhat less massive than in OCDM1. Likewise, $`\lambda `$CDM2 has $`\sigma _8`$=0.83.

We have 33 projections of $`\lambda `$CDMH simulations at each redshift, and 54 of each of the OCDM\*, $`\lambda `$CDM\*, and $`\tau `$CDM simulations. Including the fact that each of these is viewed at three redshift stages, our simulated cluster sample size is larger than any of the individual data samples studied in MCM, though not as large as the APM sample used by Plionis (2002).

Fair sampling is a consideration here. We are viewing the same overall simulation volume, and in most cases, the same clusters, at a variety of redshifts. This is obviously not true in the observed samples. If this were important, cosmological simulations as currently analyzed would be of limited relevance to understanding structure formation. Separate simulations would be needed for each redshift bin. We partially compensate for this effect (which we assume to be small) in the OCDM\*, $`\lambda `$CDM\*, $`\tau `$CDM simulations by identifying clusters separately in each redshift bin. This occasionally resulted in the selection of different objects at different stages of evolution.

In the $`\lambda `$CDMH, the mesh-refinement method requires following the same volume element, chosen by the presence of a cluster at $`z=0`$, through its evolution. Only the most massive clusters in this simulation were chosen for mesh refinement. The space density and threshold mass for the objects chosen corresponds to a richness $`R>2`$. Clusters at higher redshift are taken from this sub-volume; it is possible that in this case they may not be the most massive clusters at that $`z`$ in the large simulation volume. An additional caveat in the application of the $`\lambda `$CDMH group is therefore its correspondence to richer clusters. However, according to the study of Plionis (2002), ellipticity evolution is not weaker for richer clusters.

## 3 Determination of Eccentricity

Observers had widely varying procedures for selecting which portions of clusters to analyze, and in some cases the procedure was not specified. When major and minor projected axes were found, varying definitions of $`ϵ`$, the ellipticity, were applied. In fact, varying definitions of ellipticity appear in common use in mathematics, physics and astronomy literature. We followed one well-specified procedure and a standard definition of eccentricity, fitting ellipses to our projected clusters. This geometrical term is well-known, and not subject to the ambiguity in definitions of “ellipticity”. For an ellipse with major and minor axes $`a`$ and $`b`$ respectively, the eccentricity $`e`$ is defined by

$`e=\sqrt{1{\displaystyle \frac{b^2}{a^2}}}`$ (1)
this definition of $`e`$ is also used for $`ϵ`$ by some, but not all, work which reports on “ellipticity”. Another common definition of ellipticity, $`ϵ`$= 1-b/a, has $`e`$ and $`ϵ`$ increasing monotonically with one another inside their useful ranges from $`0`$ to $`1`$. Our statistical inferences, discussed later, are nonparametric and based on ranks; thus they would be the same whether we used $`e`$ or $`ϵ`$.

Only projected information is available for observed clusters in useful precision, so we study our simulated clusters in projection. The center of these objects was defined for most analyses here as the center of mass of the linked particles. It must be stated clearly that the procedures used in observational studies are not uniform, or even always clearly specified. This center of mass is our closest possible analog for the center of mass of the optically luminous matter–the galaxies. In the case of the $`\lambda `$CDMH simulation set, we also analyzed the synthetic X-ray emission of the clusters. In this case, we chose as a center the region of strongest simulated surface brightness in the X-ray. For the sake of fuller exploration, we repeated this procedure using the highest peak of the projected mass density as the center, and found no qualitatively different conclusions after referring our computation to the center of mass of regions drawn from around this center.

We selected the widely used “annulus method”; it emphasizes the outer regions of the cluster, which is important in trying to measure ellipticity for highly centrally concentrated X-ray images. This is the most widely used method in cluster ellipticity work. For each cluster projection, having chosen a center as described above, a circle of radius 1.5 $`h^1Mpc`$ was drawn about this point. Then another circle of smaller radius was placed so that the annulus between them had 0.2 times the entire mass enclosed within the outer circle. The purpose of using annuli is to emphasize outer regions, especially important for highly concentrated, emission-weighted X-ray images. We repeated this procedure using 2.0 $`h^1Mpc`$, with no qualitatively different conclusions. We also examined the eccentricity of our clusters when the entire mass within 1.5 $`h^1Mpc`$ was used rather than an annulus. Although eccentricity values were generally lower, our conclusions again were qualitatively the same. We further tested the robustness of our conclusions against procedure by dropping the annulus procedure, and computing the values based on the entire region connected by the HOP procedure. This also shifted the values somewhat, but did not change our conclusions about evolution. The results quoted herein are based on the 1.5 $`h^1Mpc`$ annulus.

We computed the “inertia tensor” $`𝐈`$ of this annulus with respect to its center of mass (center of emission for X-ray images). Very few observational studies mentioned that they used the center of mass frame; we hope that they did since the results are nearly meaningless without. If not computed with respect to the center of mass, the eigenvalues have no particular relationship to the axes of the object. We used the conventional definition

$`I_{ij}=\mathrm{\Sigma }m_\alpha (\delta _{ij}r^2x_{\alpha i}x_{\alpha j})`$ (2)
where we sum over regions $`\alpha `$.

The eigenvalues of $`𝐈`$ are proportional to $`a^2`$ and $`b^2`$, as defined above, for a homogeneous ellipse. In this way we define the eccentricity of the projected simulated cluster approximated as an ellipse.

## 4 Results

In Table 1 below, we show for for the $`\lambda `$CDMH simulation, the median and mean $`e`$ at each redshift bin, as well as the standard deviation of the mean. Note that the standard deviation in the estimate of the mean is smaller than the standard deviation of one measurement. Figures 1 and 2 contain histograms of the eccentricity for the OCDM1, $`\tau `$CDM, and $`\lambda `$CDMH ensembles. Note that in $`\lambda `$CDMH, we show eccentricity inferred both from the dark matter distribution, which is a stand-in for the optical studies, assuming that galaxies reasonably trace the dark matter; we also show results derived from the inferred X-ray luminosity. There are weak evolutionary trends, but none as steep as that seen in the data. Again, note that this data sample corresponds to very rich clusters.

In agreement with previous studies, we find weakly increasing eccentricity in the order OCDM1, $`\lambda `$CDM1, $`\tau `$CDM at $`z=0`$.

We also performed a Wilcoxon rank-sum test (Lehman 1998) to compare samples at different stages of evolution and with one another. A nice feature of this test is that it will give the same results for comparing measures which are nonlinearly related, but monotonically increasing functions of one another, like $`e`$ and $`ϵ`$. In this case, the null hypothesis is that the lower-$`z`$ sample does not have systematically smaller eccentricity than the high-$`z`$ sample. In the last column, we give the confidence in rejection of eccentricity evolution (reduction with lower $`z`$). In particular, the number in this column is the probability that the change in eccentricity seen in the simulations could result from random draws of a cluster population which did not have decreasing eccentricity. The first number given in the last column is the confidence resulting from comparison with the cluster ensemble on that line with the ensemble on the following line (next higher redshift datafile). The following number (in parentheses) is the confidence resulting from comparison with the ensemble two redshift bins higher (i.e. between redshift 0 and 0.2 or 0.25) Small values of $`P_W`$ correspond to the direction of evolution seen in the data.

Although there is a significant signal in a variety of both optical and X-ray data samples that galaxy clusters have undergone a morphological evolution, becoming significantly more circular since $`z`$ = 0.1, we have not found this in our ensemble of simulations. There is such evolution only in the OCDM1 and $`\lambda `$CDM1 simulation sample for $`0.1z0.2`$. OCDM3 has evolution in the opposite sense over the same redshift interval. Evolution found in data by MCM was found only for $`z0.1`$; that found by Plionis was for $`z0.18`$ and in fact appeared concentrated in lower redshift ranges. It was confirmed (Plionis, personal communication) that no significant ellipticity evolution appears in that sample for $`z0.13`$. Despite the differences in codes and lesser differences in cluster selection procedures, none of the simulated samples showed strong eccentricity evolution over the crucial period $`0z0.1`$. Evolution similar to that found in most of the data would correspond roughly to eccentricity dropping from about 0.8 to 0.7 during that period, which would stand out with high significance in our simulations.

There is a slightly different situation with the APM cluster sample. The evolution rate is slower there; it would correspond to mean eccentricity changing from about 0.76 to 0.70 during the redshift interval 0.1 to 0. (extrapolating its slope). This is a dangerous extrapolation, since this sample does not contain information for $`z<0.038`$, and its rate of evolution appears to steepen for lower $`z`$. Such a steepening would be consistent with the more rapid rate of evolution found at lower redshifts in all but one of the MCM samples.

It should be mentioned that the APM cluster sample was selected in a somewhat different way than the controlled Abell/ACO sample used elsewhere. These clusters were initially found by selecting regions of a high overdensity in a smaller region, about half the diameter typically used in the Abell catalog (Dalton et al. 1997). This procedure might produce a selection effect accounting for the difference.

To summarize, the APM sample would appear to be consistent with the rate of evolution in many of our samples at about the 2$`\sigma `$ level, if extrapolated linearly to $`z`$=0.

## 5 Discussion

In some way, the simulations as we have analyzed them do not faithfully model the observed “local” Universe. There are a variety of ways to account for this; we will mention a few.

We initially used a normalization $`\sigma _8`$=0.93. There is a hint in the OCDM1 results that a lower normalization might produce the strong evolution seen in the data in the redshift range needed. We decided to test this, and the OCDM2 and OCDM3 results provide evidence against this hypothesis (see Table 1 and Figure 3). We also wished to investigate cosmic variance beyond the simple variation of random phases in the initial conditions. Density fluctuations on the scale of the N-body box should be no more than about 10% $`rms`$. As an aggressive test of cosmic variance, we reran OCDM1 in a background with $`\mathrm{\Omega }_{m0}`$=0.2. This also would be an effective test of any possible effects due to our being located in a very large underdense region, such as that suggested by Busswell et al. (2003) and Frith et al. (2003). We found no significant change in the results of the ensemble.

The most obvious possibility is observational selection effects. The sign of the effect needed is that more eccentric clusters are included in cluster surveys when they are far away than when they are close.

The data used in MCM and Plionis (2002) are based on optical cluster catalogs. In MCM, they used subsets of the Abell/ACO catalogs that had published ellipticities (in either X-ray or optical). Plionis (2002) used a subset of the Automated Plate Machine (APM) clusters as published in Dalton et al. (1997), with ellipticities as given in Basilikos, Plionis, and Maddox (2000). Plionis (2002) also looked at Abell/ACO clusters in their APM sample and found the same ellipticity evolution as for the entire sample.

In the optical and Xray samples examined in MCM, ellipticity appears to increase with redshift for all samples. This increase exists in all samples studied in both papers, but the way in which it appears is different for different samples: (1) a lack of nearby highly elliptical clusters (compared to higher redshift); (2) an increase in the number of highly elliptical clusters at higher redshift (compared to nearby); (3) a lack of highly circular clusters at high-$`z`$ (compared to low-$`z`$); and of course any combination of these. In fact, there are observed samples consistent with all of these possibilities (MCM, Plionis 2002). It is worthwhile asking whether selection effects (both at the low-$`z`$ and the high-$`z`$ end of the cluster samples) could play a role in this correlation.

For instance, it is possible that both the APM and the Abell/ACO catalogs are incomplete at low redshifts and that those clusters missing from the sample have significant ellipticity. Likewise, an incompleteness in circular clusters at higher redshifts is also possible. This latter case seems less plausible, since such “regularly” shaped clusters should be easy to identify. Ebeling et al. (1999) claimed that high levels of substructure will lead to a decreased probability of cluster selection. The selection function for regularly shaped (circular) clusters should be stable over the range of redshifts examined in MCM and Plionis et al. In this case it is doubtful that circular clusters would be preferentially missed at higher redshifts.

However, Adami et al. (2000) conducted extensive testing on their X-ray cluster selection function (for the Serendipitous High-Redshift Cluster (SHARC) survey (Romer et al. 2000)). They reported that highly elliptical clusters were more easily found at high redshift. Adami et al. note that for a fixed cosmology, an elliptical cluster will have a higher central surface brightness than a circular cluster with the same luminosity. They found that this effect was only significant at the high-$`z`$ limit of the SHARC survey, where the selection function is very sensitive to the central surface brightness of a cluster. All redshifts simulated by Adami et al. were greater than the ones considered here. In spite of this objection, it is possible that elliptical clusters could be more easily detected at high-$`z`$ and that some highly circular clusters could be missed at high-$`z`$. Since all of the samples examined in MCM and Plionis et al. were optically selected, one would have to assume that an analogous effect could be present in optical samples.

In general, it remains unclear how ellipticity affects the selection of clusters. The argument made by Ebeling et al. (1999) that clusters with substructure are harder to detect, need not necessarily translate into elliptical clusters being harder to detect. The use of extensive simulations (as in Adami et al. 2000) is probably a requirement before any conclusions can be drawn regarding a ellipticity evolution in any given dataset. Unfortunately, the selection procedures used for the cluster samples in Plionis (2002) and MCM are so complicated that their selection functions are too difficult to model. Our best hope is to examine the next generation of cluster catalogs (e.g. from the Sloan Digital Sky Survey) whose completeness and selections should be well-measured.

It is possible that the method of defining included mass and computing eccentricity could affect these results. However, we also explored using a friends of friends algorithm, as well as various HOP parameters, and abandoning the annulus method in favor of inclusion of all mass in the clusters. None of these changed the evolution of eccentricity, though some changed its amplitude (consistently over redshifts, so as to preserve the lack of evolution).

Another possibility is that the simulations are missing crucial physics which might induce this change. Although violent relaxation may isotropize objects through the time-dependent gravitational potential, we cannot exclude processes in the hot gas which is the dominant baryonic component of the clusters. Radiative cooling may cause the gas to contract, gravitationally entraining some of the dark matter. Also, it has been shown (Narayan and Medvedev 2001) that thermal conduction may be important, especially in the turbulence following major mergers.

## 6 Acknowledgments

We thank M. Plionis and G. Evrard for useful comments. SNF and ALM gratefully acknowledge the support of the National Science Foundation through grant AST-0070702, especially a supplement for Research Experiences for Undergraduates. Invaluable computing support came from the National Center for Supercomputing Applications.