text:
# Discovery of a Cyclotron Resonance Scattering Feature in the X-ray Spectrum of XTE J1946+274

## 1 Introduction

The transient accreting X-ray pulsar XTE J1946+274 was discovered during an outburst in 1998 September with the All Sky Monitor (ASM) on board the *Rossi X-ray Timing Explorer* (RXTE) (Smith & Takeshima, 1998). At the same time, 15.8 second pulsations were observed with the CGRO/BATSE, with the pulsating source designated GRO J1944+26 (Wilson et al., 1998). The best X-ray position (BeppoSAX/MECS: 1 radius, 90% confidence, Campana et al., 1998) falls within the Ariel V error region for the transient source 3A 1946+274, making it likely that the two objects are identical (Campana, Israel & Stella, 1999).

The initial outburst (see Fig. 2), which reached $``$6 cps in the RXTE/ASM, lasted about 100 days and was followed by a series of smaller flares. Campana, Israel & Stella (1999) discussed the extended RXTE/ASM lightcurve, finding that the flaring was nearly periodic with a repetition period of $``$80 d. They interpret this as either the half or full binary orbital period. Either case is consistent with the known orbital period range of Be/X-ray binary pulsars. Given its probable orbital period and the outburst characteristics, XTE J1946+274 is most likely another example of a Be star/X-ray binary pulsar transient. This class of binaries accounts for over half of the known accreting pulsars (Liu, van Paradijs & van den Heuvel, 2000).

Several RXTE observations (see Figure 2) were made spanning the peak of the initial outburst. In this *Letter*, we report on the spectral analysis of these observations and, in particular, the discovery of a cyclotron resonance scattering feature, or “cyclotron line”, at $``$35 keV. Santangelo et al. (2001) report their independent discovery of the XTE J1946+274 cyclotron line with BeppoSAX in this volume.

Cyclotron lines result from the scattering of X-rays by electrons in quantized Landau orbits in the $``$1012 G fields near the magnetic poles of accretion powered pulsars. The characteristic energy of the Landau transition scales with the magnetic field as $`E_{cyc}=(11.6\text{ keV})(1+z)^1B_{12}`$, where $`E_{cyc}`$ is the cyclotron line energy in keV, $`z`$ is the gravitational redshift in the scattering region, and $`B_{12}`$ is the magnetic field in units of $`10^{12}`$ G. Because of this proportionality, cyclotron lines give us the only direct measure of the neutron star magnetic field. In general, harmonically spaced lines (corresponding to higher order Landau transitions) may exist (see for example Heindl et al., 1999, 2000; Santangelo et al., 1999). Depending on the temperature and geometry of the emitting and scattering material and the viewing angle with respect to the magnetic field, the line profiles may be broad and complex (e.g. Araya-Góchez & Harding, 2000; Kretschmar et al., 2000). To date, about a dozen pulsars have well-established cyclotron lines, for the most part discovered with Ginga, RXTE, and BeppoSAX (see for example Makishima et al., 1999; Heindl et al., 2000; Dal Fiume et al., 2000).

## 2 Observations

During the XTE J1946+274 outburst, 12 pointed observations were made with the Proportional Counter Array (PCA) and the High Energy X-ray Timing Experiment (HEXTE) on-board the RXTE (see Tab. 1 and Fig. 2). The PCA (Jahoda et al., 1996) is a set of five Xenon proportional counter units (PCUs) sensitive in the energy range 2–60 keV with a total effective area of $``$7000 $`\mathrm{cm}^2`$. The HEXTE (Rothschild et al., 1998) consists of two arrays (“clusters A and B”) of 4 NaI(Tl)/CsI(Na) phoswich scintillation counters (15–250 keV) totaling $``$1600 $`\mathrm{cm}^2`$. Early in the mission, the pulse height analyzer of a single cluster B phoswich failed, making the effective area of cluster B approximately 3/4 of cluster A. The HEXTE clusters alternate pointing between target and nearby blank fields in order to measure the background. The PCA and HEXTE fields of view are collimated to the same 1 full width half maximum (FWHM) region. In order to extend the life of the PCA detectors, most observations are currently made with one or more PCUs turned off. All the observations here were performed with PCUs 0, 1, and 2 on. Because PCUs 3 and 4 were sometimes off, they have been excluded from this analysis.

## 3 Analysis

For each observation, we accumulated source and background spectra from the PCA standard 2 data and the HEXTE event mode data. We used the FTOOLS 5.0.1 and the procedures detailed by Wilms et al. (1999). In order to take advantage of the best PCA calibrations, we limited the data to detector layer 1 only. We summed the spectra for the 3 active PCUs to form a single pulse height spectrum. We used the bright source background model (appropriate for the counting rate of XTE J1946+274) to estimate the PCA background. Because of the high PCA counting rate of XTE J1946+274 and the duration of the observations, the statistical errors in the energy spectra ($`1\%`$) were small compared to the uncertainties in the instrument calibration. It was therefore necessary to include systematic uncertainties on the PCA data in order to achieve a reasonable fit to the data. We determined the size of these errors through the procedures recommended by the PCA team (Jahoda, 2000) and explained in Coburn et al. (2001). In short, we fit contemporaneous observations of the Crab Nebula and Pulsar to a two power law model, using the XTE J1946+274 response matrices. This model allows for different power law indices for the pulsar and nebular components (Knight, 1982). We adjusted the systematic errors until a reduced chi-squared of one was achieved – at a level of 0.3% per channel below 20 keV. Above 20 keV, the systematic errors became unacceptably large, $``$5%, so we chose to ignore these PCA data, rather than apply large systematic errors.

We accumulated the HEXTE cluster A and B spectra separately then summed them to form joint source and background spectra for each observation. For the analysis of these summed spectra, we added the cluster A and B response matrices, weighted by the cluster effective areas (4:3, respectively). Fits with this matrix to a single power law model, appropriate in the HEXTE energy band, to the Crab have residuals of $`\pm `$1%. This is smaller than the statistical uncertainties, so we applied no additional systematic errors to the HEXTE data.

Because the observations were spread over $``$1 month and a range of $``$40% in flux, we were concerned that changes in the source spectrum could affect our analysis. To search for spectral changes, we calculated for each observation the ratio of the net PCA and HEXTE counts spectra to the total spectrum of all observations (see Figure 3). In all but observation 12, we

found that above8 keV, the ratio was nearly flat, indicating that only the flux, and not the spectral shape, had changed. In observation 12 (the last and dimmest observation – see Table 1), significant spectral changes were observed, and this pointing was eliminated from our analysis. Observations 1 to 11 showed no changes in spectral shape at a level of $``$1–2% between 8–20 keV. Below 8 keV, systematic changes of up to $``$5% were

seen, particularly in the lowest energy channels. These changes are consistent with $``$10% variations in the absorbing column and the iron line emission. Such variations are not unexpected during a Be/X-ray binary outburst, and they are likely not indicative of changes in the underlying continuum emission. The HEXTE ratios showed no evidence for spectral variability at a limiting level of $``$5–10% below 35 keV. Above 35 keV the statistics did not allow for a meaningful comparison. For these reasons, we concluded that by excluding data below 8 keV we could safely combine observations 1 to 11.

After accumulating observations 1 through 11 as described above, we attempted to fit the PCA and HEXTE spectra jointly with a set of continuum models typically used for accreting pulsars. While there is a notable lack of adequate theoretical models, pulsar spectra are heuristically well-described by a power-law at low energies which breaks to an exponentially cut-off power-law at high energies (White, Swank & Holt, 1983). Kreykenbohm et al. (1999) review several models with this asymptotic behavior – the exponentially cut-off power law (PLCUT, White, Swank & Holt, 1983), the Fermi-Dirac cut-off (FDCUT, Tanaka, 1986), and a combination of two power-laws with a positive and a negative exponent (NPEX, Mihara, 1995). None of these models provided an acceptable fit to the XTE J1946+274 spectrum. All models left significant, cyclotron line-like residuals near 35 keV (see Figure 3). We therefore added to each model a cyclotron absorption term with a Gaussian optical depth profile:

$$\tau (E)=\tau _{\text{cyc}}\times e^{(EE_{\text{cyc}})^2/2\sigma _{\text{cyc}}^2}.$$
(1)
Adding the Gaussian absorption line greatly improved the fits. The best fit was achieved with the NPEX continuum, given by:

$$\text{NPEX}(E)(E^{\mathrm{\Gamma }_1}+\alpha E^{+\mathrm{\Gamma }_2})\times e^{(E/E_{\text{fold}})}$$
(2)
where $`\text{NPEX}(E)`$ is the photon flux at energy $`E`$; $`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_2`$ are the indices of the falling and rising power law components, with $`\alpha `$ their relative normalizations; and $`E_{\text{fold}}`$ is the exponential folding energy. Our best-fit model was of the form $`\text{Flux}\text{NPEX}(E)\text{e}^{\tau (E)}`$. In this case, the reduced $`\chi ^2`$ changed (with the addition of the line) from 3.61 to 0.97 for 52 and 49 degrees of freedom respectively. The resulting $`F`$-Test probability for this to be a chance improvement is $`1.2\times 10^{14}`$. Table 2 gives the best fit parameters for this model.

## 4 Results and Discussion

Figure 3 shows the best fit model and residuals together with the inferred incident photon spectrum. It also shows residuals to the best NPEX model without a cyclotron line, which is then apparent in the residuals around 35 keV. From the cyclotron line energy, we deduce a magnetic field strength in the scattering region of $`3.1(1+z)\times 10^{12}`$ G. The line is weakly resolved with the HEXTE energy resolution of $``$8 keV (FWHM, at 35 keV), and the fitted width of the line is $`9.3_{2.1}^{+2.5}`$% of the line energy. This is within the range of typical values determined with the Gaussian absorption profile model (Coburn, 2001). No harmonic line near 70 keV was found, but as is evident in Figure 3, the falling high energy continuum provides inadequate statistics to make a sensitive search for such a feature.

Figure 4 (after Coburn, 2001) shows the spectral cut-off energy plotted against the cyclotron line energy for 8 pulsars measured with RXTE. Also shown is the (somewhat controversial) HEXE/TTM result for A0535+26 (Kendziorra et al., 1994). For uniformity, the pulse-phase average spectra of all sources were fit with the PLCUT continuum:

$$\text{PLCUT}(E)E^\mathrm{\Gamma }\times \{\begin{array}{cc}1\hfill & \text{if }EE_{\text{cut}}\text{;}\hfill \\ e^{(EE_{\text{cut}})/E_{\text{fold}}}\hfill & \text{otherwise}\hfill \end{array}$$
(3)
smoothed at $`E_{cut}`$ (see Coburn, 2001, for details). Two sources are excluded from the plot: 4U 1626-67, whose cut-off energy varies by a factor of four with pulse phase (Coburn, 2001), and the low luminosity ($``$$`4\times 10^{34}`$$`\mathrm{ergs}\mathrm{s}^1`$) object 4U 0352+309 (X Per) whose spectrum is not well fit by the usual pulsar models. Due to the complex shape of the 4U 0115+63 fundamental line (Heindl et al., 2000), its plotted energy is half of the first harmonic value.

A clear correlation exists between the line energy and spectral cut-off indicating that the cut-off is related to the magnetic field. Makishima & Mihara (1992) and Makishima et al. (1999) first noted this relationship by plotting cyclotron line energies from Ginga and other instruments, derived with a variety of continuum models, against the PLCUT cut-off energy from sometimes non-contemporaneous Ginga observations. They found that the relationship was consistent with a power law, $`E_{cyc}E_{cut}^{1.4}`$ (or, $`E_{cut}E_{cyc}^{0.7}`$), indicating a saturation as compared to a linear correlation. Figure 4 has the advantage that each point is derived from a uniform model fit to a single spectrum. And, with the exception of A0535+26, all points are from the same set of instruments. XTE J1946+274 fits nicely on the correlation, and below 30 keV, the slope is consistent with $`E_{cut}E_{cyc}^{0.7}`$. However, there appears to be a break in the relationship above $``$30 keV. This flattening is more abrupt than the smooth turnover of the $`E_{cyc}^{0.7}`$ power law suggesting that the processes that form the continuum saturate at higher magnetic fields.

With XTE J1946+274, we have added a 13th accreting pulsar to the list of objects with secure cyclotron line detections. Nearly all of these have been confirmed, and several were discovered, with RXTE and BeppoSAX. We have now begun detailed studies of these objects as a class. In particular, since ten of the thirteen sources have been observed with RXTE, we are applying uniform analyses to all these objects to further understand the line forming regions. First results of such studies, including correlations between other line parameters, are given in Coburn (2001).

ASM data are provided by the RXTE/ASM teams at MIT and at the RXTE SOF and GOF at NASA’s GSFC. This work was supported by NASA grant NAS5-30720, NSF travel grant INT-9815741, DFG grant Sta 173/25-1, and a travel grant from the DAAD.


format: text

text:
# Exchange Bias in Ferromagnetic/Compensated Antiferromagnetic Bilayers

## Abstract

By means of micromagnetic spin dynamics calculations, a quantitative calculation is carried out to explore the mechanism of exchange bias (EB) in ferromagnetic (FM)/compensated antiferromagnetic (AFM) bilayers. The antiferromagnets with low and high Néel temperatures have been both considered, and the crossover from negative to positive EB is found only in the case with low Néel temperature. We propose that the mechanism of EB in FM/compensated AFM bilayers is due to the symmetry broken of AFM that yields some net ferromagnetic components.

Exchange anisotropy was first discovered in 1956 by Meiklejohn and Beanfirst , who found that the hysteresis loop of Co/CoO after cooling in a magnetic field was no longer centered at zero field (H=0) but was shifted along the field axis. The shifted direction was found to be opposite to the applied magnetic field (negative exchange bias $`H_B<0`$) and the magnitude of this shift is known as exchange bias (EB). It was subsequently established that this might be a general phenomenon for any ferromagnet (FM)/antiferromagnet (AFM) systems cooling in an applied magnetic field (cooling field) from above the Néel temperature(TN) of the AFM, with the FM Curie temperature(TC) greater than TN. In recent years, since the phenomenon of exchange bias has become the basis for an important application in information storage technology Application , tremendous efforts have been made for exploring the mechanism Review , Review2 .

Meiklejohn and Bean originally suggested that exchange bias was a consequence of the presence of interfacial uncompensated AFM spins. In view of this argument, a natural question to ask is whether the exchange bias also exists in a FM/ compensated AFM system. Surprisingly, in a compensated $`Fe/FeF_2`$ bilayer system, Nogues et al. observed not only the usual negative exchange bias but also an unexpected positive exchange bias ($`H_B>0`$) under large cooling fields Positive .

Several important theories have existed to study the exchange bias in compensated AFM. KoonKoon  presented a microscopic explanation of EB due to a irreversible AFM domain wall, and found a perpendicular orientation between the FM/AFM axis directions, namely spin-flop state. With consideration of magnetostatic interactions in this spin-flop state, Schulthess et al. Butler  obtained the opposite results, i.e., not EB but a large uniaxial anisotropy, and attributed EB to the interfacial defects. Unfortunately they did not further show the EB-magnetostatic interactions phase diagram, in other words, EB should gradually change with magnetostatic interactions. HongHong  argued that interface spin configuration persisted after cooling below TN, and negative/positive bias respectively corresponded to parallel/perpendicular easy axes of FM and AFM. Kiwi et al.Kiwi  suggested a canted AFM spin configuration frozen into a metastable state, and proposed the incomplete FM domain wall model to explain positive exchange bias. However, the former two theories were carried out with micromagnetic calculations without consideration of the cooling field; the later two theories pointed out the cooling field without micromagnetic calculations, and were lack of much more detailed and sufficient microscopic information. Up to now, exchange bias mechanism is still controversial.

In this paper, based on the assumptions that an antiferromagnetic interface coupling between FM/AFM is responsible for exchange bias in FM/compensated AFM Positive , Postive2 \- Coupling , and that the biased hysteresis loop is basically determined by the spin configurations in the underlying antiferromagnetic layer after cooling Nature -asymmetry2 , we carry out micromagnetic calculations using spin dynamics to explain the mechanism of EB in the FM/compensated AFM systems. Physically the key point different from previous Koon’s micromagnetic calculationsKoon , addresses the cooling field during cooling process. We succeed to reproduce both the negative and positive EB effects. Qualitatively speaking, it is a competition among (i) the cooling magnetic field (ii) the interface coupling of FM/AFM, and (iii) the spin-spin interaction and anisotropy of AFM, that eventually determines the spin configurations in AFM during the cooling process. For an AFM with weak spin-spin interaction (low TN), the spin configuration of AFM at low cooling magnetic field is dominated by the AF type interface coupling of FM/AFM. Therefore, the initially compensated AFM layers especially the interface AFM layer becomes weekly uncompensated, resulting in a net ferromagnetic component opposite to the cooling field (or the magnetization in FM). The hysteresis loop is then measured at low temperature after removing the cooling field, while the spin configuration in AFM is frozen. Similar to the arguments given by Meiklejohn and Bean, it can easily be deduced that the broken symmetry of AFM in this case favors the negative exchange bias. However on the other hand, if the cooling field $`H_{CF}`$ is large and becomes to dominate, then a net ferromagnetic component along the cooling field is expected. Because of the AFM type interface coupling of FM/AFM, it turns out that the broken symmetry of AFM in this case favors the positive exchange bias. On other hand, for AFM with strong spin-spin interaction ( high TN ) only negative EB can be found in a reasonable high $`H_{CF}`$. The quantitative results are given in the following to reveal in details how these different terms affect the broken symmetry of AFM layers and its correlation with the exchange bias.

Our model Hamiltonian is

$$H=H_{AA}+H_{FF}+H_{AF},$$
(1)
where $`H_{AA}`$ is the part of AFM layers, $`H_{FF}`$ and $`H_{AF}`$ the FM and the interface coupling between AFM and FM layers. They are

$$H_{AA}=\underset{}{}J_{AA}𝐒_i𝐒_jD_A\underset{i}{}\left(S_i^x\right)^2H_{CF}\underset{i}{}\left(g_Au_B\right)S_i^x,$$
(2)
$$H_{FF}=\underset{}{}J_{FF}𝐒_i𝐒_jH_{CF}\underset{i}{}\left(g_Fu_B\right)S_i^x,$$
(3)
$$H_{AF}=\underset{}{}J_{AF}𝐒_a𝐒_f,$$
(4)
where $`g_A`$, $`g_F`$, $`u_B`$, $`D_A`$ and $`H_{CF}`$ denote AFM Lande factor, FM Lande factor, Bohr magneton, antiferromagnetic anisotropy and cooling field in parallel with AFM anisotropy, respectively. The exchange coupling among spins is considered for nearest neighbor sites only. The subscripts $`a`$ and $`f`$ are associated with AFM and FM respectively. It is noticed that the anisotropy of FM layer is neglected based on the fact that most experiments did use the soft ferromagnets. The dipole-dipole interactions in the system are not considered here, since it affects only quantitatively rather than qualitatively on the symmetry broken of AFM. As the previous models, we also assume that $`J_{AF}`$ $`J_{AA}`$ Koon , Butler . The Néel temperature increases monotonically with $`J_{AA}`$, thus the interface coupling is stronger in the FM/AFM systems with higher TN and vice versa.

Now we calculate the EB by the following spin dynamics approachButler , i.e.,the local effective field is determined from the gradient of the energy, $`𝐇_i^{\mathrm{𝐞𝐟𝐟}}=\frac{H}{gu_B𝐒_i}`$, and $`\left\{𝐒_i\right\}`$ is required to satisfy the Laudau-Lifshitz equation of motion with the Gilbert-Kelley form for the damping term: $`\frac{}{t}𝐒_i=gu_B𝐒_i\times \left(𝐇_i^{\mathrm{𝐞𝐟𝐟}}\eta \frac{}{t}𝐒_i\right)`$, where $`\eta `$ denotes the damping parameter. This damping term is phenomenological and is included to remove the energy from the system and to ensure that the magnetic system is in a stable or metastable equilibrium after sufficient iterating calculating steps. A lattice with $`50\times 50\times 2(FM)`$ and 10 layers of $`AFM`$ is used in our calculation. In the beginning, the temperature $`T`$ of the system is set at $`T_c>T>T_N`$. Therefore, the initial spin configuration in our spin dynamics calculation is such that the spins are randomly arranged for AFM but are ferromagnetic arranged for FM. A cooling field is then applied along the $`x`$ direction that is also the easy axis of AFM. Meanwhile it is known that the spins of FM will be easily aligned according to the cooling field. After taking a long step of spin dynamics calculation, a stable state of FM/AFM under the cooling field is finally approached. Then the system is cooled down to the low temperature, we switch off the cooling field and start to do simulation of hysteresis loop of FM layers while the spin configuration of AFM is fixed.

Fig.1 shows EB as a function of cooling field for a AFM/FM system with low TN such as $`FeF_2`$ $`(T_N78.4K)`$ or $`MnF_2`$ $`(T_N67.3K)`$. In doing this we set the parameters in Hamiltonian as $`g_A=g_F=2.0,`$ $`J_{FF}=10mev,`$ $`J_{AA}=0.8mev,`$ $`J_{AF}=J_{AA}/2=0.4mev`$ and $`D_A=0.4mev`$ per siteReview , Koon , Butler . As a natural output from the calculation, it is indeed observed in this figure that the exchange bias $`H_B`$ changes sign from negative to positive as the cooling field increases, and a crossover field $`H_{cross}`$ is found at about $`3.7T`$. Dashed and solid lines in the inset show the negative and positive loops at 2KOe and $`7T`$ respectively.

For systems with low TN, i.e., weak spin-spin interaction $`J_{AA}`$ in AFM, the spin configuration of AFM at low cooling magnetic field is dominated by AFM type interface coupling of FM/AFM. In this case, the symmetry broken of compensated AFM layers appears. Some net ferromagnetic component along the $`x`$ axis is expected, which means that the broken symmetry of AFM in this case favors the negative exchange bias. However for higher cooling magnetic field $`H_{CF}`$, the cooling field becomes to dominate the broken symmetry so that a net ferromagnetic component along the positive $`x`$ direction is expected, i.e., the broken symmetry of AFM in this case favors the positive exchange bias.

For a quantitative description of the FM components in AFM, We define the ferromagnetic component in the n-th AFM layer

$$S^x(n)=\underset{i}{}S_{na}^x(i)/N_n,$$
(5)
where $`S_{na}^x(i)`$ is the AFM spin at site $`i`$ of the n-th layer, $`N_n`$ is the number of lattices in the n-th AFM layer, while the first layer is defined as the interface layer of AFM. This quantity $`S^x(n)`$ describes the degree of symmetry broken in each layer of AFM. It is found from our calculation that the ferromagnetic components are layer dependent. As expected, it will be larger when the layer is near the interface, and become smaller when the layer is far from the interface. In Fig.2 the ferromagnetic component of the interface AFM layer $`S^x(1)`$ is shown as a function of $`H_{CF}`$, using the same parameters as obtaining Fig.1. Similar results as Fig.1 are found, that $`S^x(1)`$ is negative at the beginning when the cooling field is small, then reaches to zero at a critical field $`(37KOe)`$ and finally becomes positive as the field further increasing. Since the first AFM layer should be the most important one in the interface coupling, it is reasonable to see that the ferromagnetic component of the interface AFM layer should be responsible for the EB effect. The inset of Fig.2 gives the layer dependent ferromagnetic components when the cooling field is fixed at 2KOe. The oscillation is caused by the antiferromagnetic exchange interaction between the layers of AFM.

Fig.3 shows the same relationship but for a AFM/FM system with higher $`T_N`$ such as FeMn $`(500K)`$. The parameters used here are $`J_{FF}=10mev,`$ $`J_{AA}=5mev,`$ $`J_{AF}=J_{AA}=5mev`$ and $`D_A=3mev`$ per site Review . In this case, it is found that $`H_B`$ is always negative and changes a little when $`H_{CF}`$ is ranged from 2T to 7T. This result also agrees qualitatively with experimentReview . Inset(a) of Fig.3 shows the magnetization loop at 3T cooling field. In fact that $`J_{AA}`$ is large when the AFM layer of the system has high $`T_N,`$ thus the interface coupling $`J_{AF}`$ also becomes large $`\left(J_{AF}J_{AA}\right)`$, then the AFM type interface coupling controls the symmetry broken of AFM. In this case, if $`H_{CF}`$ is reasonable high $`(2T7T)`$ but not too high, $`H_B`$ is found to be always negative. This can explain that the positive EB was reported only in the FM/AFM thin films with low TNPositive , Postive2  \- Coupling .

One distinguished feature of Fig.3 is that: a tip EB is found at cooling field $`H_{tip}`$ indicated by arrow, and this tip also is clearly shown in Inset(b). As previous mention, the preceding discussions are subjected to both cooling field and applied magnetic field parallel with AFM easy axis. In fact, the orthogonal FM/AFM spin configuration similar to Koon’s conclusion can also be recovered with zero or smaller cooling field for stable spin configuration. With the increasing cooling field from zero to $`H_{tip}`$, the FM spins will gradually rotate direction from the perpendicular to parallel to the AFM easy axis during cooling process, and the interface coupling contribution to the negative bias will enlarge and nearly saturate at $`H_{tip}`$. On the other hand, with the increasing cooling field above $`H_{tip}`$, the cooling field contribution to the potentially positive bias will raise, in other words, contribution to the negative bias will lessen. Thus for cooling field parallel with AFM easy axis situation, there exists a tip EB associated with cooling field $`H_{tip}`$.

In summery, micromagnetic spin dynamics calculations are carried out to explain the mechanism of EB in the FM/compensated AFM system. Different from the previous micromagnetic calculations, we address the key role of cooling field. Some important experiment results, such as the cooling field dependent transition from negative to positive EB in the AFM/FM layers with low TN, can be reproduced. It is proposed that the symmetry broken of AFM plays a key role in explaining the exchange bias.

This work is supported by the National Natural Science Foundation of China and Shanghai Research Center of Applied Physics. One of the authors X.F. Jin would acknowledge the supports of the Cheung Kong Scholars Program, Hong Kong Qiushi Science Foundation, Y.D. Fok Education Foundation.


text:
# Strange Hadron Resonances and QGP Freeze-out

## 1 Possible hadronization scenarios

An important aspect in the study of heavy ion collisions is to determine the time scale governing hadron production and the duration of the decoupling (freeze-out) process. At present, there is no general consensus on the freeze-out time scale, and hence on how much observed particles “remember” about their primordial source. quark-gluon plasma (QGP) signals are visible in hadronic particles, when freeze-out is explosive, with sudden breakup of the fireball, and little interaction between the newly-formed particles . In principle the newly-formed hadrons could undergo a period of re-interaction in a hadronic gas phase. This re-interaction phase, the time scale of which can be hadron-specific, could significantly alter any considered QGP signal .

The study of short-lived hyperon resonances produced in heavy ion collisions, either by direct observation , or through their effect in observed heavy ion $`m_T`$ distributions , provides a very promising way to distinguish between these scenarios, since natural lifetime of hyperon resonances is generally shorter than the expected duration of the rescattering hadronic gas (HG) phase. Therefore, the resonances signal and distribution will be sensitive to the duration (or indeed the existence) of this HG phase. In particular, rescattering of resonance decay products should deplete the observable resonance signal, since short-lived resonances are detected by invariant mass reconstruction. A significant amount of rescattering would re-thermalize resonance decay products, which affects the resulting particle’s $`m_T`$ distributions. This paper analyses existing resonance data, in an atempt to distinguish between the sudden and staged hadronization scenarios.

In next section 2 we look the spectra of hyperons and show that non-rethermalized resonance decays are much better agreeing with the observed $`m_{}`$-shapes. We than study the directly observed strange hadron resonances in section 3. We describe the possibility of resonance quenching in subsection 3.1, and compare in a simple model of an opaque medium with the experimental results in subsection 3.2.

## 2 The role of resonances in fits

The precise $`m_T`$ spectra produced by the WA97 experiment  provided the first opportunity to test the sudden hadronization model. A necessary consequence of explosive hadronization is that chemical freeze-out, which fixes particle abundances, and thermal freeze-out, which fixes particle spectra, coincide. Hence, it should be possible to fit both normalisation and $`m_T`$ spectra for several particles, using the same temperature, flow and chemical potentials. The fitted parameters should also agree, within error bars, with the chemical fit parameters, obtained by using just particle ratios. Our simultaneous fit to normalized WA97 data was successful .

The sample of $`K_s,\mathrm{\Lambda },\overline{\mathrm{\Lambda }},\mathrm{\Xi },\overline{\mathrm{\Xi }},\mathrm{\Omega }+\overline{\mathrm{\Omega }}`$ (divided into four centrality bins) was fitted, using a single set of thermodynamic parameters for each bin, to a Cooper-Frye distribution , truncated to eliminate unphysical emission within a spacelike freezeout hypersurface

$$\frac{dN}{m_Tdm_T}=V\underset{i}{}\lambda _i\gamma _i𝑑\sigma _\mu p^\mu \mathrm{\Theta }(d\sigma _\mu p^\mu )f(T,v,m_T).$$
(1)
Here V is the fireball volume (i.e., an experiment-specific normalisation parameter), $`\lambda _i,\gamma _i`$ are the flavor fugacities and saturation parameters , $`\sigma ^\mu `$ and $`p^\mu `$ are, respectively the 4-momentum and the freeze-out hypersurface, and $`f(T,v,m_T)`$ is the Lorentz-transformed Boltzmann distribution, dependent on temperature, flow and transverse mass.

In such an analysis , $`T`$ and the chemical parameters agreed, within the error bars, with earlier abundance fits , flow and volume increased with centrality, as expected for a thermalize expanding fireball, while the hadronization surface velocity was consistently close to the speed of light, in agreement with an explosive hadronization scenario. $`\chi ^2`$ profiles showed a clear, statistically significant minimum in each of the quantities fitted.

The only systematic deviation between statistical hadronization model and data is present in the $`\mathrm{\Omega }`$ spectrum. It’s unusually large enhancement and anomalous slope are also interpreted as evidence for a staged freeze-out, with $`\mathrm{\Omega }`$s freezing out earlier than other particles. However, as shown in , it is only the low momentum $`\mathrm{\Omega }`$s that deviate significantly from the expected behavior in statistical hadronization. Such a low-momentum yield enhancement is explained by an additional QGP specific physical process, such as strange diquark correlations , or the formations of disoriented chiral condensates .

For the purpose of the study of the time scales in statistical hadronization an important direct constraint is provided by the prominence, in these fits, of resonance decay products. A non-negligible fraction of hyperons (more than half in case of the $`\mathrm{\Lambda }`$) within the spectrum fitted are produced not directly, but through decay of produced resonances. Yields of heavy resonances are thermally suppressed because of their higher mass, but are typically enhanced by their higher spin and isospin degeneracy. If the resonance decays with an appreciable center of mass momentum (as is the case in all hyperon resonances, except $`\mathrm{\Sigma }^0\mathrm{\Lambda }\gamma `$), the distribution of particles originating from decays differs considerably from a thermal distribution. Therefore, the shape of $`m_{}`$ spectra which includes resonance decay products will depend on whether the decay products of short-lived resonances have had the opportunity to rethermalize or not. If the freeze-out is sudden, the observed spectra are a mixture of direct and decay shapes, and if there is a phase in which the particles produced during hadronization can undergo re-interaction, the spectrum observed is thermal, but expansion flow deformed.

In order to obtain a final particle spectrum from a thermal distribution of resonances, it is necessary to transform the decay products momenta from the resonance rest frame to the fireball rest frame, and integrate over the kinematically allowed momentum space. This procedure has been implemented in several ways by different authors (  ); While numerically cumbersome, it requires no extra degrees of freedom within the fit.

The test for the reaction mechanism is obtained from the study of the $`\chi ^2`$ distribution. A consistent lowering of $`\chi ^2`$ once resonances are taken into account in a hyperons $`m_T`$ distribution confirms that resonance decay products emerge from the fireball unthermalized. In contrast, if resonance decay products re-thermalize, their effect will be a correction to normalisation (chemical potentials and volume) rather than slope parameters. In such a case, using the decay distributions , will not result in a better fit to the data.

Fig. 1 shows the temperature $`\chi ^2`$ profiles (for the four WA97 centrality bins, with bin 4 being most central), calculated using the prescription in  to obtain decay spectra (solid lines) and assuming that resonance decay products rethermalize (dashed lines). These profiles are computed by minimising every parameter except the temperature (flow, hadronization surface velocity and chemical potentials). Thin lines are results obtained allowing unphysical range of parameters, thick lines with the constraint on the light quark saturation parameter $`\gamma _q`$, corresponding to the condensation limit of the $`\pi `$: $`\gamma _\mathrm{q}^2e^{\frac{m_\pi }{2T}}`$, where $`m_\pi `$ is the pion mass. The maximum value of $`\gamma _\mathrm{q}`$ is prefered by the hadronizing high entropy QGP phase .

We see that the condensation constraint does not modify the value of the $`\chi ^2`$ minimum for the solid lines (which correspond to non-thermalized resonance decay products), but greatly reduces it’s ambiguity, and is ruling out a freeze-out temperature lower than $`145`$ MeV. As Fig. 1 shows comparing thick solid with the thick dashed lines, the hyperon $`m_{}`$ fit strongly favours a scenario in which the decay products of resonances emerge from the fireball without re-interaction. If the condensation limit on $`\gamma _\mathrm{q}`$ is not enforced (thin lines in 1), the $`\chi ^2`$ is still consistently lowered by the inclusion of unrescattered resonances. On the other hand, once condensation limit is taken into account (thick lines), a model with thermalize resonance decay products statistically fails to describe the data: $`\chi ^2`$ increases by a factor of 3, and the physically sensible minimum for thermal freeze-out disappears.

We therefore conclude that the WA97 $`m_{}`$ data strongly favours a sudden hadronization scenario in which the hadrons emerging from the fireball without much re-interaction after freeze-out. However, a more direct way to measure the hadronization time scale is necessary in order to firmly confirm this result. We hope and expect that the experimental measurement of the resonance signal will prove useful in this regard.

## 3 Direct probes of sudden hadronization

### 3.1 First results on resonance yields

The resonance signal can be measured experimentally using invariant mass reconstruction. At present, the $`K^0`$ has been observed by both NA49  and STAR , and NA49 has also measured the $`\mathrm{\Lambda }(1520)`$. It appears that while the $`K^0`$ and $`\overline{K^0}`$ are produced in the expected amounts, the $`\mathrm{\Lambda }(1520)`$ yield is severely depleted. In fact, so far it has proven to be the only hyperon suppressed with respect to a proton-proton collision.

The observation of this suppression, by itself, has not been enough to settle the freeze-out time scale controversy. In principle such a suppression is welcome by both camps. The proponents of sudden freeze-out argue that freeze-out happens at a temperature considerably lower than that of the equilibrium QCD phase transition . The yield N of a resonance of mass $`m`$ at temperature $`T`$ obeys,

$$Nm^2TK_2(m/T)(Tm)^{3/2}e^{m/T},$$
(2)
here $`K_2`$ is the modified Bessel function. We see that heavy resonances such as the $`\mathrm{\Lambda }(1520)`$ are more suppressed, since temperature is lower. The other camp notes that if one assumes a staged freeze-out, the suppression comes from rescattering of resonance decay products within HG phase. However, in this case it is difficult to see how the $`K^{}`$ is not suppressed, especially considering it’s larger width.

The actual situation is in our opinion more complicated . The $`\mathrm{\Lambda }(1520)`$ resonance is a D-wave of the ground $`\mathrm{\Lambda }`$ (Isospin=1/2) state. It is therefore especially susceptible to in-medium mixing with the lower-mass hyperons (P-wave of the $`\mathrm{\Sigma }`$, isospin $`I=3/2`$ state), through reactions such as

$$\pi +\mathrm{\Lambda }(1520)\mathrm{\Sigma }^{}\mathrm{\Lambda }\pi ,\mathrm{\Sigma }\pi .$$
(3)
We therefore believe that other resonances will need to be measured in order to extract freeze-out information from their abundance. In particular, the non-suppression of $`K^{}`$ s makes other P-wave resonances good candidates for investigation. The $`\mathrm{\Sigma }^{}(1385)`$ is the least massive resonance of this type, and for this reason, as well as it’s other characteristics which we will explore in the next subsection, we regard it as a promising candidate.

### 3.2 Model and data comparison

Within statistical hadronization picture a particles yield depends on it’s mass, the freeze-out temperature, and the chemical potentials corresponding to the particle’s quantum numbers. If experimental acceptance does not cover the full $`4\pi `$ range, flow and hadronization surface velocity also need to be taken into account to determine the partial yields.

Since we are interested here primarily in the temperature, as well as the evolution subsequent to hadronization, it is indicative to pick an observable which is independent of as many of the above quantities as possible. To make chemical potentials cancel out, we can examine the ratio of the resonance to the “ground state” hyperon which contains the same valence quarks. For example, ratios such as $`\mathrm{\Lambda }(1520)/\mathrm{\Lambda }`$$`\mathrm{\Sigma }^{}/\mathrm{\Lambda }`$, $`K^0/K^{}`$ and $`\overline{K^0}/K^+`$ are all independent of chemical potentials. Kaons require additional attention, because the $`K^{}`$’s lifetime is too short for oscillation, but the $`K_s`$ does oscillate. Therefore, the observed $`K^0`$ and the $`K^0`$ do not have the same valence quark composition. We also found that such ratios, and more generally, ratios of particles with comparable mass, are independent of flow and freeze-out surface geometry to a very good approximation.

To get a more quantitative estimate of how temperature and rescattering combine in fixing resonance ratios, a microscopic model that describes the rescattering of resonance decay products in a baryon-rich medium needs to be considered . In such a model the probability of the resonance to decay in a hadron-rich medium, as well as the interaction cross-section between the decay products and the in-medium particles is allowed for. The number of reconstructible resonance decay products is obtained solving,

$`{\displaystyle \frac{dN_i}{dt}}`$ $`=`$ $`\mathrm{\Gamma }N_N^{}N_i\mathrm{\Sigma }_j\sigma _{ij}\rho _{j0}({\displaystyle \frac{R_0}{R_0+vt}})^3,`$
$`{\displaystyle \frac{dN_N^{}}{dt}}`$ $`=`$ $`\mathrm{\Gamma }N_N^{}.`$ (4)
Here $`N^{}`$ and $`N_i`$ are, respectively, the resonance and decay product abundances, $`\mathrm{\Gamma }`$ is the resonance width, $`\rho _{j0}`$ is the density of particle $`j`$ in the hadron gas at freeze-out and $`\sigma _{ij}`$ is the thermally averaged interaction cross-section. Final observable yields of particle $`i`$ are given by integrating this model from chemical freeze-out (with a certain temperature) for an amount of time corresponding to the interacting system’s lifetime.

Because of the rapid decrease in hadron density within the fireball due to flow (commonly estimated at around 0.5c), resonance ratios are extremely sensitive to the lifetime of the interacting phase compared to the resonances width. The chemical freeze-out temperature is also important, since it controls both the initial resonance abundance (through Eq. 2) and the initial density of the hadron gas.

The model is however remarkably insensitive to it’s other parameters, such as the initial radius $`R_0`$ (which is also constrained by the entropy per baryon, $`40`$ in SPS experiments), and the cross-sections. In particular, we found that an order of magnitude change in the interaction cross-sections will result in ratios changing by 30 %. For this reason, we believe that, while our model is very simplified (it neglects inelastic scattering and resonance creation within a hadron gas, and uses thermally averaged cross-sections), the results presented here are qualitatively correct, and serve as a good indication of how the interesting quantities may be extracted from experimental measurements.

We now turn to discuss the available experimental results. Fig. 2 gives the dependence of $`K^{}/K`$ and $`\mathrm{\Lambda }(1520)/\mathrm{\Lambda }`$ on the chemical freeze-out temperature and the lifetime of the interacting hadron gas phase. It is apparent that these two quantities determine the ratios of two resonances of different masses and decay widths, as seen in the figure: an experimental point on a diagram such as Fig. 2 is just what is needed to distinguish between the explosive and gradual freeze-out scenarios.

The non-suppression of $`K^{}`$, together with the very strong suppression of $`\mathrm{\Lambda }(1520)`$, would seem to indicate instantaneous freeze-out with a temperature of less than 100 MeV. However, in obtaining this result $`\mathrm{\Lambda }(1520)`$ quenching was neglected. If 50 $`\%`$ of the $`\mathrm{\Lambda }`$ are quenched through reactions such as those in Eq. (3), existing data would be perfectly compatible with the sudden freeze-out picture, as we show in Fig. 3.

Observation of the $`\mathrm{\Sigma }^{}(1385)`$ yield should remove the uncertainty introduced by the quenching phenomenon. We have shown in , that the masses and lifetimes of the $`K^{}`$ and $`\mathrm{\Sigma }^{}(1385)`$ combine to constrain the observed ratios to a very narrow yield band. Significant deviation from this band could be a strong sign of broadening of either of these two resonances. In addition, because of a coincidence between the masses of the $`\mathrm{\Sigma }^{}`$, the $`\mathrm{\Xi }`$ and the availability of $`\mathrm{\Xi }^{}(1530)`$, the explicit dependence of the $`\mathrm{\Sigma }^{}/\mathrm{\Xi }`$ ratio on temperature is very weak . For this reason, if fugacities and saturation parameters of the light and strange quark are determined independently , comparison between $`\mathrm{\Sigma }^{}(1385)`$, $`\mathrm{\Xi }`$ and $`\mathrm{\Lambda }`$ should be sufficient to disentangle temperature and fireball lifetime conclusively.

We use these properties in Fig. 4, where we show how the $`\mathrm{\Sigma }^{}/\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }^{}/\mathrm{\Xi }`$ ratios depend on the freeze-out temperature and lifetime. We see that lines corresponding to a given freeze-out temperature depend weakly on the details of the rescattering model. This makes the $`\mathrm{\Sigma }^{}/\mathrm{\Xi }`$ ratio a probe of the hadronization temperature. A severely depleted $`\mathrm{\Sigma }^{}/\mathrm{\Xi }`$ ratio, in the absence of evidence of broadening, would unambiguously signal a staged freeze out in which resonance decay products have ample opportunity to rescatter. The $`\mathrm{\Sigma }^{}`$ measurement could therefore distinguish between a scenario based on the explosive freeze-out from a super cooled QGP, and a high T chemical freeze-out, followed by a reinteraction period.

## 4 Conclusions

We have shown considering the statistical significance of different approaches to fit the $`m_{}`$ spectra of WA97 experiment , that the strategy which does not allow for resonance products to rescatter is greatly favoured. To find this, we needed to limit the magnitude of the chemical quark-pair non-equilibrium to the maximum value compatible with the high entropy content, which is hadronizing. The $`\chi ^2`$ minima are very pronounced and thus offer strong evidence for sudden mechanism of hadronization. This finding strongly supports the picture of suddenly hadronizing QGP as the source of these particles.

A more direct study of hadronization lifespan is derived from the observation of hyperon resonance yields. We have shown that combining results of several resonances allows to derive directly information on freeze-out temperature and duration. We have also argued that P-wave resonances are more suitable, than others, for this exercise, as they are less susceptible to in-medium modifications. Measurement of the $`\mathrm{\Sigma }^{}(1385)`$ yield, and it’s ratio to hyperons such as the $`\mathrm{\Lambda }`$ and the $`\mathrm{\Xi }`$, has the potential of being just what is needed to unambiguously extract the freeze-out time scale and temperature from experimental heavy ion data.

### Acknowledgements

We thank Ch. Markert and L. Csernai for valuable comments. Work supported in part by a grant from the U.S. Department of Energy, DE-FG03-95ER40937.

## References


text:
# Vortex structure in mesoscopic superconductors

## I Introduction

Recent advances in nanoscience have demonstrated that fundamentally new physical phenomena are found when the size of the sample shrinks and becomes comparable to the length scale of the investigated phenomenon.

Superconductivity is a macroscopic quantum phenomenon and therefore it is interesting to see how this quantum state is influenced if the sample is reduced in size. For rings and hollow cylinders this leads to the well-known Little-Parks effect  which results in a periodic variation of the critical temperature as function of the applied magnetic field, the period being determined by the magnetic flux value within the tube.

Bulk superconductors are divided up into type-I ($`\kappa <1/\sqrt{2}`$) and type-II ($`\kappa >1/\sqrt{2}`$) superconductors, the distinction between them is completely determined by the Ginzburg-Landau parameter $`\kappa =\lambda /\xi `$, where $`\lambda `$ is the magnetic field penetration depth and $`\xi `$ the coherence length, which is a material parameter. The difference is clearly seen in the magnetic response of the system, where type-I has a complete Meissner effect while type-II superconductors can have a state with partial expulsion of the field in which flux lines penetrate the superconductor. The latter is a superconducting state with disconnected circular areas of normal state. For mesoscopic superconductors we find that type-I superconductors can behave like type-II or even show some mixed behaviour depending on the size of the system and therefore $`\kappa `$ is no longer the only determining parameter characterizing the vortex state of the system.

The vortex structure in superconductors will be described in the framework of the Ginzburg-Landau (GL) theory, which consists of two coupled nonlinear differential equations:

$$\left(i\stackrel{}{}_{2D}\stackrel{}{A}\right)^2\mathrm{\Psi }=\mathrm{\Psi }\left(1\left|\mathrm{\Psi }\right|^2\right),$$
(1)
$$\mathrm{\Delta }_{3D}\stackrel{}{A}=\frac{d}{\kappa ^2}\delta \left(z\right)\stackrel{}{j}_{2D},$$
(2)
where

$$\stackrel{}{j}_{2D}=\frac{1}{2i}\left(\mathrm{\Psi }^{}\stackrel{}{}_{2D}\mathrm{\Psi }\mathrm{\Psi }\stackrel{}{}_{2D}\mathrm{\Psi }^{}\right)\left|\mathrm{\Psi }\right|^2\stackrel{}{A},$$
(3)
is the density of superconducting current and $`d`$ the thickness of the sample. The superconducting wavefunction satisfies the boundary conditions $`\left(i\stackrel{}{}_{2D}\stackrel{}{A}\right)\mathrm{\Psi }|_n=0`$ normal to the sample surface and $`\stackrel{}{A}=\frac{1}{2}H_0\rho \stackrel{}{e}_\varphi `$ far away from the superconductor. The distances are measured in units of the coherence length $`\xi =\mathrm{}/\sqrt{2m^{}\alpha }`$, the order parameter in $`\mathrm{\Psi }_0=\sqrt{\alpha /\beta }`$ and the vector potential in $`c\mathrm{}/2e\xi `$. $`\kappa =\lambda /\xi `$ is the Ginzburg-Landau parameter, and $`\lambda =c\sqrt{m/\pi }/4e\mathrm{\Psi }_0`$ is the penetration length. The magnetic field is measured in $`H_{c2}=c\mathrm{}/2e\xi ^2=\kappa \sqrt{2}H_c`$, where $`H_c=\sqrt{4\pi \alpha ^2/\beta }`$ is the critical field. Temperature is included in $`\xi `$, $`\lambda `$, $`H_{c2}`$, through their temperature dependencies: $`\xi (T)=\xi (0)/\sqrt{|1T/T_{c0}|}`$, $`\lambda (T)=\lambda (0)/\sqrt{|1T/T_{c0}|}`$, $`H_{c2}(T)=H_{c2}(0)|1T/T_{c0}|`$, where $`T_{c0}`$ is the critical temperature at zero magnetic field.

Different approaches have been used to solve the GL-equations. In the lowest Landau level approximation  a linear combination of solutions of the linearized first GL-equation is used and the internal (taken homogeneous) magnetic field and the expansion coefficients are determined by minimizing the free energy. An extension of this approach beyond the lowest Landau level was given in Ref. . It was shown recently  that in order to account for demagnetization effects in the lowest Landau level approximation it is necessary to introduce an effective $`\kappa `$ which depends on the size of the superconductor. Trial functions  with several variational parameters were used to study the vortex configuration in mesoscopic cylinders with suppressed surface superconductivity. Here, we will solve the two GL equations numerically using a finite difference technique. Details of this approach can be found in Ref. .

## II Giant vortex to multi-vortex transition

First we consider a circular disk with thickness $`d`$. Because of the circular symmetry we assume that the order parameter is cylindrical symmetric: $`\mathrm{\Psi }(x,y)=F(\rho )e^{L\varphi }`$ where $`L`$ is the vorticity. In doing so we restrict our set of solutions to a subset in which the modulus of the local order parameter is axially symmetric.

Because of the non-linear term in the GL-equation the superconducting state is in general non-axial symmetric even if the sample is circular symmetric. A well-known example is the Abrikosov lattice which has triangular symmetry. We found that for our circular disks such non-axial symmetric states are found when the radius of the disk is sufficiently large and when the magnetic field is not too large. For small disks the boundary condition dominates which imposes its symmetry on the order parameter. For large magnetic fields the inner part of the disk becomes normal and the superconducting state survives only near the sample surface where the shape of it will again determine the symmetry of the order parameter. The free energy of the superconducting disk is shown in Fig. 1(a) as function of the magnetic field for different values of the vorticity. The full curves are for the giant vortex state and the dashed curves for the multi-vortex states. The transition point between them is given by the open dot. Figs. 1(c,d,e) shows a contour plot of the modulus of the superconducting density for $`L=3`$ at $`H_0/H_{c2}=0.62`$, $`0.72`$ and $`0.82`$, respectively. Notice that with increasing magnetic field the size of the vortices grows and they move closer to each other. At the transition field the single vortices coalesce into one giant vortex. This is a continuous transition and therefore of second order .

The magnetization is shown in Fig. 1(b) where the vertical lines indicate the ground state transitions. Notice that for small magnetic fields we have a linear $`MH`$ relation which is typical for an ideal diamagnet. With increasing magnetic field there is a continuous penetration of the magnetic field at the edge of the sample which leads to a smaller than linear increase of $`M`$ with $`H`$. This effect is enhanced by demagnetization effects which leads to an enhanced magnetic field at the edge of the sample (see inset of Fig. 1(a)). When we further increase the magnetic field the energy of the superconductor increases up to the point where it is energetically more favorable to transit to the $`L=1`$ state. As a consequence $`M(H)`$ exhibits a zig-zag behavior which was measured recently by A. Geim et al  and explained in Ref. . This behavior is due to the fact that with increasing magnetic field more field can penetrate into the sample in a discontinuous way in which the vorticity of the order parameter increases with one unit. Approximately (because of boundary effects), one unit of flux enters the superconductor. Therefore, by counting the number of jumps in the $`M(H)`$ curve it is possible to count the vorticity of the sample and therefore the number of vortices inside the disk.

Notice that for a given magnetic field different states are possible. Experimentally, it has been possible to drive the system into the different metastable regions and to map out the complete magnetization-magnetic field curves as shown in Fig. 1(b) by slowly ramping the field up and down. These metastable states are responsible for hysteretic behavior which was observed in the magnetic field response of superconducting disks. A detailed comparison between experimental results and our theoretical calculations can be found in Ref. . It was shown recently  that in certain samples the superconducting state can be forced very far into the metastable region such that fractional flux can penetrate and in some cases even negative flux penetration is observed for the transition $`LL+1`$.

## III Geometry dependence

One may wonder how important the exact shape of the sample is for the symmetry of the order parameter. To investigate this, we took as an example a thin flat triangular shaped superconducting sample. For the free energy we found a similar behavior as depicted in Fig. 1(a) but with the following major differences: 1) the critical magnetic field $`H_{c3}`$ at which the normal state is reached is substantially larger, i.e. for a triangle with the same surface area as the disk of Fig. 1(a) we obtained $`H_{c3}/H_{c2}=2.5`$. This is a consequence of the enhanced surface conductivity in wedge shaped samples ; 2) The multi-vortex state (see Fig. 2(a)) is substantially stabilized, i.e. the magnetic field range over which it is stable is strongly enhanced. Furthermore, for certain values of the vorticity only the multi-vortex state is stable, and the giant vortex state does not occur. The giant vortex state is now no longer circular symmetric, but it is characterized by the fact that the order parameter has only one zero. 3) A new vortex structure appears which is a combination of a giant vortex in the center of the triangle with single vortices around it. This state occurs e.g. for $`L=5`$; 4) From the contour plot of the superconducting current flow (see Fig. 2(b)) we see that the screening currents along the edge flow clockwise while the currents around the vortices flow counterclock wise. Notice that near the edges their are clockwise spiralling local vortex-like currents which are not connected to the position of a vortex (see the contour plot of Fig. 2(c) for the phase of the order parameter) but are rather back-flow currents which are well-known in hydrodynamics.

Recently, Chibotaru and co-workers  studied the linearized GL-equation for triangles and squares and found that the symmetry of the boundary imposes a similar symmetry on the superconducting state. For example, in case of the triangle and for $`L=2`$ this occurs by having three vortices sitting in a triangular arrangement and one anti-vortex in the center of this triangle. From the outset we should point out that the linearized GL-equation is only strictly valid at the superconducting/normal boundary where the modulus of the order parameter approaches zero. It was found recently  that once one moves away from this phase line the anti-vortices are quickly annihilated. Then a state is found in which a large central area of the triangle is normal and only superconductivity survives in the corners of the triangle . The anti-vortices and the vortices are found in the central area of the triangle where $`|\mathrm{\Psi }|^2`$ is extremely small. Therefore, the experimental relevance of this new state is questionable.

## IV Different arrangements of multi-vortices

When the superconducting state is in the multi-vortex state, it is possible that different vortex configurations are possible. This is illustrated in Fig. 3 where for a disk of radius $`R/\xi =6`$ the energy of the $`L=6`$ and $`L=7`$ vortex states are plotted as function of the magnetic field. In this case a multi-vortex configuration is possible with (dashed curves) or without (solid curves) a vortex in the center of the disk. The corresponding contourplots of the superconducting density in the disk is shown in the insets of Fig. 3. The states are indicated by $`(n;L)`$ where $`n`$ refers to the number of vortices in the center and $`L`$ is the vorticity of the superconducting state. Notice that for $`L=6`$ the (1;6) state has a larger energy than (0;6) when $`H/H_{c2}>0.45`$ which then becomes the ground state. This change of configuration at $`H/H_{c2}.45`$ is a first order transition because it involves a change of symmetry of the superconducting state. At $`H/H_{c2}0.6`$ the ground state transits to a higher vorticity state with configuration (1;7). This transition is also a first order transition in which the magnetization is discontinuous. Such transitions have been found experimentally .

Notice that the appearance of the different multi-vortex configurations is very similar to the classical system of interacting (repulsion) particles which are confined into a potential. It was found that the particles are situated in ring-like configurations , very similar to the vortex configurations found in Fig. 3. Also in this case different meta-stable states consisting of different particle configurations are found.

## V Flux penetration and expulsion

In the above discussion we considered only the local minima in the energy functional which determine the ground and metastable states. The transition between states with different vorticity does not always follow the ground state. The reason is that their are barriers for flux penetration and explusion. The most well-known is the Bean-Livingston model  for the surface barrier which is a result of the competition between the vortex attraction to the sample walls by its mirror image and its repulsion by screening currents. For nonelliptical samples a geometrical barrier appears because of Meissner currents flowing on the top and bottom surface . Additionally, vortex pinning by defects can play an important role in the delay of vortex expulsion or promotion of vortex penetration. These models describe the vortex formation far from the sample boundary.

As an example, we will consider here mesoscopic disks in which boundary effects are predominant and consequently previous approaches are not applicable. We consider a defect-free thin superconducting disk with a perfect circular boundary such that $`Rd/\lambda ^21`$. When calculating the free energy we not only search for the local minima but also for the saddle points which are the lowest energy barriers between two such minima. The details of the calculation can be found in Refs. .

The results are shown in Fig. 4 where the energy of the saddle point is given by the dotted curve. In contrast to known surface and geometrical barrier models, we find that in a wide range of magnetic fields below the penetration field, the saddle point state for flux penetration into a disk does not correspond to a vortex located near the sample boundary, but to a region of suppressed superconductivity (see Figs. 4(a-d)) at the disk edge with no winding of the current (Fig. 4(c)), and which is a nucleus for the following vortex creation. The height of this nucleation barrier is shown in the inset of Fig. 4 and determines the time of flux penetration.

## VI Conclusions

The vortex state of a mesoscopic superconductor is strongly determined by its size and to a lesser extend by the material parameters the superconductor is made of. For example, by increasing the radius (starting from $`R\xi `$) of a superconducting disk  it is possible to obtain a magnetic response which, as function of the magnetic field, is continuous, type-I, a type-I with multiple steps and type-II. The exact geometry of the superconductor has also a strong influence on its vortex state. This vortex state can be brought into a metastable region owing to the presence of barriers for flux motion which leads to unexpected effects like fractional flux and even negative flux entry. In mesoscopic superconductors the flux in general is not quantized . This is even more so in small superconducting rings . The lowest energy barrier between two flux states was identified and corresponds to a saddle point of the energy functional.

## Acknowledgments

This work was supported by the Flemish Science Foundation (FWO-Vl), the ”Onderzoeksraad van de Universiteit Antwerpen” (GOA), the ”Interuniversity Poles of Attraction Program - Belgian State, Prime Minister’s Office - Federal Office for Scientific, Technical and Cultural Affairs”, and the European ESF programme on ’Vortex Matter’. Discussions with S. Yampolskii, V. Schweigert, V. Moshchalkov, A. Geim and L. Chibotaru are gratefully acknowledged.


text:
# 1 Introduction

## 1 Introduction

Efforts to understand the origin of cosmic rays at any energy are greatly hampered by our lack of knowledge of the mass distribution in the incoming cosmic ray beam. The determination of the arrival direction of cosmic rays does not depend on knowledge of the mass of the primary and, using the fluorescence technique, the primary energy is obtainable with only a small systematic uncertainty because of the unknown mass. Even with the traditional ground arrays it has been possible to devise ways of deducing the primary energy that are reasonably independent of model and mass uncertainties, at least at the 30 % level. However use of the data on the energy spectrum and arrival direction distribution to decide between various origin models does require knowledge of the primary mass distribution. While there is a common assumption that protons dominate at the highest energies, hard experimental evidence is lacking, as it is at energies above 1018 eV.

Some models of cosmic ray origin lead to the conclusion that photons dominate the cosmic ray beam. This hypothesis can be more readily tested than can competing hypotheses that advocate protons or iron nuclei as the majority primaries. This is because the shower development of a photon primary is very different from that of any hadronic primary. The difference is larger than is expected between the showers created by proton and iron primaries. For example, Halzen et al.  have shown that a photon is extremely unlikely to produce a shower in the atmosphere that has a development profile that looks like the famous event of $`3\times 10^{20}`$ eV observed by the Fly’s Eye Group . A study of horizontal showers recorded at Haverah Park has shown that at 1019 eV there can be no more than 40 % of the primaries that are photons , while a search, using the traditional idea of photon showers being deficient in muons, has given a similar limit at the same energy . These analyses all assume that photons $`>`$ 1019 eV have not acquired hadron-like properties, as has sometimes been speculated.

There is no reliable evidence on the hadronic mass composition above 1019 eV. Some years ago an analysis of the Fly’s Eye data  pointed to a change from an iron-dominated composition at $`3\times 10^{17}`$ eV to a proton-dominated composition near 1019 eV. This conclusion was drawn from a study of the variation of depth of maximum with energy (the elongation rate) and from an analysis of the spread in depth of maximum at a given energy. The elongation rate behaviour was in agreement with an earlier, model-independent, analysis from Haverah Park  and with a related study using Cherenkov light by the Yakutsk group . The conclusions on mass composition from the study of the spread in depth of maximum were confirmed by  but none of these analyses extended to energies above 1019 eV.

A different conclusion has been reached by the AGASA group based on the variation of the muon content of showers with energy . Their analysis favours a composition that remains ’mixed’ over the 1018 to 1019 eV decade. Furthermore there is a difference between the conclusions reached from the data on depth of maximum from the Fly’s Eye experiment  and those reached from the HiRes prototype, operated with the MIA detector, in the range 1017 to 1018 eV . The more recent data are used to infer protons as early as 1018 eV, if the QGSJET model is correct. There is a good evidence from Haverah Park data  that the QGSJET model is satisfactory, at least to 1018 eV. However the HiRes/MIA conclusion is, in turn, in contradiction with a new analysis of Haverah Park data  that suggests that protons make up only 34 % of the flux, independent of the energy, between $`3\times 10^{17}`$ and $`2\times 10^{18}`$ eV, with iron nuclei making up the remainder.

The addition of more data in recent years and a better understanding of what might be the optimum shower model has not helped to clarify the position: rather the reverse. In figure 1 we show a compilation of data on depth of maximum, $`X_\mathrm{m}`$, as a function of energy . The sensitivity of a conclusion to the choice of model is evident. The matter is thus far from being resolved and there is clearly scope for alternative approaches.

## 2 Outline of a new approach

It seems unlikely that a single method of deriving the mass composition will be sufficient to resolve this important issue. In this paper we argue for a new approach that is based on the study of inclined showers, and which can be added to the armoury that can be used. The method is derived from the work that led to our estimate  of a limit to the photon flux at 1019 eV, an estimate that has now been confirmed independently , as mentioned above. Our approach stems from the detailed understanding of the propagation of inclined showers that has been developed  and from its application to the prediction of the rate of triggering of the Haverah Park array . Inclined showers induced by very high energy cosmic rays produce particle densities at ground level which are dominantly due to muons, in contrast to vertical showers that have an important component due to photons and electrons stemming from $`\pi ^0`$ decay. Measurements of the cosmic ray rate at different zenith gives a handle on the relative number of muons in a shower, which is dependent on composition.

To find the photon limit we assumed that the energy spectrum was known, and calibrated in a mass-independent way, through the work with the Fly’s Eye detector. Using this known ’vertical’ spectrum, a prediction was made of the number of showers expected with energies above 1019 eV for photon, proton and iron primaries. We found strong evidence to reject a photon intensity at the level predicted by several super-heavy relic models ,,. Recent reassessment of the energy spectrum in the appropriate range , and the differences reported between the AGASA group and the HiRes group ,, make it premature to be too firmly in favour of protons or iron above 1019 eV, but our conclusion about photons remains robust. An extensive description of our method used to derive the photon limit has been given recently . In this paper we show how the technique can be used to infer the mass at energies between 1018 and 1019 eV. The conclusions, as to the mass variation with energy that we can presently draw, is constrained by uncertainties in the form of the spectrum in this energy range, and by the assumed hadronic interaction model, rather than by statistical limitations. We pressent conclusions for a number of spectra, using the QGSJET model  as a reference for hadronic interactions, in the expectation that the spectrum uncertainties will soon be resolved.

## 3 The analysis of inclined showers

The analysis of showers with zenith angles $`>`$ 600 is not straightforward, even with an array like Haverah Park in which the detectors were deep-water Cherenkov detectors. These detectors present a significant area to inclined showers by comparison with thin scintillators. There are two major difficulties. Firstly the showers lose their near-circular symmetry very rapidly above 600  and, secondly, the finite size of the array makes it difficult to locate the core. During the operational phase of the Haverah Park array (1967 -1987) a decision was made to ignore events $`>`$ 600 for these reasons. More recently, with the advent of the analytical techniques of Ave et al.  and with the availability of much more computing power, it has proved possible to get accurate directions for the events, to estimate their energies (subject to assumptions about shower models) and to compute the flux. In this section we justify these claims.

The methods of finding the directions and of calculating the energy are described in detail in  for events above 1019 eV. Here we have applied the same techniques, namely an iterative process in which arrival directions are fitted to the time distributions, followed by shower energy and core position fits to the detector densities. The timing of the signals use the information on core position to account for curvature corrections. Following the same steps as described in , we have performed three quality cuts to eliminate uncertainties in the reconstruction procedure: (i) the distance from the central triggering detector to the core position in the shower plane is required to be below $`r_{\mathrm{max}}=500`$ m, (ii) the $`\chi ^2`$ probability for the energy and direction fits must be $`>`$ 1 %, (iii) the downward error in the energy determination is required to be less than a factor of 2. The chosen value of $`r_{max}`$ guarantees that the core position is always close to the triggering detectors in the array. Similar results, but with reduced statistics, was obtained with $`r_{\mathrm{max}}=300`$ m. This last cut is the only one that differs from those used in : it ensures that the smaller showers, that are the concern of this paper, are near to sufficient detectors to allow reliable density fits. After making the cuts described above we found 385 events with $`E_0>10^{18}`$ eV and 2 events with $`E_0>10^{19}`$ eV.

Again mimicking the procedure described in Ref. , we have generated a sample of artificial events to which we have applied the same reconstruction algorithms as used for the data. In figure 2 we show the normalised $`\chi ^2`$ distributions for the timing fit to the 385 events above 1018 eV and for the energy estimates made, assuming proton primaries. The match between the real data and artificial data generated using the spectrum description given in  (from here on NW), and assuming proton primaries, is reassuring. We argue that these plots show that we have developed a good understanding of the techniques needed to analyse the directions and energies of very inclined showers.

The energy resolution for events above 1018 eV is shown in figure 3 for two energy bands. The energy error for the artificial events is the difference between the reconstructed energy $`E_r`$ and that of the simulated shower $`E_0`$. The rms spread in the relative errors for the artificial events are 0.43 and 0.40: we note the tail towards higher energies in the lower energy band.

The energy error in the data is obtained by combining in quadrature the error obtained from the maximum likelihood analysis and the error associated with the uncertainty in the zenith angle reconstruction. In figure 4 we show the upward and downward error distributions for the real data as compared with the artificial events, again generated with the NW spectrum and proton primaries. Again the agreement is very satisfactory.

In figure 5 the zenith angle distribution for events above 1018 eV is compared to real data using the NW (proton) spectrum to make the prediction. This comparison is similar to that made in our earlier work  but we have now developed methods to estimate the energy, most importantly including zero density measurements in the energy reconstruction fits . Here we see a difference between the predictions and the data. Note that the data are not normalised. The fact that the predicted rate is above the artificial rate indicates that we are using (i) a spectrum with fluxes that are too high, or (ii) the wrong mass composition, or (iii) an incorrect hadronic interaction model. Of course a combination of all these uncertainties can also be involved.

## 4 Rate dependence on primary mass and spectrum assumptions

We have studied how different assumptions about the shape of primary spectrum affect the expected rates of inclined showers for both pure proton and pure iron samples.

In figure 6 we show a compilation of cosmic ray data for energies above 1017 eV together with various parameterizations of the energy spectrum. The AGASA , Monocular HiRes (HiRes 1 and HiRes 2  ), HiRes-MIA , Haverah Park  and Fly’sEye Stereo  spectral measurements are as reported by the different groups. For the purposes of this work we have taken four recent parameterizations of the spectra represented by the lines in Fig. 6: the spectra attributed to Nagano and Watson , the one obtained by Szabelski et al.  (both of these are syntheses of various data sets), the parameterization given by HiRes experiment , and a parameterization based on a recent analysis of the Haverah Park data . We have not used the recently reported AGASA spectrum  which starts above 3 x 1018 eV.

In figures 7,8,9,10 we show the sensitivity of the inferred variation of mass for reasons explained later with energy to the different input spectra. While from the Szabelski et al. input we find a mass spectrum that is iron dominated over the whole energy range (1018-1019 eV), with the NW spectrum one would argue for a proton dominated composition. However if we take the recently revised Haverah Park spectrum  we find that the mean mass lies between proton and iron near 1018 eV and possibly gets heavier as we move to energies close to 1019 eV. We note that the mean mass at the differential energy point at $`2\times 10^{18}`$ eV is $`<\mathrm{ln}A>=1.4\pm 0.4`$. A similar conclusion can be drawn with the HiRes spectrum as an input. In the energy decade below 2 x 1018 eV we have recently argued that $`<\mathrm{ln}A>=2.8\pm 0.4`$, assuming a bi-modal mass composition of proton and iron .

Comparing the measured data to the corresponding predictions assuming an iron only spectrum and a proton only spectrum, it is straight forward to obtain the fraction of protons in a dual mass compositon model. In figure 11 the variation with energy of the predicted fraction of protons in the cosmic ray beam is plotted for the four spectra. In this figure, above 1019 eV, we have included events with less stringent quality cuts ($`r_{max}=2`$ km), as used in . For sufficiently high energies, the number of detectors with signals increases substantially for any surface array and this allows the relaxation of the condition, while retaining reasonably constrained fits. There are large uncertainties, but a trend towards a light composition at energies above the ankle is apparent with all spectra used. We note also that if the true flux is like that in the recent analysis presented by AGASA , the inclined shower data would be inconsistent with even a dominantly iron composition because the predicted rate would exceed the observed rate. A heavier mass primary composition would be required.

## 5 Conclusion

We believe that we have demonstrated the potential of a new method for extracting the mass spectrum above 1018 eV by studying cosmic rays at high zenith angles. It is based on the difference in abundance of muons with respect to photons and electrons for different primaries. We are not making strong claims for the finality of our conclusions (see Fig. 11) about mass composition. A definitive statement awaits a better understanding of the differences between different energy spectra above this energy and an exploration of the sensitivity of the method to different models of particle interactions. This technique can be extended to very high energies and, assuming the interaction model is understood, we would expect to be able to extract the mass spectrum in the range 1019 to 1020 eV with data from the Pierre Auger Observatory.

Acknowledgements: This work was partly supported by a joint grant from the British Council and the Spanish Ministry of Education (HB1997-0175), by Xunta de Galicia (PGIDT00PXI20615PR), by CICYT (AEN99-0589-C02-02) and by PPARC(GR/L40892).


text:
# Electrodynamics of metallic photonic crystals and problem of left-handed materials

## Abstract

An analytical theory of low frequency electromagnetic waves in metallic photonic crystals with a small volume fraction of a metal is presented. The evidence for such waves has been obtained recently by experiments and computations. The cutoff frequency of these waves, $`\omega _0`$ is studied. An analytical expression for the permittivity $`ϵ`$ is obtained and shown to be negative below $`\omega _0`$. If the crystal is embedded into a medium with a negative $`\mu `$, there are no propagating modes at any frequency. Thus, such a compound system is not a left handed material (LHM). The recent experimental results on the LHM are discussed.

In his seminal work Veselago has shown that if in some frequency range both the permittivity $`ϵ`$ and permeability $`\mu `$ are negative, the electromagnetic waves (EMW’s) propagate but they have some peculiar properties. All these properties come from the fact that vectors k, E, H form a left handed rather than a right handed set. It follows that the Poynting vector and the wave vector k have opposite directions. The materials with these properties are called the left handed materials (LHM’s).

The idea that a metallic photonic crystal (MPC) may be a technological base for the LHM  appears as a result of the computational and experimental studies of a few groups which have found the EMW’s in the MPC propagating above a very low cutoff frequency. The MPC’s they have considered are three- or two-dimensional lattices of thin straight metallic wires. Their discovery is very interesting and important because the waves propagate under the condition $`f\sigma /ϵ_0\omega 1`$, where $`f\sigma `$ is an average conductivity, $`f`$ being the volume fraction of a metal in the system. This propagation must be due to the MPC structure, because it would be impossible in a homogeneous medium with the conductivity $`f\sigma `$.

Various groups obtained different cutoff frequencies and they understood them differently. The group of Soukoulis qualitatively interpreted the effect of propagation in terms of waveguide modes, while the group of Pendry presented a completely original physical picture based upon a new longitudinal mode, called “plasma mode”. According to Pendry et al. the resulting permittivity has a plasma-like behavior

$$\frac{ϵ}{ϵ_0}=1\frac{\omega _p^2}{\omega (\omega +i\mathrm{\Gamma })},$$
(1)
however, the “plasma” frequency contains the light velocity and has a form

$$\omega _p^2=\frac{2\pi c^2}{d^2\mathrm{ln}(d/R)},\mathrm{\Gamma }=\frac{ϵ_0d^2\omega _p^2}{\pi R^2\sigma },$$
(2)
where $`d`$ is the lattice constant, $`R`$ is the radius of the metallic wires, $`\sigma `$ is the static conductivity of the metal. The same results for $`ϵ`$ and $`\omega _p`$ have been later obtained theoretically by Sarychev and Shalaev.

The San Diego group has accepted the “plasma model” and considered the negative $`ϵ`$ at $`\omega <\omega _p`$ as one of the two crucial conditions for creation of the LHM. To obtain negative $`\mu `$ the split ring resonators (SRR’s) are added to the MPC. The first observation of the negative refraction at the interface of this compound system and vacuum has been reported recently. The negative refraction is the most important manifestation of the LHM.

We claim in this paper that the “plasma mode” in the MPC is in fact an EMW and that is why the experiments of the San Diego group cannot be simply understood using $`ϵ`$ of the MPC and $`\mu `$ of the SRR’s. Our paper is organized as follows. First we discuss the arguments of Pendry et al. in favor of the plasma model and the theoretical approach by Smith et al. . Then we derive and solve an exact dispersion equation for the cutoff frequency $`\omega _0`$ and find $`ϵ(\omega )`$. The results are different from Eqs.(1,2) but $`\omega _0`$ is in a very good agreement with all computational and experimental data we are aware of. It follows that all groups discuss the same mode. The permittivity becomes negative at $`\omega <\omega _0`$. Then we show that the MPC does not support any waves if it is embedded in the medium with negative $`\mu `$ and discuss this result together with the experimental results of the San Diego group.

The Pendry group has proposed that the plasma mode in the MPC appears as a result of an electromagnetic renormalization of the electron mass in the second order in $`1/c^2`$. Due to this renormalization the electron mass becomes $`15`$ times larger than the mass of a proton. As a result, the plasma frequency of the metal shifts down into the GHz range. It is known, however, that the fields can be excluded from the interaction energy and the Lagrangian can be written as a function of instantaneous velocities and coordinates of the interacting charges keeping terms of the order of $`c^2`$. This so-called Darwin Lagrangian does not contain any mass renormalization due to the interaction. It is also strange that one should care about both mass renormalization and plasma frequency in connection with a problem based upon Maxwell equations with a static conductivity. These equations contain neither electron mass nor plasma frequency, and the very concept of electrons is not important for them.

The interpretation of the electric permittivity \[Eqs.(1,2)\] as given by Smith et al. is in terms of conventional electrodynamics and it provides a reasonable basis for a discussion. As far as we understand, it is based upon the ansatz

$$E_{ac}=\overline{E}+i\pi R^2\omega Lj,$$
(3)
where $`E_{ac}`$ is the electric field acting on a wire, $`\overline{E}`$ is the average field, $`j=\sigma E_{ac}`$ is the current density, and $`L`$ is the self-inductance of the wire per unit length. Using the ansatz (3) one can easily get Eqs.(1,2).

Now we derive an exact dispersion relation for the s-polarized EMW under the condition $`f1`$ in the system of infinite parallel thin straight wires ordered in a square lattice. The electric field of the wave is along the wires (z-axis), while the wave vector $`𝐤`$ is in the $`x`$-$`y`$ plane. Assume that the total current in each wire is $`I_0\mathrm{exp}[i(\omega t𝐤𝐫_i)]`$, where $`𝐫_i`$ is the two-dimensional radius-vector of the wire in the $`x`$-$`y`$ plane. The external solution for electric field $`E_z`$ of one wire with $`𝐫_i=0`$ has a form

$$E_z^i=(I_0\mu _0\omega /4)H_0^{(2)}(\omega \rho /c),$$
(4)
where $`H_0^{(2)}=J_0iN_0`$ is the Hankel function which decays exponentially at $`\mathrm{Im}\omega <0`$. Neither $`J_0`$ nor $`N_0`$ has this important property. Here and below we omit the time dependent factor. The solution is written in cylindrical coordinates $`z,\rho ,\varphi `$ and it obeys the boundary condition $`B_\varphi =(i/\omega )dE_z/d\rho =I_0\mu _0/2\pi \rho `$ at $`\rho =R`$.

The electric field created by all wires is

$$E_z(𝐫)=\frac{I_0\mu _0\omega }{4}e^{i𝐤𝐫}\underset{j}{}e^{i𝐤(𝐫_j𝐫)}H_0^{(2)}(\frac{\omega }{c}\rho _j),$$
(5)
where $`\rho _j=\sqrt{(xx_j)^2+(yy_j)^2}`$, summation is over all sites of the square lattice and the sum is a periodic function of $`𝐫`$.

The dispersion equation follows from the boundary condition  that relates the total electric field at the surface of any wire $`l`$ to the total current through this wire $`E_{zl}=I_0\mathrm{exp}(i𝐤𝐫_l)/\sigma _{ef}\pi R^2`$, where $`\sigma _{ef}=2\sigma J_1(\kappa R)/\kappa RJ_0(\kappa R)`$, $`\kappa =(1+i)/\delta `$, and $`\delta `$ is the skin depth. At small frequencies, when $`\delta >R`$ one gets $`\sigma _{ef}\sigma `$. At high frequencies, when $`\delta 0.1R`$, one gets the Rayleigh formula $`\sigma _{ef}(1+i)\sigma \delta /R`$.

Note, that the EMW exists mostly if the skin-effect in the wires is strong. Using Eqs.(1,2) one can show that $`\mathrm{\Gamma }/\omega _p=\delta ^2/(R^2\mathrm{ln}d/R)`$, where $`\delta `$ is the skin depth at $`\omega =\omega _p`$ (see also Ref.). We show below that the exact solution has similar properties. That is why we mostly concentrate here on the case of the strong skin-effect.

Finally, the dispersion equation for $`\omega (𝐤)(\chi i\gamma )c/d`$ has a form

$$(\chi i\gamma )\underset{l,m}{}e^{id(k_xl+k_ym)}H_0^{(2)}(z_{lm})=\frac{4cϵ_0}{d\sigma _{ef}f},$$
(6)
where $`z_{lm}=(\chi i\gamma )\sqrt{l^2+m^2+(R/d)^2}`$, $`l`$ and $`m`$ are integer numbers, the small term $`(R/d)^2`$ under the square root is important only when $`l=m=0`$. Taking real and imaginary parts of Eq.(6) one gets two equations for $`\chi `$ and $`\gamma `$.

In the continuum approximation one can substitute the summation by integration in Eq.(6) to get

$$1\frac{c^2k^2}{\omega ^2}+i\frac{f\sigma _{ef}}{ϵ_0\omega }=0.$$
(7)
Equation (7) describes propagation of a plane wave through a homogeneous medium with the conductivity $`f\sigma _{ef}`$, which is possible if $`f\sigma _{ef}/ϵ_0\omega 1`$. However, outside the continuum approximation there are propagating modes in the low frequency range $`f\sigma _{ef}/ϵ_0\omega 1`$. For these modes the fields are strongly modulated inside the lattice cell and $`𝐄_z(𝐫)`$ is close to zero near each wire so that the absorption is small.

We begin with the frequency $`\omega _0`$ which is the solution of Eq.(6) at $`|𝐤|=0`$

$$(\chi i\gamma )\underset{l,m}{}H_0^{(2)}(z_{lm})=\frac{4cϵ_0}{d\sigma _{ef}f}.$$
(8)
This mode is an eigenmode of the system and its frequency is the cutoff frequency for the EMW’s. The numerical results for the real part of the frequency are shown in Fig. 1. One can see that $`\chi `$ is of the order of few units. The values of $`\gamma `$ are of the order of the right hand side of Eq.(8). Thus, $`\chi \gamma `$ if $`f\sigma _{ef}/ϵ_0\omega 1`$.

In addition to the numerical solution we propose an approximation valid at very small $`f`$, when $`|\mathrm{ln}f|1`$. We separate the term with $`l=m=0`$ and substitute the rest of the sum by the integral. Then

$$(\chi i\gamma )\left[\frac{2i}{\pi }\left(\mathrm{ln}\frac{2}{\chi \sqrt{f/\pi }}𝐂\right)\frac{4i}{(\chi i\gamma )^2}\right]=\frac{4cϵ_0}{d\sigma _{ef}(\chi )f},$$
(9)
where we assume that $`\gamma \chi `$. Here $`𝐂`$ is the Euler’s constant. The second term at the left hand side represents the average field $`\overline{E}_z`$ which can be found from the Maxwell equation by the following way. One can show that the average (or macroscopic) magnetic induction $`\overline{𝐁}`$ is zero. Then $`\times \overline{𝐁}=0`$, $`\overline{ȷ}_z+ϵ_0\overline{E}_z/t=0`$, and

$$\overline{E}_z=\frac{I_0}{i\omega ϵ_0d^2}.$$
(10)
The second term in the square brackets of Eq.(9) is $`4\overline{E}_zd/\mu _0I_0c=4i/(\chi i\gamma )`$. Thus, the expression in the square brackets describes deviation of the field acting on a wire from the average field and it is assumed that only one term of the sum in Eq.(8) makes this difference. This approximation is similar to the ansatz of Smith et al.. The main difference between Eq.(3) and Eq.(9) which is crucial for the imaginary part of the frequency and also important for the real part, is the frequency dependence of $`\sigma _{ef}`$ due to the skin effect. Say, in the work by Smith et al. $`\delta /R710^4`$, so that the skin effect is very strong.

Figure 1 compares our result for the real part of the cutoff frequency given by Eq.(8) with the results given by Eqs.(2,9). In these calculations we assume that $`\sigma _{ef}`$ is given by the Rayleigh formula so that the right hand side of Eqs.(8,9) has a form $`\alpha (1i)\sqrt{\chi }`$, where $`\alpha =2cϵ_0R/fd\sigma \delta `$ and $`\delta `$ is taken at $`\omega =c/d`$. One can see from Fig.1 that approximation Eq.(9) is much better than approximation Eq.(2). Both approximations coincide at small $`f`$ and are accurate at extremely low values of $`f`$ ($`10^7`$) when the logarithmic term in Eqs.(2,9) is very large. The computational and experimental data of Ref. are also shown at Fig. 1 and they are in a good agreement with Eq.(8). The results of Pendry group  (not shown) are also in a good agreement with Eq.(8). Thus, we can make the conclusion that the San Diego group, group of Pendry and Soukoulis group discuss the same mode but at different values of parameters and that our analytical theory describes the same mode as well.

Now we find the component of electric permittivity $`ϵ(\omega )=ϵ_{zz}`$, which describes the s-polarized extraordinary waves in the uniaxial crystal. It is defined by the relation

$$ϵ=ϵ_0+i\frac{\stackrel{~}{\sigma }d}{\chi c},$$
(11)
where effective macroscopic conductivity $`\stackrel{~}{\sigma }`$ relates average current density $`\overline{ȷ}_z`$ to the average electric field $`\overline{E}_z`$ by equation $`\overline{ȷ}_z=\stackrel{~}{\sigma }\overline{E}_z`$. To find $`\stackrel{~}{\sigma }`$ we introduce an external electric field $`_ze^{i\omega t}`$. The average field $`\overline{E}_z`$ is given by equation

$$\overline{E}_z=_zi\frac{I_0\mu _0c}{d(\chi i\gamma )},$$
(12)
where the second term is the discussed above average field created by the wires. The boundary condition on a wire now has a form

$$\frac{I_0c(\chi i\gamma )\mu _0}{4d}\underset{l,m}{}H_0^{(2)}\left(z_{lm}\right)+_z=\frac{I_0}{\pi R^2\sigma _{ef}}.$$
(13)
Making use of Eqs.(12,13) one can find a relation between the current and the average field which gives both $`\stackrel{~}{\sigma }`$ and $`ϵ`$. Finally one gets

$$\frac{ϵ}{ϵ_0}=1\left\{\frac{\chi }{\chi i\gamma }i\chi \left[\frac{\chi i\gamma }{4}\underset{l,m}{}H_0^{(2)}\left(z_{lm}\right)\frac{cϵ_0}{fd\sigma _{eff}}\right]\right\}^1.$$
(14)
Assuming that the medium is transparent ($`\gamma 0`$, $`\sigma \mathrm{}`$) one can get the electric permittivity from the total energy density $`U`$ of the electric field $`U=(1/2)\overline{E_z^2}d(\omega ϵ)/d\omega `$. We have checked that this method gives the same result for $`ϵ(\omega )`$. The expression in the square brackets of Eq.(14) is the dispersion equation (8). One can see that $`\mathrm{Re}ϵ`$ changes sign at $`\omega =\omega _0`$ and becomes negative at $`\omega <\omega _0`$, where $`\omega _0`$ is the root of the dispersion equation (8). The derivative $`d(\mathrm{Re}ϵ)/d\omega `$ at $`\omega =\omega _0`$ is the same as the derivative which follows from Eq.(1) at $`\omega =\omega _p`$. However, far from the root of $`ϵ`$ equations (14) and (1) differ substantially.

To find $`\omega (𝐤)`$ one should solve Eq.(6). For small $`|𝐤|`$ one can get analytical result $`\omega ^2=c^2k^2+\omega _0^2`$, which is isotropic in the $`x`$-$`y`$ plane.

Now we discuss the possibility of creation of the LHM using the negative $`ϵ`$ of the MPC. Suppose that the wires are embedded into a medium with the negative magnetic permeability $`\mu `$. One can see that in this case the propagation of any EMW is suppressed. Indeed, instead of Eq.(6) one gets equation

$$\underset{j}{}e^{i(k_xx_j+k_yy_j)}K_0\left(\frac{\omega }{c_0}\sqrt{x_j^2+y_j^2+R^2}\right)=\frac{2i}{|\mu |\omega R^2\sigma _{ef}},$$
(15)
where $`c_0=1/\sqrt{|\mu |ϵ_0}`$ and $`K_0`$ is the modified Bessel function. One can see that at $`|𝐤|=0`$ all the terms on the left hand side of this equation are positive and real if $`\mathrm{Im}\omega `$ is small. Thus, if the right hand side is small, the equation cannot be satisfied. At small values of $`\omega `$ and $`|𝐤|`$ summation in the Eq.(15) can be substituted by integration. Assuming that $`\mathrm{Re}\omega \mathrm{Im}\omega `$ one gets

$$1+\frac{k^2c_0^2}{\omega ^2}=i\frac{\sigma _{ef}f}{\omega ϵ_0}.$$
(16)
This equation does not have real solutions for $`\omega (k)`$. Thus, at negative $`\mu `$ there are no propagating modes at any frequency under the study.

This result obviously follows from the fact that $`\mu ϵ_0<0`$ in the space between the wires. Therefore we get the plus sign in the left-hand side of Eq.(16) which forbids any EMW propagation. Thus, instead of the LHM we get a material without any propagating modes.

Now we compare this result with the theoretical idea to obtain the LHM, where negative $`ϵ`$ is created by the system of wires and negative $`\mu `$ is created in some other way. This idea is based upon the assumption that the negative $`ϵ`$ at $`\omega <\omega _0`$ results from a “longitudinal plasma mode”. It is taken for granted that its frequency is independent of magnetic properties of the system, which is usually the case for plasmons. However, the mode discussed above is not a plasma mode (see also ). One can show that this mode has zero average value of the magnetic induction $`\overline{𝐁}`$ over the unit cell. In this sense this is indeed a longitudinal mode. But the average value of the magnetic energy, which is proportional to $`\overline{𝐁^2}`$, is not zero and it is large. The physics of this mode is substantially related to the magnetic energy. That is why negative $`\mu `$ completely destroys this mode. It destroys also the region of negative $`ϵ`$. In fact this could be predicted from the observation that $`\omega _0c/d`$ becomes imaginary at negative $`\mu `$. Say, one can see from Eqs.(1,2) that at $`\mu <0`$ one gets $`\omega _p^2<0`$ and $`ϵ>0`$ at all frequencies assuming that $`\mathrm{\Gamma }`$ is small.

Thus, we have shown that the simple explanation  of the negative refraction in the compound system of the MPC and SRR’s, based upon the permittivity $`ϵ`$ of the MPC and the negative permeability $`\mu `$ of the SRR’s does not work because negative $`\mu `$ blocks propagation of EMW’s in the MPC. The propagation observed by the San Diego group might be a manifestation of the remarkable conclusion of Landau and Lifshitz (See Ref. p.268) that $`\mu (\omega )`$ does not have physical meaning starting with some low frequency. Then, the explanation of the negative refraction in this particular system would be outside the simple Veselago scenario (see Ref. as an example). In this case, to explain the negative refraction one should use a microscopic equation similar to Eq.(6) but with the SRR’s included.

Finally, an analytical theory of the low frequency EMW in the two-dimensional MPC is proposed. It is shown that the propagation of the low frequency waves is possible because electric and magnetic fields in the wave are strongly inhomogeneous inside the lattice cell and that electric field is small near the wires. If the dielectric part of the MPC has a negative $`\mu `$, no waves can propagate through the system. We argue that the explanation of the experiment Ref. is much deeper than it has been supposed before.

We are grateful to A. Zakhidov and V. Vardeny who have attracted our attention to this interesting problem and to L. P. Pitaevskii and B. L. Spivak for important discussions. The work has been funded by the NSF grant DMR-0102964.