text:
# Evidence for a Long-period Planet Orbiting Epsilon Eridani 1footnote 11footnote 1Based on observations collected at McDonald Observatory, Lick Observatory, European Southern Observatory at La Silla, and the Canada-France-Hawaii Telescope

## 1 Introduction

Over the past 5 years radial velocity (RV) surveys have had stunning success at finding the first giant gaseous planets ($`M`$ = 0.5 – 10 $`M_{Jupiter}`$) in orbit around other stars. Although to date over 40 of these systems have been found, none qualify as ‘true Jupiters’, i.e. giant planets in low eccentricity orbits with semimajor axes in excess of 3 AU.

Epsilon Eri ($`=`$ HR 1084) is a bright ($`V`$ = 3.7), nearby (3.22 pc) K2 V star exhibiting a high level of chromospheric activity (e.g. Gray & Baliunas 1995) that is consistent with its relatively young age of $`<`$ 1 Gyr (Soderblom & Däppen 1989). This star also has a dusty ring 60 AU from the star (Greaves et al. 1998).

Epsilon Eri has been the subject of several RV planet searches. Walker et al. (1995) using measurements spanning 11 years found evidence for $``$ 10 yr variation with an amplitude of 15 m s-1 . These results were substantiated by Nelson & Angel (1998) using an analysis of the same data set. Cumming et al. (1999) analyzed 11 years of RV data on this star taken at Lick Observatory and found significant variations with comparable amplitude but with a shorter period of 6.9 years. Because of the high level of magnetic activity for $`ϵ`$ Eri these RV variations were largely interpreted as arising from a stellar activity cycle.

The McDonald Observatory Planet Search Program (Cochran & Hatzes 1999) has been monitoring $`ϵ`$ Eri since late 1988. In this paper these results in combination with other surveys are presented. We confirm the presence of long period RV variations and demonstrate that they likely owe to the presence of a planetary companion.

## 2  The Radial Velocity Data Sets

The RV data represent measurements from four planet search groups. Data taken from late 1980 to late 1991 using the coudé spectrograph of the Canada-France-Hawaii 3.6m telescope (“CFHT” set) utilized an HF absorption cell for the velocity metric (Walker et al. 1995). Cumming et al. (1999) presented precise radial velocity measurements for $`ϵ`$ Eri taken with an iodine gas absorption cell spanning 1987.69 – 1999.99 using the Hamilton Spectrograph of the 3m Shane Telescope at Lick Observatory (“Lick” set). Observations of $`ϵ`$ Eri spanning the time interval 1992.84–1998.02 were made as part of an ESO planet search program using the 1.4m Coude Auxiliary telescope $`+`$ CES spectrograph at La Silla. Details of this program can be found in Kürster et al. (2000) and Endl et al. (2000).

The largest data set is based on observations made at the 2.7m telescope at McDonald Observatory and comprises three subsets. From 1988.74–1994.81 observations of $`ϵ`$ Eri were made using telluric O2 as the wavelength reference (“McD-O2” set). From late 1990 to late 1998 an iodine absorption cell was used as the velocity metric (“McD-I2” set). Both of these data sets utilized the “6-ft” camera of the coude spectrograph which isolated a single order of the echelle and thus had limited wavelength coverage (9 Å), but high spectral resolution (resolving power, $`R`$ = $`\lambda `$/$`\delta \lambda `$ = 240,000). Starting in 1998 the McDonald program switched to using the 2-d coude spectrograph (“McD-2dC” set) of the 2.7m telescope (Tull et al. 1995). This instrument allowed us to use the full spectral region having usable I2 lines ($``$ 5000–6100 Å) but at a lower resolving power ($`R`$ = 60,000). Analysis of the radial velocity derived from the the iodine (I2) cell technique including modeling of the instrumental profile as outlined in Valenti et al. (1995).

Table 1 summarizes the various data sets, the RV technique employed, the time interval of the observations, the number of observations, and the rms scatter, $`\sigma `$, of each set (from the fitted orbit, see $`\mathrm{\S }`$4 for details). The number of observations represent nightly averages.

## 3 Period Search

A Lomb-Scargle periodogram (Scargle 1982) was used to search for periodic signals in the data. Each data set represents relative radial velocities that are have different zero-point offsets with respect to the others. These offsets were determined in an unbiased manner from the orbital fitting (see below) which computed simultaneously all orbital elements as well as the relative velocity offsets (there is significant temporal overlap between the different sets). The top panel of Figure 1 shows the resulting periodogram of the combined data which exhibits significant power at a period of 7.13 yr ($`=`$ 2603 d). The false alarm probability (FAP) was estimated using Eq. 18 of Scargle (1982) and this yielded a FAP $`=`$ 5 $`\times `$ 10-9. This low FAP was confirmed using a bootstrap randomization scheme (e.g. Kürster et al. 1996). For comparison false data sets with the same variance as the original observations were generated by randomly re-assigning the radial velocity values the times of the observations. This was done 5$`\times 10^5`$ times and there was no instance where maximum power in the “random” periodogram exceeded that found in the data.

## 4 The Orbital Solution

An orbital solution was performed using a GaussFit model (Jefferys et al. 1988; McArthur et al. 1994). The relative velocity offsets of the different data sets were solved simultaneously with the orbital elements and a long term trend in the RV measurements. Both standard least squares and robust estimation analyses were performed.

Table 2 lists the orbital parameters from least squares and robust estimation. The solutions in Table 2 include a slight slope (0.4 m s-1 yr-1) in the velocities. The orbital solution without the trend resulted in parameters consistent to within the errors of those listed in Table 2, but with larger standard deviations. The GaussFit results were independently verified using the ORBIT program (Forveille et al. 1999) which yielded orbital parameters to within the errors of the GaussFit solution. Assuming a mass of $`M`$ $`=`$ 0.85 $`M_{}`$ for $`ϵ`$ Eri (Drake & Smith 1993) results in a $`M`$ sin $`i`$ $`=`$ 0.86 $`M_{Jupiter}`$. Figure 2 shows the robust estimation fit of the orbital solution RV measurements including the linear trend.

The rms scatter of the various data sets about the orbital curve are listed in Table 1. The average $`\sigma `$ is about 14 m s-1, excluding the O2 data. (The telluric method suffers from larger systematic errors such as winds, pressure, and temperature changes in the Earth’s atmosphere, etc.) A large part of the scatter may be related to magnetic activity. Saar & Donahue (1997) noted that stars as active as $`ϵ`$ Eri should show activity-related RV “jitter” of 12–17 m s-1. The scatter of our data is quantitatively consistent with expectations for star with the same age and activity level as $`ϵ`$ Eri.

## 5 The Nature of the RV Variations: Ca II H& K measurements

Compared to the sun, $`ϵ`$ Eri is a modestly active star and the 7 yr RV period is well within the range of periods expected for activity cycles. The Ca II H & K S-index has proved to be a powerful technique for discerning activity cycles and rotational modulation in magnetically active stars (see Baliunas et al. 1995 and references therein). If the RV variations of $`ϵ`$ Eri were due to stellar activity then one would expect to see the same period in the Ca II S-index variations which would exclude the planet hypothesis.

Observations of the relative flux in the Ca II H and K emission cores in $`ϵ`$ Eri have been made since 1966 at Mount Wilson Observatory with its 100-inch and 60-inch telescopes. Spectrophotometric records of about 100 lower main sequence stars show good long-term precision (rms $`=`$ 1-2%) and are described in a series of papers cited in Baliunas et al. (1995). Although the program of observations began with measurements made approximately monthly during the seasonal window of accessibility of each star for observations, by 1980 observations were made more frequently in order to record modulation of the Ca II flux by axial rotation owing to the uneven longitudinal distribution of Ca II emitting regions.

The interval of Ca II spectrophotometric measurements that overlap with the RV data is approximately 1980 through 1999. The Ca II data were analyzed by one of us (SLB) for sinusoidal periodicities according to the prescription for unevenly sampled data given in Horne and Baliunas (1986), and the results given in Table 3. The periodogram analysis was made on data averaged over 30-day intervals. The lower panel of Figure 1 shows the Lomb-Scargle periodogram of the 30-day S-index averages.

Two significant periodicities, one of 20 years (FAP $`=`$ 4$`\times 10^7`$) and one close to 3 years (FAP $`=`$ 6$`\times 10^6`$) are present suggesting multi-mode variations in surface magnetism, a not uncommon class of Ca II variability in younger stars like $`ϵ`$ Eri. Both periods show strong and statistically significant normalized power in the periodogram. A third, weaker period of $``$ 3.8 years (FAP $``$ 0.0001) is also possibly present.

The peak in the Ca II periodogram closest to the RV period is at 2476 d (= 6.78 yrs; $`\nu `$ = 0.147 $`\pm `$ 0.008). However, this is only the sixth highest peak in the periodogram and it has a very high false alarm probability (FAP = 16%) indicating the period’s relative insignificance. Furthermore, a detailed investigation yielded no convincing correlation between the RV and Ca II data sets. We conclude that there is no significant periodicity in the Ca II spectrophotometric measurements near that of the Doppler RV variation.

## 6 Discussion

We have provided convincing evidence for a 6.9 yr period in the radial velocity variations of $`ϵ`$ Eri based on six independent data sets taken at four telescopes using three different measurement techniques. Although these RV variations can be well fit with a Keplerian orbit, there might be some concern that RV variability may in fact be related to the high magnetic activity of this star.

A periodogram analysis of contemporaneous Ca II H&K S-index measurements revealed no significant periods at the one found in the RV analysis. This seems to exclude stellar activity as the cause of the observed RV variations. Since RV measurements provide an indirect planet detection and the exact relationship between Ca II and RV variability is poorly understood, we cannot absolutely exclude that photospheric motion related to magnetic activity is indeed responsible for the observed RV variations. However, this motion would have to be unrelated to the Ca II variability and known magnetic timescales on this star which is counter to our current understanding of stellar activity. Furthermore, it is not clear whether activity related RV variations can mimic Keplerian motion with a high orbital eccentricity over such a long time scale. Given all the available evidence the planet hypothesis is the simplest, most likely explanation for the RV variability of this star.

This extra-solar planet is particularly interesting in light of the dusty ring around $`ϵ`$ Eri. This ring is non-uniform which may be evidence for other planets in the system (Greaves et al. 1998; Liou & Zook 1999), most likely at large orbital radii ($``$ 10–30 AU). It is not clear if the planet that we have detected can dynamically influence the dust ring, although inner planets may play some role in clearing out the inner parts of the dust disk (Liou & Zook 1999). Interestingly, we do detect a linear acceleration of the central star of about 0.4 m s-1 yr-1 which is consistent with a 2 $`M_{Jupiter}`$ planet orbiting at the inner radius of the ring (30 AU). Such a planet would produce an astrometric acceleration of 0.03 mas yr-2 which can easily be confirmed by future space-based interferometry missions. Epsilon Eri may be a good candidate star for multiple planets and is deserving of future studies.

Adopting a stellar inclination of $`i`$ = 30 which is consistent with $`V`$ sin $`i`$ measurements (Saar & Osten 1997) and the ring inclination (Greaves et al. 1998) and assuming that the orbital axis of the planet, stellar rotation axis, and ring axis are all aligned, we estimate the true mass of the planet to be 1.7 $`M_{Jupiter}`$.

The candidate planet around $`ϵ`$ Eri would qualify as the first “Jupiter analog” except for its high orbital eccentricity. Giant planets in eccentric orbits is a common phenomenon among extra-solar planets (Cochran et al. 1997; Korzenik et al. 2000; Vogt et al. 2000). Explanations include interactions with a binary companion (Holman et al. 1997), interaction with the disk (Artymowicz 1992), or mergers and scatterings between two or more giant planets (e.g. Lin & Ida 1997).

It seems unlikely that the eccentricity is due to interactions with a binary companion as there is scant evidence for $`ϵ`$ Eri being a binary star. The claim of an astrometric perturbation to $`ϵ`$ Eri with a 25 yr period by van de Kamp (1974) has been refuted by Heintz (1993). Wielen et al. (1999) compared ground-based proper motion measurements of $`ϵ`$ Eri to those of Hipparcos and found that this star was a binary “on the verge of detectability and significance”. The “best-fit” orbit does have a slight linear trend which may be due to a long-period companion, but this is not compelling given that the slope amounts to a velocity change of less than 10 m s-1 over the course of our measurements or comparable to the RV rms scatter. Data covering a longer time base is needed to confirm this. At the present time there is no strong evidence to support the hypothesis that the high eccentricity is caused by interactions with a stellar companion. Given the diversity of planets showing highly eccentric orbits (now with long orbital periods) there is probably a number of mechanisms at work producing the observed eccentricities of extra-solar planets.

With a period of nearly 7 years and a $`K`$-velocity amplitude of only 19 m s-1 , the planet of $`ϵ`$ Eri has the longest orbital period and one of the lowest velocity amplitudes (for the central star) yet found. The rms scatter of the RV measurements about the orbital solution (due to measurement error and stellar activity) is comparable to the orbital $`K`$-amplitude. In spite of this the planet was detected with a high degree of confidence. This demonstrates that RV signals with amplitudes comparable to the noise level can be detected reliably given sufficient data. Ancillary measurements such as Ca II H&K emission measures are also extremely important for confirming planet discoveries, particularly long period ones.

The stellar distance of 3.2 pc yields an angular separation between star and planet of about 1 arcsecond. There is hope that we may one day be able to detect light from the planet using ground-based imaging with adaptive optics or space-based imaging. We searched the Hubble Space Telescope archives and found several images taken at 2$`\mu `$m with the coronograph of NICMOS, but no nearby companions were seen. However, the albedo of the planet at this wavelength is expected to be zero (Sudarsky et al. 2000). Furthermore, with an estimated effective temperature $``$ 100 K the radiated flux from the planet would also be undetected by NICMOS in this bandpass. To detect the thermal radiation from the planet requires observations at wavelengths near 20$`\mu `$m where the companion would have a magnitude $``$ 20. Searches for reflected light from the planet must be done in spectral regions where the planet albedo is expected to be high. For example, the theoretical albedo at 5500 Å for a gaseous giant planet with an effective temperature near 100 K is about 0.7 (Sudarsky et al. 2000). This results in a “reflected” magnitude for the planet about 21–22 magnitudes fainter than the primary star.

The astrometric perturbation of the central star is estimated to be $``$ 2 mas. This is at the limit of ground based measurements for this star (see Wielen et al. 1999) but should be measurable using space-based astrometric measurements with the Hubble Space Telescope (McArthur et al. 1999) or with future space-based or ground-based (Keck, VLTI) astrometric measurements. These combined with the RV measurements should yield orbital inclination and thus the true mass of the companion. We note that either direct imaging or astrometric measurements would provide the best confirmation of this planet.

The discovery of the planet around $`ϵ`$ Eri has now begun to push the parameter space of extra-solar planet discoveries to the long period, solar-system giant planet analogs. As RV searches lengthen their time base more of these systems will undoubtedly be discovered.

A Hatzes and W. Cochran acknowledges the support of NASA grants NAG5-4384 and NAG5-9227 and by NSF grant AST-9808980. B. McArthur receives support through NASA grant NAS5-1603. M. Endl acknowledges support by the Austrian Fond zur Förderung der wissenschaftlichen Forschung Nr. S7302. S. Baliunas is grateful to the dedication of M. A. Bradford, R. A. Donahue, K. Palmer and W. H. Soon to the HK Project at Mount Wilson Observatory, which is operated by the Mount Wilson Institute under an agreement with the Carnegie Institution of Washington. Support of this research has been made possible through the generosity of the Richard C. Lounsbery Foundation. R.P. Butler and G. Marcy thank D. Fischer, E. Williams, S. Vogt, and C. McCarthy. Finally, we thank K. Dennerl and S. Döbereiner for the help with the ESO observations.


format: text

text:
# Replica-exchange simulated tempering method for simulations of frustrated systems

## I Introduction

In complex systems such as spin glasses and biopolymers, it is very difficult to obtain accurate canonical distributions at low temperatures by conventional Monte Carlo (MC) and molecular dynamics (MD) simulation methods. This is because simulations at low temperatures tend to get trapped in one of huge number of local-minimum-energy states. One way to overcome this multiple-minima problem is to perform a simulation in a generalized ensemble where each state is weighted by a non-Boltzmann probability weight factor so that a random walk in potential energy space may be realized. (For a review of generalized-ensemble approach in the protein folding problem, see, e.g., Ref. .) The random walk allows the simulation to escape from any energy barrier and to sample much wider phase space than by conventional methods. Monitoring the energy in a single simulation run, one can obtain not only the global-minimum-energy state but also canonical ensemble averages as functions of temperature by the single-histogram  and/or multiple-histogram  reweighting techniques.

Two of the most well-known generalized-ensemble algorithms are perhaps multicanonical algorithm (MUCA)  and simulated tempering (ST)  (the latter method is also referred to as the method of expanded ensemble ). (For reviews, see, e.g., Refs. ). A simulation in MUCA performs a free 1D random walk in potential energy space. This method and its generalizations have already been used in many applications in protein systems (see, for instance, Refs. ). A simulation in ST performs a free 1D random walk in temperature space, which in turn induces a random walk in potential energy space and allows the simulation to escape from states of energy local minima again. ST has also been applied to protein folding problem .

The generalized-ensemble algorithms are powerful, but in the above two methods the probability weight factors are not a priori known and have to be determined by iterations of short trial simulations . This process can be non-trivial and very tedious for complex systems with many local-minimum-energy states.

The replica-exchange method (REM)  alleviates this difficulty. (A similar method was independently developed earlier in Ref. . REM is also referred to as multiple Markov chain method  and parallel tempering .) In this method, a number of non-interacting copies of the original system (or replicas) at different temperatures are simulated independently and simultaneously by the conventional MC or MD methods. Every few steps, pairs of replicas are exchanged with a specified transition probability. The weight factor is just the product of Boltzmann factors, and so it is essentially known. We have developed a molecular dynamics algorithm for REM  (see, also, Ref. ). We then developed a multidimensional REM .

However, REM also has a computational difficulty: As the number of degrees of freedom of the system increases, the required number of replicas also increases, and so does the required computation time. This is why we want to combine the merits of MUCA and ST and those of REM so that we can determine the weight factors for MUCA and ST with ease and save the computation time greatly. We have presented such an example; we developed a new method for the determination of the multicanonical weight factor . The method was referred to as the replica-exchange multicanonical algorithm (REMUCA). In REMUCA, a short replica-exchange simulation is performed, and the multicanonical weight factor is determined by the multiple-histogram reweighting techniques .

In this Letter we present another example of such a combination. In the new method, which we refer to as the replica-exchange simulated tempering (REST), a short replica-exchange simulation is performed, and the simulated tempering weight factor is determined by the multiple-histogram reweighting techniques . The effectiveness of the method is tested with a penta peptide, Met-enkephalin, in gas phase.

## II Methods

We first briefly review the original simulated tempering (ST) method . In this method temperature itself becomes a dynamical variable, and both the configuration and the temperature are updated during the simulation with a weight:

$$W_{\mathrm{ST}}(E;T)=e^{\beta E+a(T)},$$
(1)
where $`\beta =1/k_BT`$ ($`k_B`$ is the Boltzmann constant), $`E`$ is the potential energy, and the function $`a(T)`$ is chosen so that the probability distribution of temperature is given by

$$P_{\mathrm{ST}}(T)=𝑑En(E)W_{\mathrm{ST}}(E;T)=𝑑En(E)e^{\beta E+a(T)}=\mathrm{const},$$
(2)
where $`n(E)`$ is the density of states. Hence, in ST the temperature is sampled uniformly. A free random walk in temperature space is realized, which in turn induces a random walk in potential energy space and allows the simulation to escape from states of energy local minima.

In the numerical work we discretize the temperature in $`M`$ different values, $`T_m`$ ($`m=1,\mathrm{},M`$); we can order the temperature so that $`T_1<T_2<\mathrm{}<T_M`$. The lowest temperature $`T_1`$ should be sufficiently low so that the simulation can explore the global-minimum-energy region, and the highest temperature $`T_M`$ should be sufficiently high so that no trapping in a local-minimum-energy state occurs. The probability weight factor in Eq. (1) is then given by

$$W_{\mathrm{ST}}(E;T_m)=e^{\beta _mE+a_m},$$
(3)
where $`a_m=a(T_m)`$ ($`m=1,\mathrm{},M`$). The parameters $`a_m`$ are not a priori known and have to be determined by iterations of short simulations (see, e.g., Refs.  for details). This process can be non-trivial and very difficult for complex systems. Note that from Eqs. (2) and (3) we have

$$e^{a_m}𝑑En(E)e^{\beta _mE}.$$
(4)
The parameters $`a_m`$ are therefore “dimensionless” Helmholtz free energy at temperature $`T_m`$ (i.e., the inverse temperature $`\beta _m`$ multiplied by the Helmholtz free energy).

Once the simulated tempering parameters $`a_m`$ are determined and the initial configuration and the initial temperature $`T_m`$ are chosen, a ST simulation is realized by alternately performing the following two steps :

1. A canonical MC or MD simulation at the fixed temperature $`T_m`$ is carried out for a certain MC or MD steps.
2. The temperature $`T_m`$ is updated to the neighboring values $`T_{m\pm 1}`$ with the configuration fixed. The transition probability of this temperature-updating process is given by the Metropolis criterion (see Eq. (3)):

$$w(T_mT_{m\pm 1})=\{\begin{array}{cc}1,\hfill & \mathrm{for}\mathrm{\Delta }0,\hfill \\ \mathrm{exp}\left(\mathrm{\Delta }\right),\hfill & \mathrm{for}\mathrm{\Delta }>0,\hfill \end{array}$$
(5)
where

$$\mathrm{\Delta }=\left(\beta _{m\pm 1}\beta _m\right)E\left(a_{m\pm 1}a_m\right).$$
(6)

In the present work, we employ Monte Carlo algorithm for Step 1. After a long ST production run, one can obtain canonical ensemble averages as functions of temperature by the multiple-histogram reweighting techniques  as described in detail below.

The replica-exchange method (REM) was developed as an extension of simulated tempering  (thus it is also referred to as parallel tempering ) (see, e.g., Ref.  for a detailed description of the algorithm). The system for REM consists of $`M`$ non-interacting copies (or, replicas) of the original system in the canonical ensemble at $`M`$ different temperatures $`T_m`$ ($`m=1,\mathrm{},M`$). Let $`X=\{\mathrm{},x_m^{[i]},\mathrm{}\}`$ stand for a state in this generalized ensemble. Here, the superscript $`i`$ and the subscript $`m`$ in $`x_m^{[i]}`$ label the replica and the temperature, respectively. The state $`X`$ is specified by the $`M`$ sets of coordinates $`q^{[i]}`$ (and momenta $`p^{[i]}`$) of the atoms in replica $`i`$ at temperature $`T_m`$ ($`i,m=1,\mathrm{},M`$):

$$x_m^{[i]}\{q^{[i]},p^{[i]}\}_m.$$
(7)
In Monte Carlo algorithm we need to take into account only $`q^{[i]}`$.

A REM simulation is then realized by alternately performing the following two steps .

1. Each replica in canonical ensemble of the fixed temperature is simulated $`simultaneously`$ and $`independently`$ for a certain MC or MD steps.
2. A pair of replicas, say $`i`$ and $`j`$, which are at neighboring temperatures, $`T_m`$ and $`T_{m\pm 1}`$, respectively, are exchanged: $`X=\{\mathrm{},x_m^{[i]},\mathrm{},x_{m\pm 1}^{[j]},\mathrm{}\}X^{}=\{\mathrm{},x_m^{[j]},\mathrm{},x_{m\pm 1}^{[i]},\mathrm{}\}`$. The transition probability of this replica-exchange process is given by the Metropolis criterion:

$$w(XX^{})=\{\begin{array}{cc}1,\hfill & \mathrm{for}\mathrm{\Delta }0,\hfill \\ \mathrm{exp}\left(\mathrm{\Delta }\right),\hfill & \mathrm{for}\mathrm{\Delta }>0,\hfill \end{array}$$
(8)
where

$$\mathrm{\Delta }=\left(\beta _{m\pm 1}\beta _m\right)\left(E\left(q^{[i]}\right)E\left(q^{[j]}\right)\right).$$
(9)
Here, $`E\left(q^{[i]}\right)`$ and $`E\left(q^{[j]}\right)`$ are the potential energy of the $`i`$-th replica and the $`j`$-th replica, respectively.

In the present work we employ Monte Carlo algorithm for Step 1. A random walk in temperature space is realized for each replica, which in turn induces a random walk in potential energy space. This alleviates the problem of getting trapped in states of energy local minima.

The major advantage of REM over other generalized-ensemble methods such as MUCA  and ST  lies in the fact that the weight factor is essentially a priori known (which is just a product of Boltzmann factors), whereas in the latter algorithms the determination of the weight factors can be very tedious and time-consuming. However, the number of required replicas (or temperatures) for REM increases like $`\sqrt{N}`$ where $`N`$ is the number of degrees of freedom of the system . We will need a huge number of replicas and a large amount of computation time to simulate a complex system such as a real protein system, while in MUCA or ST only one replica (the original system itself) is simulated. This led us to combine the merits of REM and ST (the combination of REM and MUCA was also realized ).

We finally present the new method which we refer to as the replica-exchange simulated tempering (REST). In this method we first perform a short REM simulation (with $`M`$ replicas) to determine the simulated tempering weight factor and then perform with this weight factor a regular ST simulation with high statistics. The first step is accomplished by the multiple-histogram reweighting techniques  as follows. Let $`N_m(E)`$ and $`n_m`$ be respectively the potential-energy histogram and the total number of samples obtained at temperature $`T_m=1/k_B\beta _m`$ of the REM run. The density of states, $`n(E)`$, is then given by

$$n(E)=\frac{{\displaystyle \underset{m=1}{\overset{M}{}}}g_m^1N_m(E)}{{\displaystyle \underset{m=1}{\overset{M}{}}}n_mg_m^1e^{f_m\beta _mE}},$$
(10)
where

$$e^{f_m}=\underset{E}{}n(E)e^{\beta _mE}.$$
(11)
Here, $`g_m=1+2\tau _m`$, and $`\tau _m`$ is the integrated autocorrelation time at temperature $`T_m`$. Note that Eqs. (10) and (11) are solved self-consistently by iteration to obtain the density of states $`n(E)`$ and the dimensionless Helmholtz free energy $`f_m`$.

Once the estimate of the dimensionless Helmholtz free energy $`f_m`$ are obtained, the simulated tempering weight factor can be directly determined by using Eq. (3) where we set $`a_m=f_m`$ (compare Eqs. (4) and (11)). A long ST run is then performed with this weight factor. Let $`N_m(E)`$ and $`n_m`$ be respectively the potential-energy histogram and the total number of samples obtained at temperature $`T_m=1/k_B\beta _m`$ from this ST run. The multiple-histogram reweighting techniques of Eqs. (10) and (11) can be used again to obtain the best estimate of the density of states $`n(E)`$. The expectation value of a physical quantity $`A`$ at any temperature $`T(=1/k_B\beta )`$ is then calculated from

$$_T=\frac{{\displaystyle \underset{E}{}}A(E)n(E)e^{\beta E}}{{\displaystyle \underset{E}{}}n(E)e^{\beta E}}.$$
(12)
## III Results and discussion

We tested the effectiveness of the new algorithm for the system of a penta-peptide, Met-enkephalin, in gas phase. This peptide has the amino-acid sequence, Tyr-Gly-Gly-Phe-Met. The backbone was terminated by a neutral NH2\- group and a neutral -COOH group at the N-terminus and at the C-terminus, respectively, as in the previous MC works of Met-enkephalin. The potential energy function $`E_{\mathrm{tot}}`$ (in kcal/mol) is given by the sum of the electrostatic term $`E_\mathrm{C}`$, 12-6 Lennard-Jones term $`E_{\mathrm{LJ}}`$, and hydrogen-bond term $`E_{\mathrm{HB}}`$ for all pairs of atoms in the molecule together with the torsion term $`E_{\mathrm{TOR}}`$ for all torsion angles. The parameters in the energy function as well as the molecular geometry were taken from ECEPP/2 . The computer code KONF90  was used, and MC simulations based on the replica-exchange simulated tempering (REST) were performed. The dielectric constant $`ϵ`$ was set equal to 2 as in the previous works. The peptide-bond dihedral angles $`\omega `$ were fixed at the value 180 for simplicity. With this choice of parameters our results can be directly compared with Refs. .

The remaining dihedral angles constitute the variables to be updated in the MC simulations: $`\varphi _k`$ and $`\psi _k`$ in the main chain ($`k=1,\mathrm{},5`$) and $`\chi _k^j`$ in the side chains (there are 3 $`\chi `$’s for Tyr, 2 $`\chi `$’s for Phe, and 4 $`\chi `$’s for Met). For Met-enkephalin the number of degrees of freedom is thus 19. One MC sweep consists of updating all these 19 angles once with Metropolis evaluation for each update. The simulations (of all replicas) were started from randomly generated conformations.

In Table I we summarize the parameters of the simulations that were performed in the present work. As described in the previous section, REST consists of two simulations: a short REM simulation (from which the dimensionless Helmholtz free energy, or the simulated tempering weight factor, is determined) and a subsequent ST production run. The former simulation is referred to as REM1 and the latter as ST1 in Table I. In REM1 there exist 8 replicas with 8 different temperatures ($`M=8`$), ranging from 50 K to 1000 K as listed in Table I (i.e., $`T_1=50`$ K and $`T_M=T_8=1000`$ K). The same set of temperatures were also used in ST1. The temperatures were distributed exponentially between $`T_1`$ and $`T_M`$, following the optimal distribution found in the previous simulated annealing schedule , simulated tempering run , and replica-exchange simulation . Before taking the data in REM1, we made a REM simulation of 10,000 MC sweeps for thermalization. We then performed a REM simulation of $`5\times 10^4`$ MC sweeps to obtain the weight factor for simulated tempering. After estimating the weight factor, we made a ST production run of $`10^6`$ MC sweeps (ST1), which followed additional 1,000 MC sweeps for equilibration. In REM1 and ST1, a replica exchange and a temperature update, respectively, were tried every 10 MC sweeps. Data were collected at each MC sweep in both REM1 and ST1.

We first check whether the replica-exchange simulation of REM1 indeed performed properly. For an optimal performance of REM the acceptance ratios of replica exchange should be sufficiently uniform and large (say, $`>10`$ %). In Table II we list these quantities. It is clear that both points are met in the sense that they are of the same order (the values vary between 10 % and 40 %). Moreover, in Fig. 1 the canonical probability distributions obtained at the chosen 8 temperatures from REM1 are shown. We see that there are enough overlaps between all pairs of neighboring distributions, indicating that there will be sufficient numbers of replica exchange between pairs of replicas (see Table II). We did observe random walks in potential energy space between low energies and high energies.

After REM1, we obtained the dimensionless Helmholtz free energy at the eight temperatures, $`f_m`$ ($`m=1,\mathrm{},8`$), by the multiple-histogram reweighting techniques  (namely, by solving Eqs. (10) and (11) self-consistently). This gives the simulated tempering weight factor of Eq. (3) with $`a_m=f_m`$. We remark that for biomolecular systems the integrated autocorrelation times, $`\tau _m`$, in the reweighting formulae can safely be set to be a constant , and we do so throughout the analyses in this section.

After determining the simulated tempering weight factor, we carried out a long ST simulation for data collection (ST1 in Table I). In Fig. 2 the time series of temperature and potential energy from ST1 are plotted. In Fig. 2(a) we observe a random walk in the temperature space between the lowest and highest temperatures. In Fig. 2(b) the corresponding random walk of the total potential energy between low and high energies is observed. Note that there is a strong correlation between the behaviors in Figs. 2(a) and 2(b), as there should. It is known from our previous works that the global-minimum-energy conformation for Met-enkephalin has the potential energy value of $`12.2`$ kcal/mol . Hence, the random walk in Fig. 2(b) indeed visited the global-minimum region many times. It also visited high energy regions, judging from the fact that the average potential energy is around 15 kcal/mol at $`T=1000`$ K  (see Fig. 4 below).

For an optimal performance of ST, the acceptance ratios of temperature update should be sufficiently uniform and large. In Table III we list these quantities. It is clear that both points are met (the values vary between 26 % and 57 %); we find that the present ST run (ST1) indeed properly performed. We remark that the acceptance ratios in Table III are significantly larger and more uniform than those in Table II, suggesting that ST runs can sample the configurational space more effectively than REM runs, provided the optimal weight factor is obtained.

In the previous works of ST simulations of Met-enkephalin in gas phase , at least several iterations of trial simulations were required for the simulated tempering weight determination. We emphasize that in the present case of REST (REM1), only one simulation was necessary to determine the optimal simulated tempering weight factor that can cover the energy region corresponding to temperatures between 50 K and 1000 K.

In Fig. 3 we plot the simulated tempering parameters $`f_m`$ ($`m=1,\mathrm{},8`$) and the dimensionless Helmholtz free energy (inverse temperature multiplied by Helmholtz free energy) $`f(T)`$ as a function of temperature $`T`$. The former quantities, $`f_m`$, were estimated by the multiple-histogram reweighting techniques, using the results of the first REM run (REM1). The latter quantity, $`f(T)`$, was calculated by the multiple-histogram reweighting techniques, using the results of the final ST production run (ST1). Namely, the density of states, $`n(E)`$, was obtained by solving Eqs. (10) and (11) by iteration. The function $`f(T)`$ was then calculated by taking a summation of $`n(E)\mathrm{exp}(\beta E)`$ over $`E`$ for each value of $`T`$ (i.e., replace $`f_m`$ and $`\beta _m`$ in Eq. (11) by $`f(T)`$ and $`\beta `$, respectively). The results agree very well with each other, implying that the simulated tempering parameters, $`f_m`$, that were obtained from a short REM run are already quite accurate.

To check the validity of the canonical-ensemble expectation values calculated by the new method, in Fig. 4 we compare the average potential energy as a function of temperature obtained from ST1 with that obtained from the previous MUCA simulation . The results for ST1 were obtained by the multiple-histogram method  (see Eqs. (10), (11), and (12)), whereas the single-histogram method  was used for the results of the MUCA simulation. We can see a perfect coincidence of the average values in Fig. 4.

## IV Conclusions

In this Letter we have proposed a new algorithm for configurational sampling of frustrated systems with rough energy landscape. In this method, which we refer to as replica-exchange simulated tempering (REST), the simulated tempering weight factor is determined from the results of a short replica-exchange simulation, and then a regular simulated tempering production run is made with this weight. The formulation of the method is simple and straightforward, but the numerical improvement is great, because the weight factor determination for simulated tempering becomes very difficult by the usual iterative process for complex systems. The new method was tested with the system of a small peptide in gas phase. The simulated tempering weight factor was indeed obtained by a single, short replica-exchange simulation.

Acknowledgements

The authors are grateful to Dr. Y. Sugita for useful discussions. Our simulations were performed on the Hitachi and other computers at the Research Center for Computational Science, Okazaki National Research Institutes. This work is supported, in part, by a grant from the Research for the Future Program of the Japan Society for the Promotion of Science (JSPS-RFTF98P01101).

Figure Captions

* Fig. 1. Probability distribution of potential energy obtained at the eight temperatures from REM1 (see Table I for the parameters of the simulation). The left-most distribution corresponds to the lowest temperature ($`T_1=50`$ K), and the right-most to the highest one ($`T_8=1000`$ K).
* Fig. 2. Time series of temperature (a) and potential energy (b) in ST1 (see Table I for the parameters of the simulation).
* Fig. 3. Dimensionless Helmholtz free energy as a function of temperature $`T`$ (K). The dotted curve is the result from the simulated tempering production run (ST1). The crosses are the $`f_m`$ ($`m=1,\mathrm{},8`$) that were determined from the short preliminary replica-exchange run (REM1). Both results were normalized so that the values at $`T=50`$ K agree with each other.
* Fig. 4. The average potential energy as a function of temperature. The solid curve was obtained by the multiple-histogram reweighting techniques from the results of ST1. The crosses were obtained by the single-histogram reweighting techniques from the results of the previous multicanonical MC simulation .


text:
# Pion loop fluctuations of constituent quarks and baryons

## Abstract

The implications of the pion (meson) degrees of freedom for baryon properties, taken at the constituent quark and baryon levels are compared. It is shown that there is a dramatic qualitative difference between two approaches.

The role of pion degrees of freedom for baryon observables was appreciated in the strong coupling theories of the 1950’s and later formed the basis of the Skyrme model and the chiral bag models. Consideration of the pion loop fluctuations on the baryon level, either within the traditional strong coupling approach or within modern ChPT for baryons reveals several problems:

(i) The $`\pi NN`$ coupling constant is large, $`g_{\pi NN}^2/4\pi 14`$, which indicates that the loop expansion diverges. The expansion within BChPT indicates a lack of convergence as well;

(ii) The full set of intermediate baryon resonances has to be taken into account in the calculation, which consequently involves a large number of parameters;

(iii) No information about the underlying quark structure is built into the loop integrals;

(iv) $`m_N\mathrm{\Lambda }_\chi `$.

In the region of spontaneously broken chiral symmetry, where the pion degrees of freedom are important, the structure of the elementary excitations of the QCD vacuum can be approximately reproduced by absorption of the scalar interaction between bare (current) quarks in the QCD vacuum (which is responsible for the chiral symmetry breaking) into the mass of the quasiparticles, i.e. the constituent quarks. Thus in the low-energy regime the proper chiral dynamics is due to the pion-constituent quark coupling. In this case:

(i) The $`\pi QQ`$ coupling constant is relatively small, $`g_{\pi QQ}^2/4\pi 0.60.7`$;

(ii) All possible intermediate baryon resonances are taken into account automatically by considering the $`u`$ and $`d`$ quarks in the loop amplitudes;

(iii) The role of the underlying $`SU(6)`$ quark structure is revealed in relations between different observables;

(iv) $`(m_Q/\mathrm{\Lambda }_\chi )^20.1`$.

An example of the conceptual difference between the treatment of the mesonic fluctuations of constituent quarks, on the one hand, and (directly) of baryons, on the other hand, is in order. In the latter the mesonic loops imply an infinite shift of the baryon mass in the chiral limit, which is balanced by the phenomenologically determined counter-terms. Consideration of the mesonic degrees of freedom in such approaches yields information on the corrections from the finite masses of current quarks, but not on the origin of the nucleon mass nor on the octet-decuplet splitting in the chiral limit. In the chiral constituent quark model the role of the meson degrees of freedom is broader. On the one hand meson fluctuations contribute to the self-energy of the constituent quarks (loops at the quark level), but on the other hand they imply a strong flavor-dependent interaction between the constituent quarks (meson exchange interaction between different quarks), which yields an octet-decuplet splitting in the chiral limit .

The pion-exchange interaction contains both the ultraviolet (short-range) part, which is independent of the pion mass, and the infrared (Yukawa) part. The latter one is important for the long-range nuclear force, but it does not produce any significant effect in the baryon because of its small matter radius. The short-range part of the pion-exchange interaction produces a flavor-spin dependent force between quarks and has a range $`\mathrm{\Lambda }_\chi ^1`$. While the infrared (Yukawa) part of the interaction vanishes in the chiral limit, the ultraviolet part - does not. This means that in some sense the short-range part of the pion-exchange interaction is ”more fundamental” than its Yukawa part. Note that this short-range interaction stems from the $`\gamma _5`$ structure of the pion-quark vertex (i.e. exclusively from the quantum numbers of the pion) and hence is demanded by the Lorentz invariance. This short-range part of the interaction, combined with the $`SU(6)`$ structure of the zero order baryon wave functions (which is demanded by the large $`N_c`$ limit in QCD), provides a basis for the explanation of the low-lying baryon spectrum . It should be noted that this short-range interaction cannot be obtained by considering the pion loops at the baryon level with the ground $`N`$ and $`\mathrm{\Delta }`$ states as the intermediate states within the loop . At the baryon level one needs to consider the whole infinite tower of the radially excited states within the loop in order to meet this short-range meson exchange interaction. This implies that the meson-baryon or quark models which employ only the subspace of the $`\pi N`$ and $`\pi \mathrm{\Delta }`$ states are very incomplete and therefore do not take into account the most important short-range effects of the pion (meson) degrees of freedom for baryon structure.

As another example consider the different role of the pion cloud for magnetic moments in both approaches. In hadronic models the pion loops is the only source for the nucleon anomalous magnetic moment. In contrast, within the naive quark model the nucleon magnetic moments are exclusively due to the intrinsic $`SU(6)`$ quark structure of the baryon. The small pion-quark coupling constant together with the additional small parameter $`(m_Q/\mathrm{\Lambda }_\chi )^20.1`$ guarantee that the loop corrections to the naive quark model predictions (i.e. loops at the constituent quark level) are not large . For the absolute values of the proton and neutron anomalous magnetic moments one finds contributions at the level of 5-10%. What is more important, the ratio $`\mu _n/\mu _p`$, comes out to be only 2% above the empirical ratio 0.68, while the naive quark model prediction, $`2/3`$, is 2% below.

The pion loop corrections supply only a tiny contribution to the neutron Dirac charge radius . They therefore do not perturb the satisfactory description of the negative mean square radius of the neutron, which is implied by the empirical value of the neutron magnetic moment (Foldy term).

The kaon loop fluctuations of the valence $`u`$ and $`d`$ constituent quarks induce only a very small strangeness magnetic moment of the proton , which is in the range -0.01 – -0.05, and is well consistent with the most recent data from the SAMPLE experiment, reported at this conference. This is again in contrast with the large negative strange magnetic moment that is implied by the kaon loop fluctuations considered at the baryon level.


text:
# Resistance effects due to magnetic guiding orbits

## I Introduction

Recently, there has been a growing experimental and theoretical activity directed towards an increased functionality of present day electronic devices. Previously, electrical potentials were used to modify the current, while more recently one became interested in the effects of magnetic field profiles, modulated or not, on the motion of electrons in semiconductor structures. The latter is usually a heterostructure which contains a two-dimensional electron gas (2DEG). Inhomogeneous magnetic field profiles in the 2DEG are created by depositing superconducting or ferromagnetic materials on top of the heterostructure which is then patterned in the desired shape using modern nanolithography.

These hybrid systems are important from a theoretical and technological point of view, since they open the door to new physics which might result in e.g. new magneto-electronic devices. An example of such a new device is the Hybrid Hall effect device  in which the magnetic material provides a local magnetic field which influences locally the electron transport in the underlying 2DEG. The 2DEG then acts as a detector measuring the magnetic state of the magnetic material.

The fringe field arising from a magnetic stripe forms a magnetic barrier for the electron motion in the 2DEG. Barriers can be created in which the sign of the magnetic field alters in different regions of space. Due to this magnetic gradient, electrons can be bound at the boundary line between two regions of opposite magnetic field. The spectrum and the corresponding magnetic edge states have been studied recently.

When an (1D) inhomogeneous magnetic field is applied across a quasi 1D wire, these magnetic edge states are confined electrically due to the wire confinement potential and they mix with the ordinary edge states. Such a situation was recently realized by Nogaret *et al.*, where the inhomogeneous magnetic profile was arising from a perpendicularly magnetized ferromagnetic stripe grown on top of the 2DEG. They measured the magneto- and Hall resistance as function of a background magnetic field, and observed a sharp resistance resonance effect, which they attributed to the formation and subsequent killing of magnetic edge states.

In the present work we give a detailed theoretical analysis of this experiment, using a semi-classical approach in which we consider the electrical and magnetic confinement quantum mechanically, and include scattering processes using classical arguments. Both, the measured Hall resistance and the magnetoresistance will be explained. We will show that the theoretical picture of Nogaret *et al.* only captures part of the physics which is involved and is unable to predict the correct position of the peak in the magnetoresistance and the Hall resistance.

The side and top view of the experimental setup of Nogaret et al.  are shown in Fig. 1. A Hall device consisting of a $`W=2\mu `$m wide 2DEG channel in a GaAs/AlGaAs-heterojunction was fabricated, with electron density $`n_e=1.94\times 10^{15}`$m<sup>-2</sup> and mean free path $`\mathrm{}=4.5\mu `$m at 4.5K. A narrow ($`W_f=0.5\mu `$m) $`32\mu `$m long ferromagnetic (Fe or Ni) stripe (thickness $`d_f=200\mu `$m) was grown a distance $`h=80nm`$ above the center of the electron channel.

The electron transport in the 2DEG is only influenced by the perpendicular component of the magnetic stray field. In absence of any background magnetic field the ferromagnetic stripe is magnetized along the easy axis, i.e. the $`y`$-direction, and the fringe field is situated outside the quasi 1D wire, i.e. in reservoir 1 and 2. Application of a perpendicular background magnetic field rotates the magnetization to align with the $`z`$-axis, and this will result in a stray field in the wire, which imposes a step magnetic field profile along the $`y`$-direction (see Fig. 1(b)). The actual magnetic field profile is slightly rounded (see Ref.) but we checked that our results are not influenced by this simplification. This magnetic step adds an inhomogeneous magnetic field component to the uniform applied magnetic field $`B_a`$ which induces the observed resistance effects. In the present analysis, we restrict ourselves to a Fe-stripe (saturation magnetization: 1.74 T), since this was studied most thoroughly in Ref.  and produced the most pronounced resonance effect.

This paper is organized as follows. In Sec. II we present our theoretical approach. In Sect. III we calculate the two terminal resistance as function of the applied background magnetic field. The Hall resistance is studied in Sect. IV and in Sect. V the magnetoresistance is calculated. We will discuss differences between our theoretical results and the experimental (and theoretical) results of Nogaret *et al.* Our theoretical explanation for the observed resonance effect in the magnetoresistance deviates from the one proposed in Ref. . In Sect. VI we summarize our conclusions.

## II Theoretical approach

The magneto- and Hall resistance are measured experimentally by use of a four-terminal configuration. In contrast to the theoretical study of Nogaret et al., we will retain this feature in the present discussion. The four-terminal configuration is schematically shown in Fig. 2 for (a) a Hall measurement and (b) a magnetoresistance measurement. The leads are in thermodynamical equilibrium and can be characterized by a chemical potential $`\mu _i`$. Each reservoir injects a current $`I_i`$ of electrons into the 1D wire. If several bands are occupied, we have to consider a many-channel situation, and according to Büttiker,  the current in each of the leads is given by

$$I_i=\frac{e}{h}\underset{n,n^{}}{}\left\{[\delta _{n,n^{}}T_{ii}(n,n^{})]\mu _i\underset{ji}{}T_{ij}(n,n^{})\mu _j\right\},$$
(1)
where $`T_{ij}(n,n^{})`$ is the probability for an electron in channel $`n`$ of lead $`i`$ to be scattered/transmitted to $`n^{}`$ of lead $`j`$. Current conservation requires $`N_i=R_{ii}+_{ji}T_{ij}`$ for all $`i`$, with $`T_{ij}=_{n,n^{}}T_{ij}(n,n^{})`$ and $`R_{ii}=1T_{ii}`$ and $`N_i`$ is the number of channels in lead $`i`$.

Each channel $`n`$ contributes a probability $`T_{ij}(n)`$ to the conductivity which is transmitted from probe 1 to probe 2. The total transmission from probe $`i`$ to $`j`$ then equals $`T_{ij}=_{nN}T_{ij}(n)`$, and Eq. (1) is simplified to

$$I_i=\frac{e}{h}\underset{n}{}\left\{[1R_{ii}(n)]\mu _i\underset{ji}{}T_{ij}(n)\mu _j\right\}.$$
(2)
In this type of measurement, only two probes are current carrying, i.e., $`i=1,2`$, which results in the condition $`I_1=I_2=I`$ while the other probes are voltage probes and do not carry any net current: $`I_3=I_4=I_5=0`$.

In order to calculate the four-terminal magneto- and Hall resistance, we will make another simplifying assumption that the voltage probes are weakly coupled, i.e. their influence on the net current $`I`$ is very small ($`I=(\mu _1\mu _2)/R_{12,12}`$) and the chemical potentials in each of the voltage probes can be calculated in the absence of the other voltage probes ($`\mu _i=(T_{i1}\mu _1+T_{i2}\mu _2)/(T_{i1}+T_{i2})`$ with $`i=3,4`$). The general formula for this kind of resistance measurement is then readily obtained and given by

$`R_{12,3i}`$ $`=`$ $`{\displaystyle \frac{\mu _3\mu _i}{eI}}={\displaystyle \frac{h}{e^2}}{\displaystyle \frac{1}{T_{12}}}{\displaystyle \frac{T_{31}T_{i2}T_{32}T_{i1}}{(T_{31}+T_{32})(T_{i1}+T_{i2})}}`$
$`=`$ $`{\displaystyle \frac{h}{e^2}}{\displaystyle \frac{1}{T_{12}}}F=R_{12,12}F,`$
which is the two terminal resistance $`R_{12,12}`$ multiplied with a geometrical form factor $`F`$, which is less than one.

In the following we will first calculate the two-terminal resistance $`R_{12,12}`$ and then concentrate on the geometrical form factor $`F`$ in the case of a Hall or magnetoresistance measurement.

## III The energy spectrum and the two-terminal resistance

The two-terminal resistance is given by $`R_{12,12}=(\mu _2\mu _1)/I`$. We know that in the absence of any collisions, the current which flows from reservoir 1 to 2 is determined by the number of subbands $`N`$ which are occupied at the Fermi level. Since the mean free path in the experiment of Nogaret et al. is $`\mathrm{}=4\mu `$m, which is larger than the wire width ($`W=2\mu m`$), we can, to a good approximation, neglect the influence of scatterers on the spectrum and calculate the number of channels quantum mechanically following the work of Müller for a pure quasi 1D quantum wire.

We consider a system of noninteracting electrons moving in the $`xy`$-plane subjected to a hard wall confinement $`W/2<x<W/2`$, where $`W`$ is the width of the wire. The electrons are subjected to a magnetic field profile $`\stackrel{}{B}=(0,0,B_z(x))`$. This profile equals $`B_z=B_i(B_a)+B_a`$, where $`B_a`$ is the uniform applied background field and $`B_i`$ is the induced magnetic field profile due to the magnetized stripe.

In correspondence with Ref.  we will model the shape of the induced magnetic field profile by the average magnetic field on the respective sides of the magnetic stripe edges, i.e. at saturation the magnetic field profile is given by $`B_{sat}=B_1(B_1+B_2)\theta (|x|W_f/2)`$, where $`\theta `$ is the heavyside step function and $`B_1=0.06`$ Tesla and $`B_2=0.28`$ Tesla are the modeled magnetic field strengths underneath and away from the stripe as shown in Fig. 1(c). We also performed the calculations for the exact magnetic field profile, but this resulted in negligible small quantitative differences.

We model the magnetization of the stripe by considering two limiting cases: (A) when the stripe is already magnetized at $`B_a=0`$ Tesla (as was considered by Nogaret *et al*.), i.e. $`B_i=B_{sat}sign(B_a)`$ which is the hard magnet case, and (B) when the applied magnetic field magnetizes the stripe as for soft magnets. In case (B) we assume $`B_i`$ to be linearly varying with applied background magnetic field $`B_a`$, up to $`B_a=0.05`$ Tesla, where saturation is attained according to Ref. . The induced magnetic field is then given by $`B_i=B_{sat}\{1[1\theta (|B_a|0.05)](1|B_a|/0.05)\}`$. The actual experimental behaviour is expected to be situated closer to situation (B) than to (A).

The one-particle states are described by the Hamiltonian

$$H=\frac{1}{2m_e}p_x^2+\frac{1}{2m_e}\left[p_y\frac{e}{c}A(x)\right]^2+V(x),$$
(3)
where $`V(W/2<x<W/2)=0`$ and $`V(x<W/2)=V(x>W/2)=\mathrm{}`$. Taking the vector potential in the Landau gauge $`\stackrel{}{A}=(0,A_y(x),0)`$, such that $`A_y(x)/x=B_z(x)`$, we arrive at the following 2D Schrödinger equation

$$\left\{\frac{^2}{x^2}+\left[\frac{}{y}+A_y(x)\right]^2+2[EV(x)]\right\}\psi (x,y)=0,$$
(4)
where the magnetic field is expressed in $`B_0`$, magnetic units are used for a homogeneous field of $`B_0=1`$ Tesla, i.e., all lengths are measured in $`ł_0=\sqrt{\mathrm{}c/eB_0}=0.0257\mu m`$ and energy is measured in units of $`E_0=\mathrm{}eB_0/m_ec=1.7279meV`$. $`H`$ and $`p_y`$ commute due to the particular choice of the gauge, and consequently the wavefunction becomes

$$\psi (x,y)=\frac{1}{\sqrt{(2\pi )}}e^{iky}\varphi _{n,k}(x),$$
(5)
which reduces the problem to the solution of the 1D Schrödinger equation

$$\left[\frac{1}{2}\frac{d^2}{dx^2}+V_k(x)\right]\varphi _{n,k}(x)=E_{n,k}\varphi _{n,k}(x),$$
(6)
where it is the $`k`$-dependent effective potential

$$V_k(x)=\frac{1}{2}[xB_z(x)+k]^2+V(x),$$
(7)
which contains the two dimensionality of the problem. We solve Eq. (6) numerically by use of a discretization procedure.

For given applied background magnetic field we calculated the energy spectrum for case (A) with $`W=2\mu m`$. The results are shown in Fig. 3 for $`B_a/B_0=0;\mathrm{\hspace{0.17em}0.1};\mathrm{\hspace{0.17em}0.5}`$. These energy spectra are symmetric in $`k`$ and for small $`B_a`$ consist of the superposition of two parabolic spectra. For small $`k`$-values and for energies below the intersection of the two parabolas Landau levels are present due to electrons which are bound underneath the stripe. These levels shift away from each other as the background magnetic field increases, due to the increase of the magnetic field underneath the stripe ($`B_1`$). For increasing magnetic field the two parabolas shift further away from each other, towards higher $`|k|`$ values. Due to the confinement of the wire, each parabola is infinitely dublicated, where its maximum is shifted to higher energy and to lower $`k`$-values. For higher magnetic fields ($`B_a>0.5B_0`$) Landau levels arise, due to the magnetic field away from the stripe ($`B_2<B_1`$) which is now strong enough to localize electrons into cyclotron orbits.

The classical trajectories (for $`E=4E_0`$) corresponding to the different regions in $`k`$-space are shown on top of Fig. 3. We restricted ourselves to trajectories of states at energy $`E=4E_0`$, since at zero temperature only channels with energy $`E=ϵ_F=4E_0=6.9meV`$ contribute to the conductivity.

For (B) the spectrum at $`B_a/B_0=0;\mathrm{\hspace{0.17em}0.03};\mathrm{\hspace{0.17em}0.06}`$ and the corresponding classical orbits at the Fermi energy are shown in Fig. 4. For $`B_a=0`$ the magnetic stripe is not magnetized and the spectrum consists only of the potential confined levels. One single parabola (and its duplicates due to confinement) centered around $`k=0`$ is found which splits into two and its center shifts towards higher $`|k|`$-values. Below the intersection of the two shifted parabolas Landau states are formed. Notice that some levels intersect the Fermi energy twice as much as before. These Landau states separate further away from each other for increasing magnetic field. For $`B_a0.05B_0`$ the spectra are identical to the ones of (A).

The number of conducting channels is given by the number of energy levels intersecting the Fermi energy $`N`$, and hence the current is

$$I=\frac{e}{h}T_{12}(\mu _i\mu _j)$$
(8)
with $`T_{12}=_{n<N}T_{12}(n)=N`$, since in absence of any collisions $`T_{12}(n)=1`$.

Nevertheless, the mean free path measured by Nogaret et al. is smaller than the length of the wire $`L_y=16\mu m`$ and also smaller than the distance between the probes. Thus, scattering will play an important role in electron transport and consequently $`T_{12}(n)`$ will be less than 1. In order to account for this, we will estimate the transmission coefficient for every channel using classical arguments. Since we consider the voltage probes as weakly coupled, they result in a very weak perturbation of the electron-current path, and scattering due to the voltage probes will be neglected. The only scattering we consider is due to collisions with impurities and other imperfections in the 1D-channel.

The rate at which these collisions occur depends classically on the velocity in the $`y`$-direction, the length of the wire and the scattering time. The lower the velocity in the $`y`$-direction, the longer it takes to overcome the distance between probe 1 and probe 2, and the more probable it will be to experience a scattering event. Because of this we consider the transmission probability of every channel to be proportional to its velocity $`v_y=E_n(k)/k|ϵ_F`$ and the scattering time $`\tau `$, and inversely proportional to the length $`L_y`$ of the wire, i.e.,

$$T_{12}(n)\frac{v_y(n)\tau }{L_y}.$$
(9)
So finally, we arrive at the two terminal resistance

$$R_{12,12}=\frac{1}{\alpha }\frac{\mathrm{}}{e^2}\frac{1}{_nv_y(n)},$$
(10)
where $`n`$ runs over all the $`N`$ electron states with positive velocity (or negative velocity) at the Fermi energy $`ϵ_F`$, and $`\alpha `$ is a function of $`L_y`$ and $`\tau `$.

First we will discuss the change of the two-terminal resistance $`R_{12,12}`$ with respect to the situation in absence of the ferromagnetic stripe $`R_{12,12}^0`$, which we will call the *induced* resistance $`R_{12,12}/R_{12,12}^0`$. This property was also calculated and discussed by Nogaret *et al.*, and is plotted in Fig. 5 as function of the applied magnetic field $`B_a`$ for the approach of Ref.  (dashed curve) and ours (dotted and solid curves). The zero temperature result is shown in the inset.

We notice that at zero temperature, many discontinuous jumps are present. As we will see further on, their position is very sensitive to the Fermi-energy and they disappear at 4.2 K. The energy distribution function $`f(E,T)=\{\mathrm{exp}[(Eϵ_F)/k_BT]+1\}^1`$ is not a stepfunction for nonzero temperature and consequently also electrons with energy different from $`ϵ_F`$ will contribute to the conductivity. This will smoothen out these oscillations, as is shown in Fig. 5. Also the broadening of the energy levels due to e.g. potential fluctuations will have such a smoothening effect on the resistance curves. Hence, we will only show these smooth curves in the next figures.

The curves (A) and (B) differ only for $`B_a<0.05B_0`$. In case (A) the resistance for $`B_a=0`$ is larger than in the absence of the magnetic stripe. Increasing the background magnetic field results in a slight overall increase of the induced resistance. At $`B_a=0.02B_0`$ the induced resistance reaches its maximum, then it decreases rapidly.

The induced resistance in case (B) starts at 1 for $`B_a=0`$, increases more rapidly and attains its maximum at a slightly higher $`B_a`$-value, i.e., $`B_a=0.0375B_0`$. Then it decreases rapidly up to $`B_a=0.2B_0`$. We again notice oscillations at zero temperature (see inset), but fewer than for case (A). For larger $`B_a`$-values the oscillations disappear and the scaled resistance increases ultimately to one.

Nogaret *et al.* obtained theoretically a somewhat similar behaviour, as is indicated by the dashed curve, except for the peak which was situated at a slightly higher value $`B_a=0.06B_0`$. They made the assumption that the stripe was already fully magnetized at $`B_a=0`$ like in our case (A). Moreover they considered the magnetoresistance for a homogeneous magnetic field profile with magnetic field strength $`B_a`$, and considered the effect of the magnetic stripe profile by adding classical trajectories of states which arise due to the presence of the stripe. In order to simplify the problem, they only considered states which do not reach (classically) the edge of the sample. They attributed the initial positive magnetoresistance to *snake orbits* (see situation “$``$” in Fig. 3(a)) which are killed with increasing magnetic field and therefore no longer contribute to the conductivity for larger fields. At $`B_a=0.06B_0`$ all snake orbits have vanished and it is due to this, they inferred, that the resistance reaches its maximum. However if the magnetic field is larger than $`0.06B_0`$, the magnetic field has the same sign over the whole sample but has different strength under the stripe and away from it, and a new type of magnetic edge states, so called *cycloidlike* states (see the states indicated by “$``$” in Fig. 3(b)), arises, which again enhances conductivity and thus lowers the resistance. The fact that the influence of the latter orbits vanishes for larger $`B_a`$-values is due to the decrease of the velocity of these states with decreasing relative difference between the two neighboring magnetic fields.

In case (A), when the saturation magnetization is already attained at $`B_a=0`$, we cannot attribute the existence of the (small) peak to the creation or annihilation of a certain classical state. Fig. 3 shows the energy spectrum and the corresponding classical states at the Fermi-energy $`ϵ_F=4E_0`$ for (a) $`B_a=0`$, (b) $`B_a=0.1B_0`$ and (c) $`B_a=0.5B_0`$. From this figure we see that the enhancement of the resistance is a pure quantum mechanical effect and involves many different types of states with different velocities. Therefore an explanation based on the appearance or disappearance of only snake orbits as done by Nogaret *et al.* is not possible, at least for small $`B_a`$. The discontinuous behaviour for small $`B_a`$ (see the inset in Fig. 5) is due to edge states at the Fermi-level whose energy moves through the Fermi level and then no longer contribute to the conduction. They have nonzero velocity and hence this is also reflected in the resistance. For larger $`B_a>0.05B_0`$ the curve coincides with the one of case (B).

In case (B), the initial magnetoresistance can be understood more easily. For $`B_a=0`$, the stripe is not yet magnetized and thus there is no effect of the magnetic stripe. Fig. 4(a) shows subbands formed due to the quasi 1D confinement ($`N=70`$ subbands contribute to the conduction). Already for a small applied magnetic field, a relative large magnetic field is induced in the wire due to the magnetization of the ferromagnetic stripe. Whereas, for $`B_a=0`$ the only confinement was due to the edges of the sample, the magnetic field tends to localize electrons into cyclotron orbits and thus forces them in Landau levels, which separate further in energy with increasing magnetic field. As a consequence less channels will intersect the Fermi-level and consequently less channels contribute to electron transport and the resistance increases.

But there is a competing effect due to the presence of the magnetic stripe which tends to lower the induced resistance: new edge states arise (see the states indicated by “$``$” and “$``$” in Figs. 3 and 4) which travel in the opposite direction of the normal edge states. In terms of energy spectra, this effect is contained in subbands which contribute twice to the conductivity, as can be seen in Figs. 3(b,c). It overcomes the previous one at $`B_a=0.0375B_0`$, and the induced resistance decreases. This effect contributes even for $`B_a>0.06B_0`$, when the magnetic field has the same sign in the whole wire, and the previously mentioned cycloidlike orbits appear. For increasing magnetic field their influence decreases, although their number with respect to normal edge states increases. This is due to their velocity, which decreases for increasing $`B_a`$ for reasons given by Nogaret *et al.*

In the following sections, we will try to reproduce the experimental results obtained by Nogaret *et al.* First we will concentrate on the Hall resistance, before we focus on the magnetoresistance.

## IV The Hall resistance

In order to calculate the magneto- and Hall resistance we will further simplify the problem, by making the assumption of symmetrical probes. In case of the Hall resistance, the two voltage probes, i.e. probe 3 and 4, are in front of each other as is clear from Fig. 2, and due to this symmetry the transmission probabilities can be written as $`T_{31}=T_{42}=T`$ and $`T_{32}=T_{41}=t`$. The Hall resistance then attains a very simple form

$$R_{12,34}=\frac{1}{\alpha }\frac{h}{e^2}\frac{1}{T_{12}}\frac{\left(T/t1\right)}{\left(T/t+1\right)},$$
(11)
as was already derived in Ref. . Note that there is additionally one parameter $`\alpha `$ and one function $`t/T`$ which describe the behaviour of the Hall resistance. In the absence of any magnetic field $`t/T=1`$, and consequently the Hall resistance reduces to zero. As the cyclotron radius decreases, electrons will be localized closer to the edge (in edge states) and consequently the probability for an electron bouncing on one edge to be transmitted in a probe on the other side of the wire decreases drastically, i.e. exponentially, with increasing magnetic field, as can be inferred from Ref.. Consequently $`t/T`$ will decrease rapidly and ultimately for already a small applied magnetic field the geometrical form factor $`F=\left(T/t1\right)/\left(T/t+1\right)`$ will be 1 in which case the Hall resistance equals the two-terminal resistance $`R_{12,12}`$.

In order to obtain qualitative agreement with experiment, we follow Ref.  and take the following functional form for $`t/T=\mathrm{exp}[25B_a(35B_a)^2]`$ with $`B_a`$ expressed in Tesla. If we take $`R_{12,34}^0(B_a)=(3669.4B_a)\mathrm{\Omega }`$ as the functional form of the Hall resistance in absence of the magnetic stripe, which we obtained from a linear fit of the experimental result by Nogaret *et al.*, we obtain the induced Hall resistance $`R_{12,34}/R_{12,34}^0`$ as shown in Fig. 6 for case (A) and (B) (dotted and solid curves, respectively) which is compared with the experimental result (dashed curve).

We notice that the induced Hall resistance for $`B_a=0`$ approaches 0.5 in both cases (A) and (B), and increases rapidly until $`B_a=0.0325B_0`$ in case (A) and $`B_a=0.0375B_0`$ for (B). The experimental peak position of the Hall resistance $`B_a=0.04B_0`$ is very close to these values. For larger $`B_a`$ the curve almost coincides with the one in Fig. 5, which is due to the exponential form of $`t/T`$.

Notice that for case (B) the peak is very close to the experimental position and the qualitative behaviour of the Hall resistance is reproduced, the experimental curve differs quantitatively with the theoretical one only for $`B_a>0.04B_0`$. Our Hall resistance in this magnetic field range is smaller than measured experimentally. From this comparison it seems that the cycloidlike electron trajectories do not contribute much for large applied magnetic field $`B_a`$. This might be due to a large concentration of scatterers underneath the magnetic stripe edge, which might arise from the fabrication process of the sample. This would not only result in an increase of the resistance, but would especially hamper/kill the cycloidlike states propagating underneath it.

Due to the fact that the Hall resistance for large $`B_a`$ equals the two terminal resistance, it is possible to estimate $`\alpha `$ by comparison of the theoretical curve with the experimental one. In order to obtain reasonable agreement with the experimental curve, we have to assume $`\alpha =1.59`$. This value is now fixed and will be used to compare our theoretical results on the magnetoresistance with the experimental results of Ref. .

## V The magnetoresistance

In order to measure the magnetoresistance, the voltage probes are on the same side of the wire and separated a distance from each other along the 1D wire as shown schematically in Fig 2(b). If the probes are situated on the same side (which is the case for the curves under study), we can approximate the transmission probabilities by: $`T_{31}=T`$, $`T_{32}=(1\beta )t`$, $`T_{42}=(1\beta )T`$, $`T_{43}=t`$, where $`\beta <1`$ is defined as the fraction of the current $`I_{12}`$ which is reflected due to collisions between probe $`3`$ and $`4`$. The magnetoresistance is then given by

$$R_{12,34}=\frac{h}{e^2}\frac{1}{T_{12}}\frac{\left(2\beta \right)\beta }{\left[\frac{t}{T}+\frac{T}{t}\right]\left(1\beta \right)+\beta ^22\beta +2}.$$
(12)
Note that in this case there are two parameters $`\alpha `$ and $`\beta `$ and one function $`t/T`$ which determine the Hall resistance.

In absence of any magnetic field, $`t/T=1`$ and we obtain for the form factor $`F\beta /(2\beta )`$. For increasing magnetic field $`B_a`$, $`t/T`$ will decrease rapidly and for $`B_a0`$, $`t/T<1`$, which results in $`F=(t/T)(2\beta )/(1\beta )`$.

Theoretically $`\alpha `$ and $`t/T`$ are identical to those of the previous section and we have only to determine the parameter $`\beta `$. It is clear that also this parameter depends on the distance between the two voltage probes, the scattering time and the velocity of the electron states at the Fermi-level. For simplicity we will consider $`\beta (B_a)\beta `$ to be independent of the magnetic field, which is justified since the function $`t/T(B_a)`$ changes more drastically than $`\beta (B_a)`$. We will show that this approximation already results in good qualitative agreement with the experimental curves.

If we insert $`\alpha =1.59`$ and $`t/T`$ identical to the ones obtained from the Hall resistance, we arrive at the magnetoresistance shown in Fig. 7. In this figure we took $`\beta =0.95`$ and plotted $`R_{12,35}`$ for case (A) (dotted curve) and case (B) (solid curve) together with the experimental result of Nogaret *et al.* (dashed curve), which is plotted with respect to the right hand axis. We find a peak in the resistance at $`B_a=0.02B_0`$ in case (A), and $`B_a=0.0275B_0`$ for (B). The experimental peak position ($`B_a=0.03B_0`$) (dotted curve) is very close to our theoretical result for case (B). Notice that the peak position occurs for smaller $`B_a`$ then for the induced Hall resistance, which is in correspondence with the experimental results.

In contrast to the experimental results, we notice that for large $`B_a`$-values the magnetoresistance is zero. This is due to the fact that we have assumed that for large $`B_a`$-values, $`t/T=0`$. But due to scattering there is always a possibility for an electron to be scattered from an edge state localized on one side of the sample to an edge state on the other side (and traveling in the other direction). If we assume that this effect results in a constant remainder $`t_0=0.005`$ and corresponding $`t/T=(1t_0)\mathrm{exp}[25B_a(35B_a)^2]+t_0`$, we obtain the positive magnetoresistance as measured experimentally. This background does not change the Hall resistance qualitatively: the slope decreases, but this can easily be compensated with a larger $`\alpha `$ in order to have good agreement with the experiment.

It is very hard to reproduce the experimental results quantitatively, as is obvious from the need for a different left and right axis. The magnitude of the experimental result is larger, and an additional background is present. Due to the approximations made in our simple approach, we underestimated the magnetoresistance. Moreover, the experiment also suffers from other effects, like backscattering etc., which also influence the resistance but which we did not take into account in this paper. Nevertheless, we were able to reproduce the position and the magnitude ($`150\mathrm{\Omega }`$) of the peak.

## VI Conclusions

In this paper we studied electron transport in a quantum wire subjected to an abrupt magnetic field gradient arising from a ferromagnetic stripe fabricated at its surface, as was investigated experimentally by Nogaret et al. We were able to reproduce the main qualitative features of the magnetic field dependence of the Hall and magnetoresistance. In particular, the position of the peak in both resistances was correctly explained. This peak is due to two competing effects, i.e. the increase of the separation between subbands for increasing magnetic field, which decreases the number of conducting channels, and the killing of snake orbits and the creation of states which travel in the opposite direction of ordinary snake orbits, the so called cycloidlike states. Two models for the magnetization of the ferromagnetic stripe were considered corresponding to the extreme cases of a hard (A) and a soft (B) magnet. Model (B) gives the closest agreement with experiment which agrees with the observation by Nogaret *et al.* that almost no hysteresis was observed.

In comparison with the theoretical approach of Nogaret *et al.*, ours differs essentially in two ways: (1) the magnetic field profile is the one created by a soft magnet while Nogaret *et al.* assumed a magnetic barrier which is already present for zero applied magnetic field, and (2) we calculated the Hall and magnetoresistance for a four probe measurement with particular geometry, by use of a semi-classical theory based on the Landauer-Büttiker formula, while Nogaret *et al.* made use of a semi-classical drift diffusion model.

###### Acknowledgements.

This work was partially supported by the Inter-university Micro-Electronics Center (IMEC, Leuven), the Flemish Science Foundation (FWO-Vl), the “Onderzoeksraad van de Universiteit Antwerpen”, and the IUAP-IV. J.R. was supported by “het Vlaams Instituut voor de bevordering van het Wetenschappelijk & Technologisch Onderzoek in de Industrie” (IWT) and F. M. P. is a research director with the FWO-Vl. We acknowledge fruitful correspondence with A. Nogaret.