text:
# Density of states of a binary Lennard-Jones Glass

## I Introduction

The transition from a liquid to an amorphous solid that sometimes occurs upon cooling remains one of the largely unresolved problems of statistical physics Ediger et al. (1996); Debenedetti and Stillinger (2001). At the experimental level, the so-called glass transition is generally associated with a sharp increase in the characteristic relaxation times of the system, and a concomitant departure of laboratory measurements from equilibrium. At the theoretical level, it has been proposed that the transition from a liquid to a glassy state is triggered by an underlying thermodynamic (equilibrium) transition Mezard and Parisi (1999); in that view, an “ideal” glass transition is believed to occur at the so-called Kauzmann temperature, $`T_K`$. At $`T_K`$, it is proposed that only one minimum-energy basin of attraction is accessible to the system. One of the first arguments of this type is due to Gibbs and diMarzio Gibbs and DiMarzio (1958), but more recent studies using replica methods have yielded evidence in support of such a transition in Lennard-Jones glass formers Mezard and Parisi (1999); Coluzzi et al. (2000a); Grigera and Parisi (2001). These observations have been called into question by experimental data and recent results of simulations of polydisperse hard-core disks, which have failed to detect any evidence of a thermodynamic transition up to extremely high packing fractions Santen and Krauth (2000). One of the questions that arises is therefore whether the discrepancies between the reported simulated behavior of hard-disk and soft-sphere systems is due to fundamental differences in the models, or whether they are a consequence of inappropriate sampling at low temperatures and high densities.

Different, alternative theoretical considerations have attempted to establish a connection between glass transition phenomena and the rapid increase in relaxation times that arises in the vicinity of a theoretical critical temperature (the so-called “mode-coupling” temperature, $`T_{MCT}`$), thereby giving rise to a “kinetic” or “dynamic” transition Götze and Sjögren (1992). In recent years, both viewpoints have received some support from molecular simulations. Many of these simulations have been conducted in the context of models introduced by Stillinger and Weber and by Kob and Andersen Kob and Andersen (1995); such models have been employed in a number of studies that have helped shape our current views about the glass transition Sastry et al. (1998); Sciortino et al. (1999); Donati et al. (1999); Coluzzi et al. (2000a, b); Yamamoto and Kob (2000). The particular model considered here consists of a binary mixture of Lennard-Jones particles, with composition 80% $`A`$ and 20% $`B`$. A total of 250 particles is employed in our calculations. The interaction parameters between particles of species $`A`$ and $`B`$ are $`ϵ_{AA}=1.0`$ and $`\sigma _{AA}=1.0`$, $`ϵ_{BB}=0.5`$ and $`\sigma _{BB}=0.88`$, and $`ϵ_{AB}=1.5`$ and $`\sigma _{AB}=0.8`$. The density is $`1.204\sigma _{AA}^3`$. Recently, a crystal structure at extremely low energies has been reported for a variant of this system Middleton et al. (2001).

High-precision data are available for the thermodynamic properties of this model at intermediate to high temperatures Kob and Andersen (1995); Yamamoto and Kob (2000); Coluzzi et al. (2000b). A series of careful simulations have placed the mode coupling temperature at $`T_{MCT}=0.435`$ and the Kauzmann temperature somewhere in the range $`T_K=0.260.31`$ Yamamoto and Kob (2000); Sciortino et al. (1999); Coluzzi et al. (2000b). Note, however, that literature studies have generally avoided direct simulations below $`T_{MCT}`$; available estimates of $`T_K`$ have been produced after making several assumptions regarding the potential energy landscape and by extrapolation (to low temperatures) of liquid-state data generated at higher temperatures (above $`T_{MCT}`$). An exception is provided by a recent report for a related model Grigera and Parisi (2001), where simulations of small systems, directly at low temperatures, suggest that an anomaly in the heat capacity $`c_v`$ arises at $`T_K`$; $`c_v`$ is reported to increase with decreasing temperature, and to exhibit a sharp drop at $`T_K`$. The drop becomes more pronounced as the system size is increased.

Simulations near a glass transition are notoriously difficult, and their results must be considered with caution. On the one hand, the relevant time scales below $`T_{MCT}`$ are too long to be sampled by conventional molecular dynamics simulations. Monte Carlo techniques, on the other hand, have been used only rarely to simulate glass formers; furthermore, it has been difficult to establish to what extent available studies have succeeded in sampling relevant regions of phase space, particularly at low temperatures and elevated densities. In this work, we use a novel Monte Carlo sampling technique to arrive at direct estimates of the thermodynamic properties of a model glass former down to temperatures well below $`T_{MCT}`$.

## II Simulation Methods

Recently, Wang and Landau have proposed an iterative method to estimate the density of states of a Potts lattice system from a Monte Carlo simulation Wang and Landau (2001a, b). The random-walk algorithm is based on the idea of entropic sampling, with a self-consistent update of the density of states. It has proven to be remarkably efficient for lattice systems, simple liquids Yan et al. (2002), proteins Rathore and de Pablo (2002); Rathore et al. (2003), and liquid crystals Kim et al. (2002); it is tempting to apply it in the context of a glass-forming liquid. In this contribution we combine it with biased sampling techniques, and we use it to generate direct estimates of the density of states of the glass-former described above.

In a conventional canonical-ensemble simulation, different states of the system are visited with probability $`\mathrm{\Omega }(E)e^{E/k_BT}`$, where $`\mathrm{\Omega }(E)`$ is the density of states (or degeneracy) of the system, $`k_B`$ is Boltzmann’s constant, and $`T`$ is the temperature. In contrast, in the random-walk scheme adopted here, the density of states $`\mathrm{\Omega }(E)`$ is estimated directly by producing a uniform, or “flat” histogram of energies, i.e. by coercing the system to visit all energy states with equal probability. In this study we have chosen to maintain a constant density and constant number of particles; extensions to other physical ensembles and to expanded ensembles have also been pursued recently Yan et al. (2002); Kim et al. (2002); Calvo (2002). Trial moves are generated by means of simple translations of the particles and by identity interchanges Grigera and Parisi (2001). The acceptance of such interchanges is enhanced using configurational bias. The resulting trial configurations are accepted with probability Jain and de Pablo (2002)

$$p=\text{min}\{1,\frac{W_{new}^R\mathrm{\Omega }^{}(E_{old})}{W_{old}^R\mathrm{\Omega }^{}(E_{new})}\}$$
(1)
where the prime indicates that this is a transient, momentary “best estimate” of the density of states. Biased moves are performed according to a Rosenbluth type algorithm; $`W^R`$ is the Rosenbluth weight of the corresponding state. In the original version by Wang and Landau $`W^R`$ was set to unity for all states. The configurational bias identity swap consists of the following steps: First a pair on unlike particles is chosen at random. After direct interchange of their positions, the smaller particle $`B`$ will fit in the cavity formerly occupied by the larger particle $`A`$. The opposite is only seldom true, thereby leading to negligible acceptance rates at low energies and high densities. In order to enhance the acceptance rate, $`N_{CCB}=12100`$ trial positions are explored for the $`B`$ particle around the position formerly occupied by the $`A`$ particle. We then apply configurational bias ideas to these $`N_{CCB}`$ positions. There are now two possible ways to calculate the Rosenbluth factors. In the first of these, the energy of the states can be used directly (as is done in standard configurational bias Monte Carlo) de Pablo et al. (1992); Frenkel and Smit (1996). To this end, a fictitious temperature $`T_F`$ is introduced for calculation of the Rosenbluth factor $`W_i^R`$ for state $`i`$

$$W_i^R=\frac{\mathrm{exp}(\beta _FE_i)}{_j\mathrm{exp}(\beta _FE_j)}.$$
(2)
Note that this fictitious temperature is not the temperature of the system, although in conventional configurational bias simulations it is set to the system temperature. Alternatively, one can avoid using a fictitious temperature by calculating a set of $`W_R`$ using the density of states itself as a bias.

$$W_i^R=\frac{\mathrm{\Omega }_i}{_j\mathrm{\Omega }_j}.$$
(3)
Both variants (using $`T_F=0.51`$) behave similarly. The acceptance rate is extremely small (it drops to less than $`10^6`$ as the effective temperature approaches $`T_g`$), but it is sufficient to perform simulations over extended amounts of computer time. The effective temperature $`T_{\text{eff}}`$ of a state with energy $`E`$ is defined as the temperature where $`E_{pot}_{T_{\text{eff}}}=E`$.

The density of states is not known á priori; it is initially set to unity throughout the entire energy range. The calculations begin by defining an energy range in which to determine $`\mathrm{\Omega }(E)`$. Whenever an energy state $`E`$ is visited, the density of states $`\mathrm{\Omega }^{}(E)`$ corresponding to that energy is multiplied by a constant $`f`$, i.e. $`\mathrm{\Omega }^{}(E)f\mathrm{\Omega }^{}(E)`$. Since the density of states varies over many orders of magnitude, it is convenient to work with its logarithm, which corresponds to the entropy as a function of energy $`S(E)=k_B\mathrm{ln}\mathrm{\Omega }(E)`$. The entropy is updated by adding a constant $`\mathrm{ln}f`$.

A histogram of energies $`h(E)`$ is also constructed and updated after every trial move. The density of states is updated continuously throughout the simulation, until the recorded energy histogram $`h(E)`$ is sufficiently flat. Note that in actual practice it is not possible to generate a perfectly flat histogram; in this work, flatness is considered to be attained if the minimum of $`h(E)`$ is at least 0.9 times the average value $`\overline{h(E)}`$. Having reached a “flat-histogram” condition, the simulation sequence is repeated: the energy histogram is erased, and the new value of $`\mathrm{ln}f`$ is set to half the “old” value. Note that this choice is arbitrary, and any monotonically decreasing function should work. The initial value of $`\mathrm{ln}f`$ was set to unity, and the final value was $`\mathrm{ln}f=10^5`$, which corresponds to 18 iterations. The factor $`\mathrm{ln}f`$ controls the convergence of $`\mathrm{\Omega }^{}(E)`$ to the true value, $`\mathrm{\Omega }(E)`$; as $`\mathrm{ln}f`$ decreases (i.e. as the simulation proceeds) the calculations per iteration become increasingly long.

To improve efficiency, it is useful to conduct multiple simulations in overlapping energy ranges. These ranges must be relatively narrow; otherwise, given the rapidly varying nature of $`\mathrm{\Omega }(E)`$, the calculations can be prohibitively long. The energy ranges employed here correspond to $`E_{\text{max}}E_{\text{min}}=200400`$, or $`(E/N)_{\text{max}}(E/N)_{\text{min}}=1.02.0`$ per particle. In order to generate $`\mathrm{\Omega }(E)`$ over a wide energy range, neighboring energy windows were constructed in such a way as to overlap by half their width; every region of the energy axis is covered by at least two independent simulations. Two sets of independent simulations with energy windows of 200 or 400 units wide were pursued here. At low energies, convergence was only possible with relatively narrow energy ranges. The lowest energy range employed here starts at $`E=1860`$ ($`E/N=7.44`$), which is slightly above the range of estimated intrinsic energies for this system Sciortino et al. (1999) and well above the crystal energies Middleton et al. (2001). After the density of states converged to within a certain accuracy the update of $`\mathrm{\Omega }(E)`$ was stopped. A simple multicanonical run was then employed, where the inverse of the DOS was used as the weighting function. After a long run (25 million steps) the density of states was corrected by adding the logarithm of the final histogram. This was necessary as we did not run the simulations to update factors of $`f=\mathrm{exp}(10^8)`$ as in the original work by Wang and Landau Wang and Landau (2001a, b). It was verified during the course of the simulation that the logarithmic energy histograms kept increasing homogeneously (Figure 1a). This was monitored by observing that the logarithms of the histograms taken at different points in time during the final run are filled homogeneously. In a theoretically ideal situation, the different (logarithmic) histograms would be parallel to each other. However, the stochastic nature of the simulation leads to departures from parallel behavior. Nonetheless, we see that the system never gets “trapped” in a certain energy region without visiting the others. This ensures that the system readily moves back and forth between high and low energy regions of phase space, which is analogous to moving through temperatures in other types of simulations. In the case of continuous degrees of freedom this additional monitoring is important.

In order to provide an assessment of the correlation time of the random-walk method, the mean squared displacement of the particles was also measured. It is observed that in the final multicanonical run approximately 2 million cycles are necessary for the mean-squared displacement to be comparable to the box length (in one MC cycle each particle is moved once). The simulations presented here are at least 10 times that length (cf. Figure 2). A positional autocorrelation function can be defined as

$$C_{pos}(t)=(\stackrel{}{x}_i(t_0+t)\stackrel{}{x}_j(t_0+t))(\stackrel{}{x}_i(t_0)\stackrel{}{x}_j(t_0))_{i,j,t_0}.$$
(4)
This measure of relaxation is more stringent than the mean-square displacement as it eliminates the possibility that blocks of particles might be moving together without too much mutual rearrangement. This function decays to zero within the simulation lengths considered here. Figure 2 also compares this function to that obtained from molecular dynamics runs at $`T=0.45`$. To the best of our knowledge, the longest runs reported in the literature Yamamoto and Kob (2000) for the system considered here have lasted $`10^5\tau `$ ($`\tau `$ is the usual dimensionless Lennard-Jones time, which typically corresponds to 100 timesteps). In $`10^5\tau `$, this function decays to about 70% of its initial value. The scale employed in Figure 2 assumes that the computational requirements for one MD timestep are comparable to those of one MC cycle.

In regions of overlap, the density of states corresponding to each window can only differ by a constant, which depends on the (arbitrary) number of histogram entries. The density of states over the entire energy range of interest is constructed by shifting local estimates of $`\mathrm{ln}\mathrm{\Omega }(E)`$ (corresponding to individual windows) until they coincide, in the middle of the overlap region.

The global density of states is therefore known to within a constant. Since internal energies are known exactly, the excess free energy of the system can be calculated as a function of $`T`$ according to

$$\beta F(T)^{\text{ex}}=\mathrm{ln}\underset{E}{}\mathrm{\Omega }(E)e^{\beta E},$$
(5)
where the brackets denote an ensemble average and where $`\beta =1/k_BT`$. Similarly, the average internal energy of the system is given by

$$\beta E(T)^{\mathrm{ex}}=\frac{_EE\mathrm{\Omega }(E)e^{\beta E}}{_E\mathrm{\Omega }(E)e^{\beta E}},$$
(6)
and the entropy can simply be determined from $`S=(UF)/T`$. Note that the total entropy also comprises an ideal-mixing contribution of the form $`1\mathrm{ln}\rho x_A\mathrm{ln}x_Ax_B\mathrm{ln}x_B`$, where $`x_A`$ represents the mole fraction of species $`A`$ in the mixture. There is a second way to access properties, directly from the microcanonical ensemble. For example, the internal energy as a function of temperature can be derived by differentiation of the entropy (i.e. $`k_B\mathrm{ln}\mathrm{\Omega }(E)`$) with respect to energy and and subsequent inversion of that curve. This relies on the microcanonical definition of temperature, namely $`T=\frac{dS}{dE}^1`$.

## III Results

Figure 3a shows the average internal energy of the binary Lennard-Jones glass former. Results by Yamamoto et al. Yamamoto and Kob (2000) obtained by replica exchange Molecular Dynamics are also shown in that figure. The agreement between the two sets of data is quantitative. At temperatures below $`T_{MCT}`$ we have performed additional, extensive simulations using biased-sampling ideas and parallel-tempering techniques. More specifically, the algorithms developed for these additional calculations use two-dimensional parallel tempering in temperature Marinari and Parisi (1992); Tesi et al. (1996); Hansmann (1997); Wu and Deem (1999); Yan and de Pablo (1999) and Hamiltonian Bunker and Dünweg (2001), and identity swap moves Grigera and Parisi (2001) augmented with configurational bias sampling de Pablo et al. (1992). These techniques permit simulations at temperatures below $`T_{MCT}`$, but become increasingly sluggish as temperature is decreased. Still, the agreement for the energy generated by those simulations and the Density of States technique is also good, thereby providing further consistency tests for the results presented here.

The entropy of the binary Lennard-Jones glass former is shown in Fig. 3b as a function of temperature. The points in the figure represent literature data generated by thermodynamic integration Coluzzi et al. (2000b); the entropies simulated in this work have been shifted by a constant to make them coincide with those data in the range $`0.5<T<0.8`$. The agreement between our results and those of Coluzzi et al. is excellent. In order to arrive at an estimate of the Kauzmann temperature, these authors extrapolate the liquid phase entropy below $`T=0.5`$ using an expression of the form $`S=aT^{0.4}+b`$. We find that such a functional form is in good agreement with our simulations; at lower temperatures, however, minor but systematic departures from our results are observed. This would be indicative of a slightly higher Kauzmann temperature than that reported in the literature ($`T_K0.3`$Sciortino et al. (1999); Coluzzi et al. (2000b), as the simulated entropy decays more rapidly than that anticipated by extrapolation.

## IV Discussion

The density of states as a function of temperature does not show any unexpected behavior over the entire energy region considered in this work. Its logarithm simply becomes steeper with decreasing energy (see Figure 1b), reflecting the fact that the number of accessible states becomes smaller. The system could conceivably undergo a gas-liquid (or gas-glass) phase transition at very low temperatures. To address this point we have also determined the pressure. Our results suggest that such a transition can occur at $`T0.2`$, where the pressure of the glass becomes equal to that of the gas $`p0`$. For the same system (but without cutoff corrections), Coluzzi and Parisi estimated such a transition at $`T0.3`$. The occurrence of a demixing transition can be ruled out by the shape of the various pair distribution functions (not shown here), which remain qualitatively unchanged in the range $`0.2<T<0.6`$. A very similar system has been found to crystallize at such low temperatures Middleton et al. (2001). However, crystallization is avoided in our calculations by restricting the simulations to energy ranges above those of the crystal. Moreover, crystallization has only been found in constant pressure simulations. The volume is kept constant in this work, thereby leading to frustration of crystallization.

The results presented in this work suggest that a Monte Carlo technique based on the concept of entropic sampling is capable of generating high-accuracy estimates of the equilibrium density of states of a binary Lennard-Jones glass former, down to temperatures below the mode coupling temperature, a region of temperature that previous studies of glass-forming liquids have avoided. With further refinement of the algorithm discussed in this work Yan and de Pablo (2003), we expect that reliable simulations in the near vicinity of the reported Kauzmann temperature will become possible.

## Acknowledgments

RF wants to thank the Emmy-Noether Program of the Deutsche Forschungs-Gemeinschaft for financial support.


format: text

text:
# Pathways to non-sequential multiple ionization in strong laser fields

## I Introduction

Laser pulses with peak intensities of the order of $`10^{14}`$ W/cm2 and wavelengths in the infrared can ionize several electrons simultaneously luk83 ; huillier83 ; boyer84 . Growing experimental evidence walker94 ; moshammer00 ; weber00prl ; weber00jpb ; weber00n ; wechenbrock01 ; feuerstein01 ; moshammer02  supports a rescattering process corkum93 ; kulander93  as the most likely mechanism for the rapid transfer of a sufficient amount of energy to the electrons. The interaction between the rescattered electron and the others takes place close to the nucleus, leading to the formation of a short lived, highly excited compound state with electrons near to the nucleus. This compound state will decay by emitting a single electron (the most likely decay mode) or several electrons. The highly excited electrons move so quickly that the changes in the field can be neglected: the electrons see an essentially static electric field. Experimental studies of the final momentum distribution for double ionization show that the electron momenta parallel to the field are strongly correlated weber00n ; wechenbrock01 ; feuerstein01 ; moshammer02 . Building on the picture of an intermediate compound state we have proposed a classical origin for the electron correlations eckhardt01 ; eckhardt01pra1 ; eckhardt01silap ; eckhardt01pra2 ; eckhardt01epl . If the energy of the compound state is above the zero field threshold energy no constraints on the outgoing dynamics are imposed. If the energy is below double ionization can occur only if the field is still present. Then the electrons will aim for the Stark saddle, but the electron repulsion will prevent them from crossing the saddle in the same place at the same time. Furthermore, if one electron is ahead it will push the others back, thus preventing multiple ionization. Therefore, the electrons have to escape in a configuration that balances the mutual repulsion. For two electrons this means that the electrons move side by side with respect to the field axis. Within this symmetric subspace it is easy to identify a saddle which the electrons have to cross on their way to double ionization: it is a stationary configuration of two electrons side by side close to the Stark saddle for a single electron (assuming, as mentioned before, that the field is stationary). The number of saddle configurations increases with the number of electrons: there are two configurations for three electrons, with different energies and different threshold exponents eckhardt01pra2 .

We here extend our earlier study of such saddle configurations for two and three electrons to four or more. The aim is to identify the possible saddles, their critical energies and their threshold exponents. The number of configurations increases rapidly with the number of electrons. Beginning with five electrons also several nearby almost degenerate saddles appear. Most interestingly, for four and fewer electrons all electrons are equivalent in the dominant mode, but this is no longer true for six or more.

We conclude this introduction with notation and the formulation of the Hamiltonian. Since the electrons escape from the highly excited compound state so quickly, we may assume that the field is constant during the escape. Let the electric field point in the $`z`$-direction and let the electrons be labeled $`i=1,\mathrm{},N`$ with positions $`𝐫_i=(x_i,y_i,z_i)`$. Then the Hamiltonian is, in atomic units,

$$H_N=\underset{i=1}{\overset{N}{}}\frac{𝐩_i^2}{2}\underset{i=1}{\overset{N}{}}\frac{Z}{|𝐫_i|}+\underset{i<j}{\overset{N}{}}\frac{1}{|𝐫_i𝐫_j|}F\underset{i=1}{\overset{N}{}}z_i.$$
(1)
Simple scaling shows that the field strength $`F`$ can be absorbed, so that without loss of generality we can take $`F=1`$ eckhardt01pra1 ; eckhardt01pra2 .

The aim now is to determine the stationary points of this Hamiltonian, where the derivatives with respect to momenta and positions vanish. For the most part this is not possible analytically. The numerical method of choice is the Newton-Raphson iteration numericalrecipes . The required calculation of the matrix of second derivatives can be turned to an advantage since it allows for an immediate determination of the stability properties of configuration. The main disadvantage of the method, but also of all other numerical methods, is that they do not provide any measure by which to judge whether all solutions have been found. On the practical side we can make sure that we have hit at least all the configurations with large domains of attraction in the Newton-Raphson method by starting with a sufficiently large number of randomly selected initial conditions. We used up to $`10^6`$ initial conditions. The fact that fairly similar but definitely different configurations could be detected may be interpreted as strong evidence that all configurations have indeed been found.

The stability analysis of the stationary point in the full phase space allows us to determine the behavior of the cross-section for non-sequential escape close to the threshold eckhardt01epl . Among all Lyapunov exponents of the saddle there is one, $`\lambda _r`$, whose eigenvector components point in the same direction along the field axis: it corresponds to a simultaneous motion of all electrons in the same direction away from the saddle. Borrowing terminology from chemical reactions, we call this subspace the reaction coordinate. Because of the repulsion between electrons all configurations have additional unstable eigenspaces, which enter in the threshold exponent. If the initial energy of the system equals the saddle energy only a trajectory living in the symmetric subspace can lead to a simultaneous escape of all electrons. This reduces the dimensionality of the problem and the cross-section vanishes. For higher energy some deviations from the symmetric motion are possible, giving a finite volume of initial conditions and a non-vanishing cross-section. The dependence of the cross-section on energy $`\sigma (E)`$ close to the saddle energy $`E_N`$ can be obtained in the spirit of the Wannier analysis wannier53 ; rau84 ; rost98 ; rost01phe , resulting in

$$\sigma (E)(EE_N)^\mu ,$$
(2)
with an exponent

$$\mu =\frac{1}{\lambda _r}\underset{i=1}{\overset{n_u}{}}\lambda _i.$$
(3)
The $`\lambda _i`$’s are the positive Lyapunov exponents of the saddle except for the reaction coordinate exponent $`\lambda _r`$. The cross section is large if the exponent is small, i.e. if the saddle is crossed quickly (large $`\lambda _r`$) or if the differences from the symmetric motion grow slowly (small $`\lambda _i`$). These cross section exponents are an additional characteristic of the multiple ionization process.

The remainder of the paper is organized as follows: in the next section we discuss a few obvious, and highly symmetric configurations. This is followed by our numerical results for up to eight electrons in section III. We conclude with a few remarks in section IV.

## II Symmetric configurations

### II.1 All electrons on a ring

We assume that all electrons are situated in a plane perpendicular to the field and that they obey a $`C_{Nv}`$ symmetry. The reflection symmetry limits the momenta to be parallel to the symmetry planes and thus confines the motion to a dynamically allowed subspace in the high-dimensional $`N`$-body phase space. That is, if in the full phase space of the $`N`$-body problem initial conditions are prepared in this subspace they will never leave it. In the symmetry subspace we can use cylindrical coordinates for the electrons, $`z_i=z`$, $`\rho _i=\rho `$ and $`\phi _i=2\pi i/N`$. Then, for zero total angular momentum along the field axis and in the scaled variables where $`F=1`$, the Hamiltonian (1) can be reduced to

$$H(p_\rho ,p_z,\rho ,z)=\frac{p_\rho ^2+p_z^2}{2N}\frac{N^2}{\sqrt{\rho ^2+z^2}}+\frac{W}{2\rho }Nz,$$
(4)
with

$$W=\underset{i<j}{\overset{N}{}}\frac{1}{\mathrm{sin}\left[\frac{\pi }{N}(ji)\right]}=\underset{k=1}{\overset{N1}{}}\frac{Nk}{\mathrm{sin}\left[\frac{\pi k}{N}\right]}.$$
(5)
The Hamiltonian (4) possesses a stationary point at $`p_\rho =0`$, $`p_z=0`$, $`\rho =\rho _N`$ and $`z=z_N`$ where

$`\rho _N`$ $`=`$ $`\sqrt{{\displaystyle \frac{W}{2N}}}\left[\left({\displaystyle \frac{2N^2}{W}}\right)^{2/3}1\right]^{1/4}`$ (6)
$`z_N`$ $`=`$ $`\sqrt{{\displaystyle \frac{W}{2N}}}\left[\left({\displaystyle \frac{2N^2}{W}}\right)^{2/3}1\right]^{3/4}.`$ (7)
and

$`E_N`$ $`=`$ $`\left[2N^2\left({\displaystyle \frac{2}{NW}}\right)^{1/6}\sqrt{2NW}\right]`$ (9)
$`\left[\left({\displaystyle \frac{2N^2}{W}}\right)^{2/3}1\right]^{1/4}.`$
The first few saddles with electrons on a ring are shown in Fig. 1.

For $`kN`$ the repulsion energy $`W`$ behaves like $`W=(N^2/\pi )(1/k)`$. This suggests that $`W`$ increases like $`N^2\mathrm{ln}N`$, and indeed numerical evidence supports

$$W(0.3N^2+0.3N3.1)\mathrm{ln}N.$$
(10)
Since $`W`$ increases faster than quadratically, the configurations exist for finite $`N`$ only. Specifically, we find that the energy of the configuration first decreases but then increases again, and that there is no configuration for $`N=473`$ or more particles 111This discussion corrects the one given in eckhardt01 , where a much smaller number of possible configurations was obtained using $`W=N(N1)/\mathrm{sin}(\pi /N)`$, with the large $`N`$ behavior $`WN^3`$. This overestimate of the repulsion between the electrons is responsible for the smaller number of possible configurations.. The properties of the configurations, including stability and cross section exponents, are listed in table 1.

### II.2 One electron in the center

The simplest extension from the symmetric ring is to place one electron on the field axis and the others on a ring as shown in Fig. 2. Such configurations have three characteristic lengths, the radius of the ring $`\rho _{N1}`$, the position of the ring $`z_{N1}`$ and the position of the electron on the axis $`z_c`$. No analytical solutions have been found. The numerically determined characteristics are collected in table 2.

### II.3 All electrons on a line

A third class of configurations has all electrons in a plane through the field axis. The electrons are placed like beads on a string. With the confinement of electrons to an almost linear arrangement the saddle energies are higher than in the previous cases and the cross section exponents are also much larger. The properties of the saddles are listed in table 3, and their positions are shown in Fig. 3.

## III All configurations

A summary of the symmetric configurations discussed in the previous section is provided in Fig. 4, which shows the saddle energies for different $`N`$. For the range of $`N`$ included in the figure the saddle energies for all configurations decrease, but as the example of the ring configurations shows, they can be expected to increase for sufficiently large $`N`$. There is a crossover in the configuration that gives the minimum, from the ones with all particles on a ring for $`N<5`$ to the ones with a one electron in the center for $`N>5`$. The case $`N=5`$ is marginal (see below). In order to see what other configurations one has to consider, we turn to a discussion of configurations found with the Newton-Raphson method for electron numbers between four and eight. The cases of two and three electrons have been analyzed previously eckhardt01 ; eckhardt01pra1 ; eckhardt01silap ; eckhardt01pra2 ; eckhardt01epl .

### III.1 Four electrons

For four electrons there are 4 configurations, shown in Fig. 5. Their properties are listed in table 4. The dominant non-sequential ionization takes place in the vicinity of the saddle corresponding to all electrons on a ring in a plane perpendicular to the field axis, as in the cases for two and three electrons. This saddle possesses both the lowest energy and the smallest critical exponent and thus the most favorable dependence of the cross-section on energy.

The next two configurations are almost degenerate in energy and have similar critical exponents. They form a triangle with an additional electron placed in the center. The saddle $`\nu =2`$ possesses a $`C_{3v}`$ symmetry, but the state $`\nu =3`$ has a $`C_v`$ symmetry only. The last $`\nu =4`$ configuration, with all electrons in a plane, has significantly higher energy and larger critical exponent compared to the previous states.

### III.2 Five electrons

For five electrons we find 5 configurations shown in Fig. 6, with properties given in table 5. The five electron problem constitutes actually a marginal case where the configuration of all electrons on a ring starts loosing its dominance. The energy of this configuration is only slightly lower than the energy of the state $`\nu =2`$ with one electron in the center and the others on a ring. The state $`\nu =2`$ possesses, however, a smaller critical exponent than the $`\nu =1`$ configuration. Therefore, unless there is a wide disparity in prefactors in the cross section law or residual correlations in the compound state that prefer one state over the other, we expect that both will contribute significantly to the process of non-sequential five electron escape.

The configurations $`\nu =3`$ and $`\nu =4`$ have the $`C_v`$ symmetry and reveal higher energy and larger critical exponents than the previous ones. The last $`\nu =5`$ configuration with all electrons on a line is the least significant because of its high energy and its high critical exponent.

### III.3 Six electrons

For six electrons we find 11 saddles, with properties listed in table 6. The first five configurations are shown in Fig. 7. The state with all particles on a ring is no longer dominant, it appears in the middle of table 6. The saddle with an electron in the center and the others on a ring dominates the non-sequential ionization process. It has both the lowest energy and the smallest critical exponent. The states $`\nu =2`$ and $`\nu =3`$ look similar, but differ in their symmetries: state $`\nu =2`$ has a $`C_v`$ symmetry, state $`\nu =3`$ has a $`C_{3v}`$ symmetry.

### III.4 Seven electrons

For seven electrons we have 14 states, with the properties listed in table 7. In Fig. 8 the first five states are shown. The configuration with one electron in the center and the others on a ring dominates the simultaneous electron escape, just as for the six electron problem.

### III.5 Eight electrons

For eight electrons we have 26 configurations, listed in table 8. The dominant configuration has two electrons in the center and the other six around it, see Fig. 9. The $`\nu =2`$ saddle also has two electrons in the center, but different arrangements of the outer electrons.

The state with one electron in the center and the others on a symmetric ring is the fifth configuration, the one with all electrons on a ring the 17th configuration. The saddles $`\nu =3`$ and $`\nu =4`$ are similar to the $`\nu =5`$ one but they only have a $`C_v`$ symmetry, not the $`C_{7v}`$ one as for state $`\nu =5`$. All of these configurations possess similar characteristics, see table 8.

The state with all particles on a line has the highest energy and the largest critical exponent. This is the case not only for the eight electron problem but for all $`N`$ values analyzed here.

## IV Concluding remarks

The analysis presented here confirms a trend anticipated from our analysis of triple ionization: that with an increasing number of electrons more and more saddle configurations and thus channels for non-sequential multiple ionization have to be considered. But the analysis also shows that, except for the somewhat marginal case of fivefold ionization, there is a well separated lowest lying saddle that dominates non-sequential multiple ionization. Table 9 collects the symmetries and the threshold exponents of this lowest lying state.

A few experiments with multiple ionization of atoms have already been performed. However, momenta of the ionization products have been measured for double and triple ionizations only. Currently, measurements of momenta of more than two electrons are not available, not even of a single component, say the one along the field axis. But already the measurement of parallel momenta for just two electrons would suffice for a preliminary test of the ionization pathways considered here. In particular, the breaking of symmetry between the outgoing electrons should be accessible. As long as the electrons are interchangeable, the distribution of parallel electron momenta should be like that for double ionization, with a clear preference for similar momenta. For $`N=5`$ and more escaping electrons additional structures in the momentum distribution should appear, since the lowest lying saddle has non-equivalent electrons, see Figs. 5 to 9. The momentum distributions for the two $`N=3`$ saddles are discussed in eckhardt01pra2 , and similar modifications in shape can be expected for more electrons.

The analysis presented here assumes simultaneous escape of all electrons in a multiple ionization event. It is of course also possible that first one or a few electrons escape, then the next batch of electrons leaves and so on, before one arrives at the final multiply charge ion. Sequential ionization becomes significant when the field intensity is high enough to ionize electrons from a highly excited state of the ion left behind when the first electrons are gone. Nevertheless, close to the threshold for non-sequential multiple ionization the characteristics of the non-sequential process as described here should be identifiable.

## V acknowledgment

We thank J. M. Rost for a remark that triggered this investigation. This work was partially supported by the Alexander von Humboldt Foundation, the Deutsche Forschungsgemeinschaft and the KBN through grant 5 P03B 088 21.


text:
# Natural Coronagraphic Observations of the Eclipsing T Tauri System KH 15D: Evidence for Accretion and Bipolar Outflow in a WTTS1footnote 11footnote 1Based on Ultraviolet-Visual Echelle Spectrograph observations collected at the European Southern Observatory’s Very Large Telescope via Director’s Discretionary Time, within the observing program P267.C-5736.

## 1 Introduction

KH 15D is a unique system in which a pre-main sequence star is periodically occulted by extended, non-luminous matter, presumably part of a circumstellar disk (Herbst et al., 2002). Every 48.36 days the star fades over 2-3 days by $``$3.5 mag and remains near minimum light for $``$20 days. The eclipse duration has been increasing with time, by $``$1-2 days/yr. Hamilton et al. (2001) obtained low resolution spectra in and out of eclipse and concluded that the star was a K7 weak-line T Tauri star (WTTS). The WTTS classification, based on the equivalent width (EW) of the H$`\alpha `$ line out of eclipse, is supported by the lack of any IR111H-K = 0.14, K-L = 0.03 out-of-eclipse, K. Haisch, private communication. or UV excess emission and the absence of substantial photometric variability outside of eclipse. Based on its membership in NGC 2264 (Sung et al., 1997), the estimated mass of the K7 star is 0.5-1 M and its age is 2-4 Myr.

In an attempt to learn more about this unique system we obtained high resolution spectra during the eclipse of December 2001. We hoped to determine whether KH 15D is a radial velocity variable, whether there was evidence for additional light in the system beyond that of the K7 star, and whether the occultation had any effect on the spectrum. In fact, we found dramatic changes in the line profiles of H$`\alpha `$ and H$`\beta `$, as well as weak forbidden emission lines that become much more visible during eclipse. Evidently, the KH 15D system behaves like a “natural coronagraph”, allowing us to see details of its close circumstellar environment during eclipse. Our evidence suggests that this WTTS is actively accreting gas and driving a bipolar outflow, although probably not at the rate of a typical CTTS. This calls into question the common practice of associating WTTS characteristics with the absence of an accretion disk.

## 2 Observations and the Absorption Spectrum

High resolution echelle spectra of KH 15D were obtained on the nights of UT 2001 Nov. 29, when it was in its bright state just prior to eclipse, on UT 2001 Dec. 14, just past mid-eclipse, and again on UT 2001 Dec. 20 during egress (see Figure 1). These data were collected with the Ultraviolet-Visual Echelle Spectrograph (UVES) on the European Southern Observatory’s Very Large Telescope (VLT), at Mount Paranal, Chile. The wavelength range is $``$4780 to 6810 Å. A 50 Å gap centered on 5800 Å is present due to use of the red arm which employs a mosaic of two 4096 x 2048 CCDs (D’Odorico et al., 2000). With a 1″ slit, the spectral resolution is $``$44,000. The spectra presented here were reduced via the UVES pipeline. Information regarding the reduction process can be found at http://www.eso.org/instruments/uves/. As a check on this, the spectra were also reduced in the manner described by Valenti (1994). Both procedures make use of a sky subtraction algorithm. No significant differences were found between the data reduced by these techniques and the UVES pipeline reduction is adopted here.

The UVES spectra confirm that, out of eclipse, the primary light source in the KH 15D system is a K7 WTTS. The 29 Nov. spectrum was visually compared in detail to 61 Cyg B, a typical K7 V, to look for any evidence of a second source of light in the system. No additional features were found. This spectrum was also compared to a rotationally broadened, synthetic K7 V template to determine its $`vsin(i)`$, which we estimate to be $`<`$ 5 km/s. The EW of the LiI 6707 feature is 0.401 Å $`\pm `$ 0.001, consistent with the value of 0.47 Å $`\pm `$ 0.05 measured on a low resolution spectrum by Hamilton et al. (2001).

The out-of-eclipse spectrum was cross-correlated against HD 55999, a standard observed on the same nights, to determine the radial velocity of KH 15D. This was done with the IRAF222Image Reduction and Analysis Facility, written and supported by the IRAF programming group at the National Optical Astronomy Observatories (NOAO) in Tucson, Arizona. NOAO is operated by the Association of Universities for Research in Astronomy (AURA), Inc. under cooperative agreement with the National Science Foundation. task FXCOR. Since the spectra were obtained on two CCDs covering different wavelength regions, a cross-correlation was performed on each, avoiding the emission line regions of H$`\alpha `$, and H$`\beta `$. A heliocentric radial velocity of +9.0 km/s $`\pm `$ 0.2, was found by this procedure for KH 15D on 29 Nov 2001.

The egress spectrum was cross-correlated against the out-of-eclipse spectrum in the same manner to look for any radial velocity variation between 29 Nov. and 20 Dec. A difference in heliocentric radial velocity of 3.3 $`\pm `$ 0.6 km/s was found. The radial velocity of the star on 20 Dec. was, therefore, +12.3 km/s $`\pm `$ 0.6. Whether this detection of radial velocity variation means that the K7 star is a spectroscopic binary remains to be seen. Photospheric line profile variations are expected in a partially eclipsed star and the likely importance of scattered radiation near minimum light further complicates the interpretation. Clearly a more extensive radial velocity study of the system is needed and is underway (J. Johnson & G. Marcy, private communication).

## 3 The Emission Line Spectrum

In Figures 2-4, we show the H$`\alpha `$, H$`\beta `$, and \[OI\] 6300 emission-line profiles for KH 15D. Each spectrum has been flux calibrated to a relative scale using an R magnitude for the date, as given in Table 1. The quoted uncertainties on the magnitudes reflect a small degree of non-simultaneity in the spectral and photometric data as well as the need to transform from I (where the data were more numerous) to R. Since spectral and photometric data were taken within an hour of each other, except at maximum light (when the brightness is nearly constant with time) and since the color variation is quite small ($``$0.1 mag in R-I over the full brightness range), the uncertainty in R is only $``$0.1 mag. Julian dates for each observation, measured I magnitudes, derived R and V magnitudes, measured EWs, and derived H$`\alpha `$ and \[OI\] fluxes are listed in Table 1. The EW measurements refer to the entire profile, both red- and blue-shifted components combined. On Figs. 2 and 3, the arrows indicate where HI nebular emission could have affected the profile, whereas on Fig. 4 the arrow indicates where \[OI\] $`\lambda `$6300 from the night sky may not have been removed completely.

Figure 2 shows that dramatic changes occurred in the H$`\alpha `$ emission line profile of KH 15D during eclipse and egress. Out of eclipse, the star appears to be a WTTS with an EW(H$`\alpha `$) $``$2 Å, but with a double-peaked profile, having a central absorption and faint, but clearly detectable broad wings. A comparison with the H$`\alpha `$ profiles of other WTTSs (Hartmann, 1982; Mundt et al., 1983; Finkenzeller & Basri, 1987; Edwards et al., 1994; Reipurth et al., 1996) shows that only 3 out of 19 such stars exhibited double-peaked profiles similar to KH 15D, with UX Tau A being the most similar (see Reipurth et al. (1996)). Most WTTSs show narrow, single-peaked emission lines. It is also interesting that the central absorption feature appears to extend below the stellar continuum. This is unusual for any TTS, weak or classical. During eclipse and egress, the “natural coronagraph” effect is clear. Near mid-eclipse, the EW of H$`\alpha `$ grows to $``$40 Å, while the relative flux drops by $``$50%. During egress the EW rises to $``$70 Å, as the flux increases, exceeding even its out-of-eclipse value. In all cases, a double-peaked profile is observed. The emission line profile during mid-eclipse extends to $`\pm `$ 200 km/s with significantly less flux in the extended wings, as compared to the profile during egress, which extends to $`\pm `$ 300 km/s.

Figure 3 shows the corresponding H$`\beta `$ line profiles. Again, dramatic EW and line profile changes are evident. Although the out-of-eclipse H$`\beta `$ line profile is heavily affected by the underlying stellar absorption spectrum, it is similar to H$`\alpha `$. In the egress spectrum, however, the emission occurs primarily on the blue side with a wing extending to -300 km/s. The asymmetry in this line is much more striking than H$`\alpha `$. During mid-eclipse, although little H$`\beta `$ flux is present, the shape of the line profile appears similar to egress. The weak features visible in this profile at - 25 km/s are probably due to an improper background subtraction of the H$`\beta `$ emission from NGC 2264.

Figure 4 shows the emission-line profiles for the \[OI\] 6300 Å line. The \[OI\] line, which is weak, but clearly discernable at maximum light, becomes prominent at mid-eclipse and during egress. The line flux seems to be about the same out of eclipse and during egress, although slightly higher at mid-eclipse. We caution against any extreme interpretation of this measurement. The profile during mid-eclipse was disturbed by an improper background subtraction, and we feel that our errors are most likely underestimated. However, this behavior would indicate that none of the \[OI\] emitting zone suffers variable occultation at the phases of our observations. This is fully expected given the spatial extent of the (bipolar) forbidden line emitting regions in CTTSs. The peak of the high velocity component of the \[OI\] emission in typical CTTSs originates at about 30 AU from the star (Hirth et al., 1997).

The EW of the \[OI\] $`\lambda `$6300 line in the out-of-eclipse spectrum is about 0.17 Å, which is quite large for a WTTS. Only one of the ten WTTSs in Table 3 of Hartigan et al. (1995) has a detectable \[OI\] $`\lambda `$6300 line (with an EW = 0.5 Å). The others have upper limits of about 0.06 Å. Most CTTSs in that table, on the other hand, have EWs of 0.5 - 3 Å. The profiles obtained during mid-eclipse and egress suggest that we have two emission peaks, at about -20 and +18 km/s with emission wings extending -60 to +50 km/s, respectively. These profiles are quite different from the \[OI\] $`\lambda `$6300 profile seen in most strong-emission CTTSs (with strong veiling), which often have a high-velocity component at -100 to -150 km/s (resulting from the jet) and a low-velocity component at about -20 km/s (Hartigan et al., 1995; Hirth et al., 1997). In addition, the profile is quite different from that of a CTTS with small veiling and IR excess (see the bottom profile in Figure 11 of Hartigan et al. (1995)). These latter profiles are usually single-peaked and unshifted in velocity. Although the profile shape for KH 15D is not fully clear due to imperfect subtraction of the emission of \[OI\] at 6300 Å in the night sky, it is most likely double-peaked. Such a profile can be most easily explained by a bipolar jet moving nearly perpendicular to the line of sight, quite consistent with our expectation that the disk associated with KH 15D is viewed nearly edge on, resulting in a very small radial velocity separation between the two jets. The profile catalog of Hartigan et al. (1995) contains several examples of bipolar jet sources (RW Aur, AS 353A, DD Tau), where the two jet components are clearly separated due to favorable inclination angles.

## 4 Discussion

It is important to keep in mind that we do not yet know the geometry of the KH 15D system. In particular, it is not certain whether the occultation proceeds along a line perpendicular to or parallel to the rotation axis of the K7 star (or neither) or to what extent the rotation axis, the presumed stellar magnetic axis and the orbital plane are aligned. It is not even known whether, relative to the system’s center of mass, it is the K7 star or the occulting matter which is primarily in motion. Detailed modeling of the system must obviously await clarification of these basic issues. However, a preliminary qualitative interpretation of the emission line variations is possible and provided in this section.

We believe that the observed spectral variations of KH 15D can be understood qualitatively in terms of a weakly accreting, “scaled down” classical T Tauri star (see e.g., Muzerolle et al. (2001)) whose photosphere and magnetosphere are periodically occulted by a relatively sharp knife-edge (as described in Herbst et al. (2002), their Fig. 5). We expect such a star to have a bipolar jet region, which we assume, by analogy with the CTTS, is revealed by the forbidden emission line radiation. Assuming a jet velocity of 200 km/s (which is the average jet velocity for CTTSs derived by Hirth et al. (1994)) and adopting a radial velocity of $`\pm `$ 20 km/s for the \[OI\] $`\lambda `$6300 peaks, we derive an inclination angle for the jet to the line of sight of 84$`\mathrm{°}`$. This is consistent with a general picture for the system in which we view the K7 star close to the orbital plane of its circumstellar (or circumbinary) disk and the jets emerge roughly perpendicular to the disk plane, as for the star associated with HH 30 and other examples imaged by the Hubble Space Telescope (Ray et al., 1996). The absence of any significant variation in the profiles or flux of the forbidden line radiation during eclipse is consistent with the expectation, based on the CTTS analogy, that it arises at distances of tens of AU’s from the star, beyond the region variably occulted.

The behavior of the hydrogen lines is complex because some components do arise close to the star and, therefore, suffer variable occultation effects along with the photosphere. Since we expect that the H$`\alpha `$ line has a much higher optical depth than the H$`\beta `$ line, it is obvious that any emission line region will be more extended in H$`\alpha `$ than H$`\beta `$. If we assume that the magnetic axis of the K7 star is tilted toward us at $``$ 5-10$`\mathrm{°}`$, which is a reasonable assumption supported by the \[OI\] line profiles, at mid-eclipse, one can qualitatively understand the H$`\alpha `$ line profile as resulting from low velocity material in the outer, more extended H$`\alpha `$ emission region while the star and the H$`\beta `$ line forming region are obscured by the occulting disk material. This would also explain why there is almost no H$`\beta `$ flux during mid-eclipse. During egress, we expect most of the H$`\alpha `$ emission line region, as well as some of the H$`\beta `$ emission line region, which is closer to the star, to be visible. The H$`\alpha `$ line profile during egress has a low velocity emission peak with a central absorption feature similar to what is seen in the out-of-eclipse profile. Additionally, two “shoulders” appear along the profile at about $`\pm `$ 150 km/s extending out to $`\pm `$ 300 km/s. These “shoulders” could be due to material rotating in the outer parts of the magnetosphere. However, given that the $`vsin(i)`$ is measured to be $`<`$ 5 km/s, and that the magnetosphere only extends out to about 5-6 stellar radii, the maximum rotational velocity is about 5-6 $`vsin(i)`$ or 25-30 km/s, making it difficult to explain the H$`\alpha `$ emission-line profile with a rotating magnetopshere.

A more attractive hypothesis is that the high velocity “shoulders” on the H$`\alpha `$ line, so prominent in the egress spectrum, arise from material falling along magnetic accretion columns. This interpretation can also qualitatively account for the H$`\beta `$ emission line profile during egress. Adopting values for the mass and radius of the K7 star from Table 1 of Hamilton et al. (2001), a free-fall velocity of about 380 km/s can be associated with material at the surface of the star. Since H$`\beta `$ is produced much closer to the star, the blue-shifted wing extending to nearly -350 km/s could be representative of material accreting along magnetic field lines near the pole. The asymmetry seen in the H$`\beta `$ emission line profile is most likely due to the fact that the star is slightly inclined toward our line of sight. The reappearance of both a blue and red wing to the H$`\alpha `$ line during the early part of egress, when most of the stellar photosphere is still occulted, shows that both red-shifted and blue-shifted gas is present along the same line of sight towards the small portion of the photosphere and magnetosphere that is being uncovered first. This also supports an accretion interpretation as opposed, say, to a rotation interpretation for the high velocities.

## 5 Conclusions

It appears that KH 15D, although clearly a WTTS by the usual criteria of weak H$`\alpha `$ emission, absence of UV or IR excess emission, and relative photometric stability, is still undergoing active accretion and driving a bipolar outflow. It provides a cautionary example against assuming that all stars with WTTS characteristics no longer have accretion disks. The unique geometry of this system, in which a relatively sharp-edged occulting mask crosses the photosphere has created a “natural cornagraph” which enhances the visibility of the star’s magnetosphere during eclipse. The occultation also evidently crosses the inner portion of the star’s magnetosphere, where high velocity gas motions arise, probably from magnetically channeled accretion. Synoptic studies of this star may ultimately allow us to reconstruct aspects of the structure of its magnetosphere with spatial resolution that will be unobtainable in other objects for decades to come.

We thank U. Bastian, M. Ibraghimov, J. Johnson, G. Marcy, F. Vrba, and J. Aufdenberg for helpful conversations and useful data on this star. W.H. acknowledges support from NASA through its Origins of Solar Systems program.


text:
# 1 Introduction

## 1 Introduction

A non-monotonic energy dependence of the $`K^+`$/$`\pi ^+`$ ratio with a sharp maximum close to 30 A$``$GeV is observed in central Pb+Pb collisions . Within a statistical model of the early stage , this is interpreted as a sign of the phase transition to a QGP, which causes a sharp change of the energy dependence of the strangeness to entropy ratio. This observation naturally motivates us to study the production of multistrange hyperons ($`\mathrm{\Xi }`$, $`\mathrm{\Omega }`$) as well as the function of the beam energy. Although this is the main motivation for our study of $`\mathrm{\Omega }`$ production in p+p and Pb+Pb collisions, there are two other points of interest. In fact, $`\mathrm{\Omega }^{}`$ and $`\overline{\mathrm{\Omega }}^+`$ production in p+p collisions could give quite detailed information on the dynamics of strangeness production in those collisions. From string-hadronic models, it is expected that the $`\overline{\mathrm{\Omega }}^+`$/$`\mathrm{\Omega }^{}`$ ratio in p+p interactions at SPS energies is higher than 1 whereas it is predicted to be smaller than 1 in a hadron gas model .
Furthermore it was suggested that the kinematic freeze-out of $`\mathrm{\Omega }`$ takes place directly at hadronization . If this is indeed the case, the transverse momentum spectra of the $`\mathrm{\Omega }`$ directly reflect the transverse expansion velocity of a hadronizing QGP .
In this report we show preliminary results on $`\mathrm{\Omega }^{}`$ and $`\overline{\mathrm{\Omega }}^+`$ production in central Pb+Pb collisions at 30, 40 and 158 A$``$GeV and compare them to data from p+p interactions at 158 GeV. We also show fully corrected rapidity and $`m_T`$ spectra of $`\mathrm{\Omega }^{}`$ and $`\overline{\mathrm{\Omega }}^+`$ in central Pb+Pb collisions at 158 A$``$GeV. Integration of the rapidity spectra gives the total yield of $`\mathrm{\Omega }^{}`$ and $`\overline{\mathrm{\Omega }}^+`$ in full phase space.

## 2 The NA49 experiment

The NA49 detector  (see Fig. 1) is a large acceptance hadron spectrometer at the CERN SPS, consisting of four TPCs. Two of them, the Vertex TPCs (VTPC), are inside a magnetic field for the determination of particle momenta. The ionisation energy loss (dE/dx) measurement in the two Main TPCs (MTPC), which are outside the magnetic field, is used for particle identification.

Central collisions were selected by a trigger using information from a downstream calorimeter (VCAL), which measures the energy of the projectile spectator nucleons.

## 3 Analysis and data sets

The $`\mathrm{\Omega }^{}`$ production is analyzed using its decay channel : $`\mathrm{\Omega }^{}`$ $``$ $`\mathrm{\Lambda }`$ $`+`$ $`K^{}`$ (BR=67.8 % ) and $`\mathrm{\Lambda }`$ $``$ $`\pi ^{}`$ $`+`$ $`p`$ (BR=63.9 % ). For $`\overline{\mathrm{\Omega }}^+`$ the charge conjugated channel is used. For the $`\mathrm{\Omega }^{}`$ ($`\overline{\mathrm{\Omega }}^+`$) analysis all $`\mathrm{\Lambda }`$ ($`\overline{\mathrm{\Lambda }}`$) candidates with a reconstructed invariant mass close to the nominal value ($`|`$$`\mathrm{\Delta }`$$`M_\mathrm{\Lambda }`$$`|`$ $`<`$ 0.005 GeV/$`c^2`$) are combined with the $`K^{}`$ ($`K^+`$) candidates. In order to identify the secondary vertex, both are extrapolated back to the target and the point of the closest approach is calculated . The resulting combinatorial background is reduced by applying various cuts. The contribution of false $`\mathrm{\Lambda }`$ ($`\overline{\mathrm{\Lambda }}`$) candidates can be reduced by selecting the decay (anti-)proton candidates using their energy loss in the TPCs. The same procedure is used to reduce background in the sample of $`K^{}`$ ($`K^+`$) candidates. The reconstructed $`\mathrm{\Omega }^{}`$ ($`\overline{\mathrm{\Omega }}^+`$) candidates have to point to the interaction vertex, while the $`K^{}`$ ($`K^+`$) candidate and the $`\mathrm{\Lambda }`$ ($`\overline{\mathrm{\Lambda }}`$) candidate should miss it. Fig. 2 shows the invariant mass distribution of $`\mathrm{\Omega }^{}`$ ($`\overline{\mathrm{\Omega }}^+`$) candidates in central Pb + Pb at 158 A$``$GeV (left), 40 A$``$GeV (middle) and p + p at 158 A$``$GeV (right). Only pairs in the kinematical regions specified in table 1 were used in figures.

A clear signal at the expected mass position ($`M_\mathrm{\Omega }`$ = 1.67245 $`GeV`$$`/`$$`c^2`$ ) is seen for all studied reactions. To estimate the number of $`\mathrm{\Omega }^{}`$ ($`\overline{\mathrm{\Omega }}^+`$) decays, the number of entries was integrated in the invariant mass interval \[$`M_\mathrm{\Omega }`$ $``$ $`\mathrm{\Delta }`$ $`M`$, $`M_\mathrm{\Omega }`$ $`+`$ $`\mathrm{\Delta }`$ $`M`$\] with $`\mathrm{\Delta }`$ $`M`$ = 0.010 $`GeV`$$`/`$$`c^2`$ and the background contribution was subtracted using a polynominal interpolation.

For the 158 A$``$GeV Pb+Pb events the raw signals were corrected for the branching ratios, the geometrical acceptance and reconstruction efficiency .

## 4 Results

### 4.1 Energy dependence of the $`\overline{\mathrm{\Omega }}^+`$/$`\mathrm{\Omega }^{}`$ ratio

The procedure presented in the previous section allows us to estimate the uncorrected $`\overline{\mathrm{\Omega }}^+`$/$`\mathrm{\Omega }^{}`$ ratios in the acceptance given in table 1. The results are given in table 2.

The errors are statistical only. In Fig. 3 (left) the NA49 $`\overline{\mathrm{\Omega }}^+`$/$`\mathrm{\Omega }^{}`$ ratios as a function of the collision energy ($`\sqrt{s_{NN}}`$) are shown and compared with results of WA97/NA57  and STAR . The NA49 and WA97/NA57 results measured at the same energies are consistent. The data show a clear increase of the $`\overline{\mathrm{\Omega }}^+`$/$`\mathrm{\Omega }^{}`$ ratio from a value of about 0.5 at SPS energies to about 1 at RHIC energies. The ratio measured for p+p interactions is similar to the ratio for central Pb+Pb collisions.

In Fig. 3 (right) the antibaryon/baryon ratios are shown as a function of the beam energy in the SPS energy domain. The energy dependence of $`\overline{B}`$$`/`$$`B`$ ratios gets weaker with increasing strangeness content.

### 4.2 $`\mathrm{\Omega }^{}`$ ($`\overline{\mathrm{\Omega }}^+`$) in Pb+Pb at 158 A$``$GeV

At 158 A$``$GeV, a high statistics data sample of central Pb+Pb collisions is available, which allows us to obtain fully corrected spectra of $`\mathrm{\Omega }^{}`$ and $`\overline{\mathrm{\Omega }^+}`$. The transverse mass $`m_t`$ $`=`$ $`\sqrt{p_t^2+m_\mathrm{\Omega }^2}`$ spectra are shown in Fig. 4. They are fitted by the exponential function

$$\frac{1}{m_t}\frac{d^2N}{dm_tdy}=Ce^{m_t/T}$$
(1)
where the fit parameters are a normalization factor C and the inverse slope parameter T. The first data point ($`m_tm_0`$ $`<`$ 0.25 GeV), was excluded from the fit. The slope parameter is similar for $`\mathrm{\Omega }^{}`$ and $`\overline{\mathrm{\Omega }}^+`$ : $`T`$($`\mathrm{\Omega }^{}`$) = 276 $`\pm `$ 23 MeV and $`T`$($`\overline{\mathrm{\Omega }}^+`$) = 285 $`\pm `$ 39 MeV . Our values agree with those measured by the WA97 collaboration .

The large acceptance of the NA49 experiment allows us to measure the $`\mathrm{\Omega }^{}`$ ($`\overline{\mathrm{\Omega }}^+`$) spectra in a large rapidity interval. Fig. 5 shows the rapidity distributions for $`\mathrm{\Omega }^{}`$ (left) and $`\overline{\mathrm{\Omega }}^+`$ (right) extrapolated to $`p_t`$ = 0 using the exponential shown in Fig. 4.

Both spectra were fitted by a Gaussian. The width of the $`\mathrm{\Omega }^{}`$ distribution ($`\sigma `$($`\mathrm{\Omega }^{}`$) = 1.0 $`\pm `$ 0.2) seems to be larger than the one of the $`\overline{\mathrm{\Omega }}^+`$ ($`\sigma `$($`\mathrm{\Omega }^{}`$) = 0.7 $`\pm `$ 0.1). Mean multiplicities in full phase-space were estimated as integrals over measured points corrected for the missing rapidity coverage using the Gaussian parametrisations. The resulting yields are $`<`$$`\mathrm{\Omega }^{}`$$`>`$ = 0.47 $`\pm `$ 0.07 and $`<`$$`\overline{\mathrm{\Omega }}^+`$$`>`$ = 0.15 $`\pm `$ 0.02, where the errors are statistical only.

## 5 Conclusions and outlook

Preliminary results on $`\mathrm{\Omega }^{}`$ and $`\overline{\mathrm{\Omega }}^+`$ production in Pb+Pb collisions and p+p interactions at CERN SPS energies were presented. The $`\overline{\mathrm{\Omega }}^+`$/$`\mathrm{\Omega }^{}`$ ratios increase with energy from about 0.5 at SPS to about 1 at RHIC energies for both central Pb+Pb (Au+Au) and p+p interactions.
The energy dependence of the antibaryon/baryon ratio is weaker for particles with a higher strangeness content. Fully corrected spectra of $`\mathrm{\Omega }^{}`$ and $`\overline{\mathrm{\Omega }}^+`$ were presented at 158 A$``$GeV.

In 2002, Pb+Pb collisions at 20 A$``$GeV and 30 A$``$GeV were taken in the interesting range between top AGS and low SPS energies. Fig. 6 shows the first $`\mathrm{\Omega }^{}`$ \+ $`\overline{\mathrm{\Omega }}^+`$ invariant mass spectrum at 30 A$``$GeV. Analysis of the data at 20, 30 and 80 A$``$GeV is planned for the near future.

## Acknowledgments

This work was supported by the Director, Office of Energy Research, Division of Nuclear Physics of the Office of High Energy and Nuclear Physics of the US Department of Energy (DE-ACO3-76SFOOO98 and DE-FG02-91ER40609), the US National Science Foundation, the Bundesministerium für Bildung und Forschung, Germany, the Alexander von Humboldt Foundation, the UK Engineering and Physical Sciences Research Council, the Polish State Committee for Scientific Research (2 P03B 130 23, SPB/CERN/P-03/Dz 446/2002-2004, 2 P03B 02418, 2 P03B 04123), the Hungarian Scientific Research Foundation (T032648, T14920 and T32293), Hungarian National Science Foundation, OTKA, (F034707), the EC Marie Curie Foundation, and the Polish-German Foundation.

## References


text:
# ”Magnetoscan”: A Modified Hall Probe Scanning Technique for the Detection of Inhomogeneities in Bulk High Temperature Superconductors

## 1 Introduction

Melt textured bulk superconductors are promising candidates for technical applications such as permanent magnets, motors or flywheels. These applications require macroscopically homogeneous properties in the largest possible specimens. However, melt processed materials are strongly inhomogeneous, at least on a microscopic scale , since the formation of defects, e.g. 211 particles, microcracks, twin planes, subgrain boundaries, stacking faults, dislocations, oxygen deficient regions or the boundaries between the growth sectors, cannot be avoided during the growth process. These defects influence the transport properties in various ways. Large normal conducting inclusions simply reduce the superconducting cross section, cracks impede the current flow locally, (sub)grain boundaries reduce the critical current density across them , but some defects can also act as pinning centers and, therefore, improve the transport properties. The defect structure is related to the growth conditions . Differences were found between the a-growth and the c-growth sectors  and the defect concentration generally increases with increasing distance from the seed.

The experimental characterization of the transport properties in large samples is commonly made by scanning the remnant flux profile. This technique is rather insensitive to material inhomogeneities , since all parts of the sample volume contribute to the magnetic field at each measured point. Even if the current distribution can be estimated under certain assumptions , it is generally impossible to assess the correlation between the critical currents and the local microstructure, since $`J_c`$ does not only depend on the microstructure, but also on the magnetic field, which changes within the sample. Therefore, a new technique - the ”magnetoscan” \- has been developed. Here, the currents are induced only locally over a small sample area, where the magnetic field remains nearly constant. Moving the magnet enables us to compare the local transport properties at different positions of the sample with a spatial resolution of less than one millimeter. The magnetoscan, therefore, represents an ideal tool to relate the local $`J_c`$ to microstructural investigations, in order to optimize the growth process in the end.

## 2 The principle: Local shielding currents induced by a permanent magnet

Although the details of the response of a superconductor to a nonuniform magnetic field induced by a permanent magnet are quite complex , shielding can be understood qualitatively, if the permanent magnet is much smaller than the superconductor and if the magnet is placed close to the surface, but not close to the edge of the superconductor.

The following experiment was made. A cylindrically shaped permanent magnet was moved vertically to the surface of a large bulk sample ($`R=12.9mm`$, $`h=10.8mm`$), which had been zero field cooled in liquid nitrogen. The magnet with a radius $`R=3mm`$ and a height $`h=18mm`$ had a magnetic induction of $`100mT`$ at the center of its top surface. After removing the magnet, the remnant flux profile of the sample was scanned by a conventional Hall probe scanner at a distance of $`0.2mm`$ above the surface (Fig. 1a).

Flux penetrated the sample near the edge of the magnet (white circle), whereas no flux is trapped in the center (the inner profile is built up by removing the magnet). Flux penetrates and leaves (bright spots) the sample very inhomogeneously, which must be related to the local properties of the sample.

This observation represents the basis for the method presented in this paper. Flux penetrates more easily in regions with ”bad” transport properties, the response of the sample is smaller there. In regions with high shielding currents, the feedback field of the sample is larger. The flux penetration at the edge of the sample (black circle in Fig. 1a) is not important for our purpose and results from the large demagnetization there. If the magnet is moved horizontally over the sample surface, the flux profile could either move simultaneously or flux could (partly) be trapped in the superconductor. This can be checked by mapping the remnant flux after moving the magnet parallel to the sample surface along the typical path of a Hall scan (Fig. 1b). We find that flux is trapped at the edge of the sample, but (almost) not inside. This indicates that hysteretic memory effects due to the scanning process are absent, apart from the sample edge.

The penetration depth of the magnetic field is roughly estimated from the relation $`\delta =B/\mu _0J`$. With a typical critcal current density of $`10^8Am^2`$, the penetration depth of the field of this permanent magnet is expected to be slightly smaller than $`1mm`$.

## 3 Magnetoscan: Experimental Setup and Tests

Fig. 2 shows the experimental setup schematically. A Hall probe (AREPOC HHP-VPO) with an active area of $`50\times 50\mu m^2`$ is fixed on the edge of the small permanent magnet used above. The permanent magnet and the Hall probe remain fixed, while the sample is moved together with the dewar by a x-y sledge. The gap between the active area of the Hall probe and the sample surface is about $`0.2mm`$, between the permanent magnet and the sample about 1 $`mm`$. Measurements of the Hall voltage were taken on a 0.25 or $`0.5mm`$ grid.

A preliminary scan was made with the Hall probe positioned at the center of the magnet’s surface, not at the edge. In this case the signal nearly did not change apart from the sample edge. This behavior is expected from the above considerations about the magnet’s ”fingerprint”, since the magnetic field was found to be zero at the center of the magnet. Only the transition from nearly perfect screening above the sample to no screening beside it is observed.

Figure 3a shows a typical example of a magnetoscan (with the Hall probe at the edge of the magnet). The arrow indicates the direction of the relative movement of the magnet and the Hall probe. The Hall probe always follows the magnet (when measurements are taken) in order to minimize the remaining hysteretic effects. The change of the signal is much faster, when the Hall probe reaches the sample, since the magnet is already above the sample in this case, while on the opposite edge, the magnet has already left the sample. Due to the rapidly changing geometrical conditions, the signal cannot be evaluated in terms of current transport properties there. A further magnetoscan was made after rotating the sample by 180°. As expected, the geometrical effects near the sample edge remained unchanged, whereas the details of the magnetoscan caused by inhomogeneities were also rotated by 180°.

It should be mentioned that the signal is very sensitive to the gap between the Hall probe and the sample surface. Only samples with a flat surface can be investigated. If the (flat) sample suface is not oriented completely horizontally, the signal increases with decreasing gap, but the inhomogeneities are still detectable, because the systematic influence of such a misalignment and of local changes caused by inhomogeneities can be easily distinguished.

## 4 Results

Figure 3b shows the details of the magnetoscan discussed above (Fig.3a) just for the inner part of the sample. The difference in the Hall voltage between two contour lines is now smaller by a factor of 7 compared to Fig. 3a. The usual trapped flux profile ($`0.2mm`$ above the sample surface) is plotted for comparison in Fig. 3c ($`B_T^{max}=0.54T`$).

Three interesting features can be observed in this magnetoscan: the irregularly shaped contour lines, the smaller response at the position, where the seed was located, and the highest values near the boundaries between the different a-growth sectors (white lines in Fig. 3b).

In order to confirm the relation of these findings to the local transport properties, the conventional trapped field profile was measured $`0.2mm`$ and $`0.4mm`$ above the sample surface. The difference of these two profiles should be approximately proportional to the derivative of $`B_z`$ with respect to the vertical coordinate, i.e. $`B_z/z`$. This quantity is strongly influenced by the surface currents as follows from the Biot-Savart law:

$$𝐁(𝐫)=\frac{\mu _0}{4\pi }\underset{V}{}\frac{𝐉(𝐫^{})\times (𝐫𝐫^{})}{|𝐫𝐫^{}|^3}d^3r^{}$$

The integrand does not depend on $`z`$ alone, but always on $`zz^{}`$. Therefore, in the integral for $`𝐁(𝐫+\delta 𝐞_z)`$ (in the difference quotient) $`𝐫`$’ can be substituted by $`𝐫^{}\delta 𝐞_z`$, and, if $`𝐉/z`$ exists, the partial derivative of $`𝐁`$ with respect to $`z`$ can be expressed by

$`{\displaystyle \frac{𝐁}{z}}(𝐫)`$ $`=`$ $`{\displaystyle \frac{\mu _0}{4\pi }}{\displaystyle \underset{F}{}}{\displaystyle \frac{𝐉(𝐫^{})\times (𝐫𝐫^{})}{|𝐫𝐫^{}|^3}}|_{z^{}=h/2}dx^{}dy^{}`$
$`{\displaystyle \frac{\mu _0}{4\pi }}{\displaystyle \underset{F}{}}{\displaystyle \frac{𝐉(𝐫^{})\times (𝐫𝐫^{})}{|𝐫𝐫^{}|^3}}|_{z^{}=h/2}dx^{}dy^{}`$
$`+{\displaystyle \frac{\mu _0}{4\pi }}{\displaystyle \underset{V}{}}{\displaystyle \frac{𝐉}{z}}|_{𝐫=𝐫^{}}\times {\displaystyle \frac{(𝐫𝐫^{})}{|𝐫𝐫^{}|^3}}d^3r^{}`$
The first and the second term are integrations over the bottom ($`z^{}=h/2`$) and top ($`z^{}=h/2`$) surface of the sample, respectively. If the field profiles are measured close to the top surface, the first term can be neglected. The influence of the third term is a priori unknown. If the variation of the currents along the sample thickness is smooth (at least near the top surface), which ensures the existence of $`𝐉/z`$, its contribution is small, otherwise it could become dominant. If it is negligible, the difference of two field profiles measured at slightly different heights (Fig. 3d), is proporional to the field profile of the surface currents (second term). The qualitative agreement with the magnetoscan is satisfactory. Note that the two profiles are not expected to be identical, since currents all over the surface contribute to the differential profile (plus a contribution of the whole sample), while the currents are induced only locally during the magnetoscan. For instance, a completely homogeneous current density results in a constant signal during the magnetoscan (apart from the edges), but the corresponding differential profile has a sharp peak at the center of the sample (field of a thin sample). Nevertheless inhomogeneities can be detected in both cases. The regions with weak currents near the sample edge (indicated by arrows in Fig. 3d), the higher currents at the boundaries of the a-growth sectors and the reduction of the currents directly below the seed are clearly observed with both methods.

Further examples of magnetoscans are plotted in Fig. 4. The dimension of the sample was $`R=13.1mm`$ and $`h=9.8mm`$, its maximum trapped field after dc activation was $`0.79T`$.

The profile of the top side (a) is nearly quadratic, although the sample is cylindrically shaped. The best shielding is again observed near the growth boundaries (oriented horizontally and vertically) and a depression of the shielding currents is found at the former position of the seed. The weaker shielding near the seed (indicated by the white arrow) corresponds to a bright spot, which is visible at the sample surface. At the bottom side (Fig. 4b) the transport properties are inhomogeneous on a larger scale. By optical inspection of the sample, two types of differently looking regions are observed at the bottom side, i.e. crystalline regions with a pronounced subgrain structure and rather dull black looking areas. The shape of these areas is nicely displayed by the magnetoscan with higher shielding currents in the crystalline regions.

Although many inhomogeneities can be ascribed to optically visible features of the samples (position of the seed, growth sector boundaries etc.), there remain peculiarities in the results, which cannot be correlated with visible inhomogeneities, e.g. the reasons for the irregular (cf. 3b) or quadratic (cf. 4a) shape of the magnetoscan profiles.

## 5 Conclusions

Magnetoscanning represents a simple method to investigate changes of the transport properties of large bulk superconductors with a spatial resolution of around one millimeter. In contrast to the usual scanning Hall probe method, shielding currents are induced only locally. The range of magnetic field penetration can be chosen by the size and the strength of the permanent magnet or could even be varied by using an electromagnet.

## References


text:
# Tunnel magnetoresistance in double spin filter junctions

## I Introduction

In the past few years, tunnel magnetoresistance (TMR) in magnetic tunnel junctions consisting of two ferromagnetic metal (FM) electrodes separated by a thin insulator (FM/I/FM) has attracted much attention due to its promising applications in magnetic field sensors and magnetic random access memory . The TMR in FM/I/FM junctions depends on the degrees of electron spin polarization of the two FM electrodes , but the lack of nearly perfectly spin-polarized current, and temperature dependence of the polarization, have limited the TMR in these junctions . If we use half-metallic electrodes in the magnetic tunnel junctions, which are completely spin-polarized at $`T`$=0 K due to complications such as disorder and surface effects, half-metallicity is destroyed and thus the spin polarization and the TMR are decreased . However, by exploiting the spin filtering phenomenon, one can create 100% spin polarization and obtain a giant TMR.

The spin filtering effect in a ferromagnetic semiconductor (FMS) has been demonstrated dramatically in spin-polarized tunnelling experiments . In these experiments, using an Al/EuS/Al tunnel junction, Moodera et al obtained 85% spin polarization for tunnel electrons. In another study with EuSe barrier junctions  and in the presence of an external magnetic field, they reported near 100% spin polarization. More recently, LeClair et al  from the combination of a spin filter (SF) barrier and a FM electrode, obtained a large TMR in an Al/EuS/Gd tunnel junction, which is a new method for spin injection into semiconductors . When a FMS layer such as EuS is used as a tunnel barrier, due to the spin splitting of the conduction band in the FMS, tunnelling electrons see a spin-dependent barrier height. Thus the probability of tunnelling for one spin channel will be much larger than the other, and a highly spin-polarized current may result. Although the TMR using other SF barriers has also been studied , the results have shown only little success.

The TMR in single and double SF junctions has also been investigated theoretically in recent years. Based on the two band model and free-electron approximation , Li et al  studied the tunnelling conductance and magnetoresistance of FM/FMS/FM junctions. The results showed that a decrease or increase in tunnelling current strongly depends on the magnetization orientations both in the electrodes and in the FMS layer. In another theory with double SF junctions, Worledge et al  studied the TMR of the NM (nonmagnetic metal)/FMS/FMS/NM junctions, using a simple model, and obtained a large magnetoresistance. In a recent paper, Wilczynski et al  investigated tunnelling in NM/FMS/NM/FMS/NM junctions in a sequential tunnelling regime. They found a strong enhancement of the TMR with increasing spin splitting of the barrier height in the ferromagnetic barriers. However, they have not investigated the effect of electric field and temperature on the tunnel currents and spin polarization. Thus, other aspects of this structure, such as the voltage dependence of TMR, remain to be explained.

In the previous paper , using the single-site approximation for the NM/FMS/NM/FMS/NM double SF junction, we studied the effect of the thickness of the central layer on the spin polarization of tunnel electrons at $`T`$=0 K. We showed that the tunnelling spin polarization has an oscillatory behavior with the thickness, which is due to the spin asymmetry of the reflection at the FMS/NM interfaces.

In this paper, using the transfer matrix method and the free-electron approximation, we have extended our previous work  to investigate the effect of the thickness of the central layer, temperature and applied bias on the TMR for tunnelling through a double SF junction. The paper is organized as follows. In section 2, the transfer matrix approach of the spin-polarized tunnelling through the double SF junctions is described. In section 3, the numerical results for the TMR and the spin polarization of the tunnel currents are discussed. The conclusions are summarized in section 4.

## II Model and formalism

The system we consider here is a NM/FMS/NM/FMS/NM sandwich structure in the presence of an applied bias $`V_a`$, which is depicted in Fig. 1. For simplicity, we assume the FMS layers, which act as SFs, are made of the same material, while the outer NM electrodes, which are semi-infinite, and the central layer are made of the same metal. This structure is grown in the $`x`$ direction. In this case, in a free-electron approximation for the spin-polarized tunnelling electrons, the longitudinal part of the effective one-electron Hamiltonian can be written as

$$H_x=\frac{\mathrm{}^2}{2m_j^{}}\frac{d^2}{dx^2}+U_j(x)+V_j^\sigma ,$$
(1)
where $`m_j^{}`$ ($`j`$=1-5) is the electron effective mass in the $`j`$th layer and

$$U_j(x)=\{\begin{array}{cc}0,& x<0,\\ eV_ax/(Lc)+U_2,& 0<x<b,\\ eV_ab/(Lc),& b<x<b+c,\\ eV_a(xc)/(Lc)+U_4,& b+c<x<L,\\ eV_a,& x>L,\end{array}$$
(2)
where $`U_2`$ and $`U_4`$ are, respectively, the barrier heights of the left and right FMS layer at above $`T_C`$, $`b`$ and $`d`$ are the barrier widths, $`c`$ is the width of the middle NM layer and $`L=b+c+d`$. The third term in Eq. (1) which is a spin-dependent potential, denotes the $`sf`$ exchange coupling between the spin of tunnelling electrons and the localized $`f`$ spins in the $`j`$th FMS layer. Within the mean field approximation, $`V_j^\sigma `$ is proportional to the thermal average of the $`f`$ spins, $`S_z`$ (a 7/2 Brillouin function), and can be written as $`V_j^\sigma =I\sigma S_z`$. Here, $`\sigma =\pm 1`$, which corresponds to $`\sigma =,`$, respectively and $`I`$ is the $`sf`$ exchange constant in the FMS layers.

The Schrödinger equation for a barrier layer under the influence of an applied bias can be simplified by a coordinate transformation whose solution is a linear combination of the Airy function Ai\[$`\rho (z)`$\] and its complement Bi\[$`\rho (z)`$\] . Considering all five regions of the NM/FMS/NM/FMS/NM junction shown in Fig. 1, the eigenfunctions of the Hamiltonian (1) with eigenvalue $`E_x`$ have the following forms:

$$\psi _{j,\sigma }(x)=\{\begin{array}{cc}A_{1\sigma }e^{ik_1x}+B_{1\sigma }e^{ik_1x},& x<0,\\ A_{2\sigma }\text{Ai}[\rho _{2\sigma }(x)]+B_{2\sigma }\text{Bi}[\rho _{2\sigma }(x)],& 0<x<b,\\ A_{3\sigma }e^{ik_3x}+B_{3\sigma }e^{ik_3x},& b<x<b+c,\\ A_{4\sigma }\text{Ai}[\rho _{4\sigma }(x)]+B_{4\sigma }\text{Bi}[\rho _{4\sigma }(x)],& b+c<x<L,\\ A_{5\sigma }e^{ik_5x}+B_{5\sigma }e^{ik_5x},& x>L,\end{array}$$
(3)
where the coefficients $`A_{j\sigma }`$ and $`B_{j\sigma }`$ are constants to be determined from the boundary conditions and

$$k_1=\sqrt{2m_1^{}E_x}/\mathrm{},$$
(4)
$$k_3=\sqrt{2m_3^{}(E_x+eV_ab/(Lc))}/\mathrm{},$$
(5)
$$k_5=\sqrt{2m_5^{}(E_x+eV_a)}/\mathrm{},$$
(6)
$`\rho _{j,\sigma }(x)={\displaystyle \frac{x}{\lambda _j}}+\beta _{j,\sigma },j=2,4,`$ (7)
with

$$\lambda _j=\left[\frac{(Lc)\mathrm{}^2}{2m_j^{}eV_a}\right]^{1/3},$$
(8)
$`\beta _{j,\sigma }=\{\begin{array}{cc}\frac{(Lc)[E_xU_jV_j^\sigma ]}{eV_a\lambda _j},& j=2,\\ \frac{(Lc)[E_xU_jV_j^\sigma eV_ac/(Lc)]}{eV_a\lambda _j},& j=4.\end{array}`$ (11)
Although the transverse momentum $`𝐤_{}`$ does not appear in the above notations, the effects of the summation over $`𝐤_{}`$ will be considered in our calculations.

### A Transmission coefficients

By using the boundary condition such that the wavefunctions and their first derivatives are matched at each interface point $`x_j`$, i.e. $`\psi _{j,\sigma }(x_j)=\psi _{j+1,\sigma }(x_j)`$ and $`(m_j^{})^1[d\psi _{j,\sigma }(x_j)/dx]=(m_{j+1}^{})^1[d\psi _{j+1,\sigma }(x_j)/dx]`$, we obtain a matrix formula that connects the coefficients $`A_{1\sigma }`$ and $`B_{1\sigma }`$ with the coefficients $`A_{5\sigma }`$ and $`B_{5\sigma }`$ as follows:

$`\left[\begin{array}{c}A_{1\sigma }\\ B_{1\sigma }\end{array}\right]=M_{total}\left[\begin{array}{c}A_{5\sigma }\\ B_{5\sigma }\end{array}\right],`$ (16)
where

$`M_{total}`$ $`=`$ $`{\displaystyle \frac{k_5}{k_1}}\left[\begin{array}{cc}ik_1& \frac{1}{\lambda _2}\frac{m_1^{}}{m_2^{}}\\ ik_1& \frac{1}{\lambda _2}\frac{m_1^{}}{m_2^{}}\end{array}\right]\left[\begin{array}{cc}\text{Ai}[\rho _{2\sigma }(x=0)]& \text{Bi}[\rho _{2\sigma }(x=0)]\\ \text{Ai}^{}[\rho _{2\sigma }(x=0)]& \text{Bi}^{}[\rho _{2\sigma }(x=0)]\end{array}\right]`$ (38)
$`\times \left[\begin{array}{cc}\text{Ai}[\rho _{2\sigma }(x=b)]& \text{Bi}[\rho _{2\sigma }(x=b)]\\ \frac{1}{\lambda _2}\frac{1}{m_2^{}}\text{Ai}^{}[\rho _{2\sigma }(x=b)]& \frac{1}{\lambda _2}\frac{1}{m_2^{}}\text{Bi}^{}[\rho _{2\sigma }(x=b)]\end{array}\right]^1`$
$`\times \left[\begin{array}{cc}\text{cos}(k_3c)& \frac{m_3^{}}{k_3}\text{sin}(k_3c)\\ \frac{k_3}{m_3^{}}\text{sin}(k_3c)& \text{cos}(k_3c)\end{array}\right]`$
$`\times \left[\begin{array}{cc}\text{Ai}[\rho _{4\sigma }(x=b+c)]& \text{Bi}[\rho _{4\sigma }(x=b+c)]\\ \frac{1}{\lambda _4}\frac{1}{m_4^{}}\text{Ai}^{}[\rho _{4\sigma }(x=b+c)]& \frac{1}{\lambda _4}\frac{1}{m_4^{}}\text{Bi}^{}[\rho _{4\sigma }(x=b+c)]\end{array}\right]`$
$`\times \left[\begin{array}{cc}\text{Ai}[\rho _{4\sigma }(x=L)]& \text{Bi}[\rho _{4\sigma }(x=L)]\\ \text{Ai}^{}[\rho _{4\sigma }(x=L)]& \text{Bi}^{}[\rho _{4\sigma }(x=L)]\end{array}\right]^1`$
$`\times \left[\begin{array}{cc}ik_5& \frac{1}{\lambda _4}\frac{m_5^{}}{m_4^{}}\\ ik_5& \frac{1}{\lambda _4}\frac{m_5^{}}{m_4^{}}\end{array}\right]^1\left[\begin{array}{cc}e^{ik_5L}& 0\\ 0& e^{ik_5L}\end{array}\right]^1.`$
Since there is no reflection in region 5, the coefficient $`B_{5\sigma }`$ in Eq. (3) is zero. In this case the transmission coefficient of the spin $`\sigma `$ electron which is defined as the ratio of the transmitted flux to the incident flux, for the double SF structure shown in Fig. 1, can be written as

$$T_\sigma (E_x,V_a)=\frac{k_5m_1^{}}{k_1m_5^{}}\left|\frac{1}{M_{total}^{11}}\right|^2,$$
(39)
where $`M_{total}^{11}`$ is the left-upper element of the matrix $`M_{total}`$ which is defined in Eq. (38). Note that the transmission coefficient depends on the longitudinal energy $`E_x`$, the applied bias $`V_a`$ and the spin orientation.

### B Spin polarization and TMR

The spin-dependent current density for single or double SF junctions at a given applied bias $`V_a`$ can be calculated within the free-electron model :

$$J_\sigma =\frac{em_1^{}k_BT}{4\pi ^2\mathrm{}^3}_0^{\mathrm{}}T_\sigma (E_x,V_a)\mathrm{ln}\left\{\frac{1+\mathrm{exp}[(E_FE_x)/k_BT]}{1+\mathrm{exp}[(E_FE_xeV_a)/k_BT]}\right\}𝑑E_x,$$
(40)
where $`k_B`$ is the Boltzmann constant, $`T`$ is the temperature, and $`E_F`$ is the Fermi energy.

The degree of spin polarization for the tunnel current is defined by

$$P=\frac{J_{}J_{}}{J_{}+J_{}},$$
(41)
where $`J_{}`$ $`(J_{})`$ is the spin-up (spin-down) current density. For the double SF junction, this quantity can be obtained when the magnetizations of two FMS layers are in parallel alignment.

The tunnel conductance per unit area is given by $`G=_\sigma J_\sigma /V_a`$. In this case, the TMR can be described quantitatively by the relative conductance change as

$$\text{TMR}=\frac{G_{}G_{}}{G_{}},$$
(42)
where $`G_{}`$ and $`G_{}`$ correspond to the conductances in the parallel and antiparallel alignments of the magnetizations in the FMS layers, respectively.

## III Numerical results and discussion

The numerical calculations have been performed for a NM/EuS/NM/EuS/NM double SF junction in which, for simplicity, we assume that the EuS layers have the same thickness $`b=d`$=0.5 nm. The appropriate parameters for EuS which have been used in this paper are: $`T_C`$=16.5 K , $`S`$=7/2, $`I`$=0.1 eV , $`m_{\mathrm{EuS}}^{}=1.5`$ me  and $`U_2=U_4=E_F`$+0.75 eV. In the NM layers, the electron effective mass and the Fermi energy are taken as $`m_{\mathrm{NM}}^{}=m_e`$ and $`E_F`$=1.25 eV. Here $`m_e`$ is referred to the free-electron mass. In this study, we calculate the spin currents, tunnelling spin polarization and the TMR by using Eqs. (40)-(42), respectively. In our considered system, the magnetization orientation (i.e. the $`f`$ spins’ direction) in the left EuS layer stays fixed but the other EuS layer is free and may be switched back and forth by an external magnetic field (see Fig. 1). Thus, when the magnetizations of two FMS layers are parallel, spin-up and spin-down electrons see a symmetric structure, while for the antiparallel alignment these electrons see an asymmetric structure. This structural asymmetry results from the difference in the two barrier heights for each spin channel.

Figure 2 shows the TMR as a function of the thickness of the central NM layer at $`T`$=0, 0.45 and 0.9 $`T_C`$, when the bias voltage $`V_a`$=50 mV is applied to the junction. It is obvious that the TMR oscillates with increasing thickness $`c`$ and have well-defined peaks in which the TMR can reach 99% in some structures. The height of these peaks decreases with increasing $`c`$. The oscillatory behavior is related to the quantum-well states of the central NM layer and the spin-polarized resonant tunnelling. On the other hand, due to the temperature dependence of spin splitting in the EuS layers, the barrier heights become spin-dependent, so that with decreasing temperature, this spin splitting, and thus the TMR, increases.

To understand the physical origin of the TMR and the oscillations, we study the energy dependence of the transmission coefficients through the double SF junction. In Fig. 3, we have plotted the spin-dependent transmission coefficients at $`T`$=0 K for $`c`$=0.50 nm, which corresponds to a flat area between the peaks, and for $`c`$=0.72 nm, which corresponds to a local maximum in the TMR. Because of the quasibound states in the central NM layer, the transmission coefficients reach unity at the resonance peaks which become sharper in the low incident energy region, since in this energy region the resonance levels are more strongly quantized. The results for $`c`$=0.50 nm show that there is one resonance level in the quantum well for both spin orientations and magnetic alignments. All these resonance levels are far from the Fermi energy. However, the transmission coefficient for spin-up electrons in the parallel alignment is higher than the antiparallel alignment, and for spin-down electrons in both alignments. Thus, there is a small difference in the current density $`J`$ (=$`J_{}+J_{}`$) in both magnetic configurations, which gives rise to relatively small TMR at this thickness of the central layer. For $`c`$=0.72 nm the resonance states shift to the lower energy side and a new resonance level slightly below $`E_F`$ appears, which is active only for spin-up electrons in the parallel alignment. Therefore, there is a large difference in the current density $`J`$ in both alignments and consequently gives rise to large TMR. It is clear from the Fig. 3(b) and 3(d) that, for each thickness $`c`$, the resonance levels for spin-up and spin-down electrons in the antiparallel alignment coincide. The reason is that, for the antiparallel alignment, electrons with up (down) spin are easy (difficult) to tunnel into the central NM layer, and difficult (easy) to tunnel out of it; thus, both the spin-up and spin-down electrons see the same effective height of the barriers during the tunnelling process through the whole system.

In Fig. 4 we show the TMR as a function of the applied bias at $`T`$=0 K for the thicknesses $`c`$ correspond to Fig. 3. With increasing bias voltage, the TMR for $`c`$=0.50 nm decreases very slowly because, in this range of the applied voltage, the discrepancy between the conductance for the parallel alignment and that for the antiparallel alignment increases only slightly. On the other hand, the results show that, for $`c`$=0.72 nm at the beginning, the TMR slowly decreases with increasing the bias voltage. However, for the voltages higher than $`V_a`$=80 mV, it quickly decreases. The reason is that, at higher voltages, one of the resonance levels becomes active for both magnetic alignments which drastically reduces the TMR.

It should be noted that, for very low values of the applied bias and the incident energy, a numerical instability is occurred in some of our calculations, which is due to the use of exact Airy functions. This instability is overcome by using the asymptotic forms of Airy functions  and numerical analytical methods.

In Fig. 5 we have shown the spin polarization of the tunnelling current versus normalized temperature $`T/T_C`$ for single and double SF junctions to reveal the SF effects of the FMS layers from another point of view. At high temperatures $`T>T_C`$, there is no spin splitting $`E_{ex}`$ in the conduction band of the EuS layers and the transmission coefficients for two spin channels coincide. Thus there is no TMR and spin polarization effect. As the temperature decreases from $`T_C`$, the barrier heights for spin-up electrons are lowered, while they are raised for spin-down electrons. This temperature dependence of the barrier height, which is attributed to the exchange splitting of the EuS conduction band, greatly increases the tunnelling probability for one spin channel and reduces it for the other. In the parallel alignment, the tunnel current for spin-up electrons is much higher than the spin-down ones, which gives rise to the TMR and spin polarization effect. On further decreasing the temperature, this spin splitting and hence the difference in the barrier heights increases. Therefore, the TMR and spin polarization reach the highest values at $`T`$=0 K. The highest value of the spin polarization for the single SF junction can reach 77%, which is qualitatively in agreement with the experimental measurements  and the theoretical results , while for the double SF junction (in the parallel alignment), it approaches 66% for $`c`$=0.50 nm and 99% for $`c`$=0.72 nm , which is due to the change in the positions of the spin-polarized quasibound states in the quantum well. Therefore, one can see that, for the double SF junctions, the TMR and the spin polarization of the tunnelling current can be controlled by the thickness of the central NM layer, the temperature, and the applied bias.

In this study, the numerical results obtained for the TMR and the degree of spin polarization can be compared with the result of Wilczynski et al . As we discussed above, the enhancement of spin splitting $`E_{ex}`$ in the conduction band of the FMS layers can be obtained, when the temperature decreases. In this case, the TMR and the spin polarization increases, and this behaviour is in agreement with the result obtained by Wilczynski et al in the sequential tunnelling regime.

## IV Summary

Based on the free-electron model, we presented a transfer matrix method for spin-polarized tunnelling through the double SF junctions. The effect of the quantum size, applied bias and temperature on the spin filtering and the spin transport process in the FMS/NM double junctions are examined. Numerical results indicate that the TMR oscillates as the thickness of the central NM layer increases. It is further confirmed that, for some thicknesses of the central NM layer, the spin-polarized resonant tunnelling can gives rise to large values for the spin polarization and the TMR, even at high temperatures ($`T<T_C`$). Therefore, in the system presented here, it is able to select an appropriate applied bias, temperature and the thickness to achieve a maximum TMR and spin polarization.

Although the temperature for observing a 99% TMR is very low and the findings of the paper are not directly applicable to spintronics technology, the results may be useful in designing future spin-polarized tunnelling devices .

In the present model, we have neglected the generally important complications, such as the interface roughness, electron-electron interaction, magnetic-domain wells, $`f`$ spin correlation, etc. These effects can play important roles in the spin transport process and reduce the efficiency of spin filtering.


text:
# Zeroth principle of thermodynamics in aging quasistationary states

## Acknowledgments

We thank C. Anteneodo for useful remarks and discussions. We have benefitted from partial support by CNPq, CAPES, PRONEX and FAPERJ (Brazilian agencies).