text:
# Symmetry broken motion of a periodically driven Brownian particle: nonadiabatic regime.

## Abstract

We report a theoretical study of an overdamped Brownian particle dynamics in the presence of both a spatially modulated one-dimensional periodic potential $`U(x)`$ and a periodic alternating force (AF). As the periodic potential $`U(x)`$ has a low symmetry (a ratchet potential) the Brownian particle displays a broken symmetry motion with a nonzero time average velocity. By making use of the Green function method and a mapping to the theory of Brillouin bands the probability distribution $`P(x,t)`$ of the particle coordinate $`x`$ is derived and the nonlinear dependence of the macroscopic velocity on the frequency $`\omega `$ and the amplitude $`\eta `$ of AF is found. In particular, our theory allows to go beyond the adiabatic limit ($`\omega =0`$) and to explain the peculiar reversal of the velocity sign found previously in the numerical analysis.

The model of a particle moving in a periodic one-dimensional potential $`U(x)`$ has been used to explain diverse fascinating phenomena in various fields of physics, chemistry and biology . Well known examples of such phenomena are Brillouin bands and Bloch oscillations in the solid state theory , dynamics of Josephson junctions , energy transport in various biological systems , etc.

The peculiar property of a particle motion in the periodic potential is the presence of two different states, a static state and a whirling (dynamic) state. In the dynamic state, the particle has a non zero value of the time averaged velocity $`v=<\dot{x}(t)>`$. As the system is subject to an external noise, the ”Brownian” particle motion becomes more complex (it displays both random damped oscillations and random jumps between the potential wells ), the difference between libration and rotation states disappears.

The theory of Brownian particle motion in a specific periodic potential of $`\mathrm{cos}`$-like shape and in the presence of a constant driving force (DF) has been developed in Refs. . As a particular result it was found that in the absence of DF the mean particle velocity is zero. Moreover, it is well known  that equilibrium noise only can not lead to the fluctuation induced transport, and correspondingly, $`v=0`$ for the particle motion in an arbitrary periodic potential $`U(x)`$. In this case the diffusive motion of particle occurs. The probabilities of fluctuation induced particle escape from the potential well to the left and to the right are identical. These probabilities are determined by the amplitude of potential $`U(x)`$ and do not depend on the symmetry of potential.

The situation changes drastically when a Brownian particle is subject to both a periodic spatially modulated potential (PP) and a periodic alternating force (AF). It was shown in Ref.  by numerical analysis of the corresponding Fokker-Planck equation that in this case the reflection symmetry of an overdamped particle motion ($`x`$ to $`x`$) can be broken, and directed transport occurs. The specific condition of such transport is the low symmetry of the periodic potential $`U(x)`$ (a ratchet potential). Moreover, by making use of the symmetry arguments and the numerical analysis it was shown that directed transport occurs even in the presence of symmetric PP but as the periodic AF has a low symmetry in time .

An analytical description of a Brownian particle motion in an anisotropic periodic potential has been carried out in an adiabatic regime as the frequency of AF $`\omega `$ is rather small . In this limit it was found for arbitrary values of the noise strength and the amplitude of AF $`\eta `$, the particle moves in the direction of the slower rate of potential change . However, in Ref.  the peculiar reversal of the sign of $`v`$ has been found as the AF frequency $`\omega `$ becomes relatively large. In this nonadiabatic limit and as the amplitude of AF $`\eta `$ increases, the dependence $`v(\eta )`$ displays oscillations. Moreover, these effects do not appear in the overdamped deterministic regime where the noise strength is zero . Notice here, that the directed transport of Brownian particle is more complex than the one in the deterministic case. In latter case, the non zero value of $`v`$ appears as a simple consequence of two different values of a depinning force, and the directed transport is absent as the amplitude of AF is small. Thus, in order to explain various peculiarities of the overdamped Brownian particle motion we develop a theoretical analysis that goes beyond an adiabatic approximation.

In this Letter we present a consistent analytical approach to the dynamics of an overdamped Brownian particle in the presence of both PP $`U(x)`$ with the period $`a`$ and harmonic AF, $`\eta \mathrm{cos}(\omega t)`$ . By making use of a particular diagrammatic technique that is valid in the limit of a large time, the dependence of the mean velocity on the amplitude $`\eta `$ and frequency $`\omega `$ of AF will be calculated. We find that in the presence of AF the average velocity is determined by various relaxation processes inside the potential well, and the directed transport of Brownian particle occurs as the potential $`U(x)`$ has a low symmetry. We also obtain the particular range of $`\omega `$ where the interplay between noise and AF leads to the reversal of the sign of $`v`$.

The dynamics of an overdamped particle is described by the Langeven equation  :

$$\alpha \dot{x}(t)+U^{}(x)=\eta \mathrm{cos}(\omega t)+\xi (t),$$
(1)
where the white noise function $`\xi (t)`$ has a zero mean and the correlation function $`<\xi (t)\xi (t^{})>=2\alpha T\delta (tt^{})`$. Here, $`\alpha `$ is the damping coefficient and $`T`$ is the effective temperature describing the strength of the fluctuations. Next, we introduce the time-dependent probability density $`P(x,t)`$ that satisfies the Fokker-Plank equation :

$$\alpha \frac{P(x,t)}{t}=T\frac{^2P(x,t)}{x^2}+\frac{}{x}\{[U^{}(x)\eta \mathrm{cos}(\omega t)]P(x,t)\}$$
(2)
with the initial condition

$$P(x,0)=\delta (xx_0).$$
(3)
Using the Laplace transformation $`P_\lambda (x)=_0^{\mathrm{}}𝑑tP(x,t)e^{\lambda t}`$ and the standard substitution $`P_\lambda (x,x_0)=e^{\frac{U(x)U(x_0)}{2T}}G_\lambda (x,x_0)`$ we obtain the integral equation for the Green function $`G_\lambda (x,x_0)`$

$$G_\lambda (x,x_0)=G_\lambda ^0(x,x_0)+\frac{d\lambda _1}{2\pi i}\eta \frac{2\lambda _1}{\omega ^2+\lambda _1^2}$$

$$\frac{1}{\alpha }𝑑yG_\lambda ^0(x,y)e^{\frac{U(y)}{2T}}\frac{d}{dy}[e^{\frac{U(y)}{2T}}G_{\lambda \lambda _1}(y,x_0)],$$
(4)
where $`G_\lambda ^0(x,x_0)`$ is the Green function of the equation:

$$T\frac{d^2}{dx^2}G_\lambda ^0(x,x_0)+G_\lambda ^0(x,x_0)(Te^{\frac{U(x)}{2T}}\frac{d^2}{dx^2}e^{\frac{U(x)}{2T}}\alpha \lambda )=$$

$$=\alpha \delta (xx_0).$$
(5)
The mean value of particle velocity $`v`$ is determined by the probability $`P(x,t)`$ in the limit of large time:

$$v=\underset{t\mathrm{}}{lim}𝑑xP(x,t)\frac{(xx_0)}{t}.$$
(6)
In the absence of AF ($`\eta =0`$) it immediately follows from Eq. (5) and the property  $`G_\lambda ^0(x,x_0)=G_\lambda ^0(x_0,x)`$ that the mean value of the velocity is zero for an arbitrary periodic potential $`U(x)`$. However, in the presence of AF the probability density $`P(x,t)`$ is determined by the product of the Green functions with different arguments $`\lambda `$. It leads to the breaking of the symmetry and, therefore, to the directed transport.

To obtain the fluctuation induced transport in the presence of AF we need to know the long time behaviour of the probability distribution $`P(x,t)`$ (see, the Eq. (6)). This limit corresponds to small values of $`\lambda `$ in the integral equation (4). Physically it may be interpreted as the particle return to the stable state after the multiple interactions with AF. Moreover, we are interested in a non oscillating part of $`P(x,t)`$ only. It allows greatly to simplify the Eq. (4) just keeping most singular contributions with small ”lambda” arguments of the Green functions $`G_\lambda ^0(x,y)`$. These terms can be presented in the diagrammatic form (see, Fig. 1). By summing these terms and calculating the intermediate integrals over $`\lambda _n`$ we obtain the expression for $`G_\lambda (x,x_0)`$ in the form:

$$G_\lambda (x,x_0)=𝑑pe^{ip(xx_0)}G_\lambda ^0(p;x,x_0)\frac{Z_\lambda (p)}{1Z_\lambda (p)}.$$
(7)
Here, we introduce the periodic Green function $`G_\lambda ^0(p;x,x_0)`$ written in the mixed momentum -coordinate representation, and $`Z_\lambda (p)`$ is determined by expression:

$$Z_\lambda (p)=\frac{\eta ^2}{\alpha ^2}𝑑y_1𝑑y_2P_\lambda ^0(p;y_2,y_1)$$

$$\frac{d}{dy_1}\frac{d}{dy_2}ReP_{i\omega }(y_1,y_2).$$
(8)
To calculate the function $`Z_\lambda (p)`$ and correspondingly, to obtain the conditions when the directed transport occurs, we use the general properties of the Eq. (5). This equation can be mapped to a well known problem of an electron motion in a periodic potential . Thus, the relaxation times spectrum $`\frac{\alpha }{\tau _n(p)}`$ of the Eq. (5) contains an infinite number of bands and is determined by the wave vector $`p`$. The lowest band starts from the zero value, and in the region of small wave vectors $`p`$ has a form:

$$\frac{1}{\tau _n(p)}=Dp^2,p\frac{2\pi }{a},$$
(9)
where the diffusion coefficient $`D`$ depends on the strength of the potential and temperature. In the limit of small fluctuations as the amplitude of $`U(x)`$ is large and the temperature $`T`$ is small, we obtain that the value of $`De^{\frac{U_{max}U_{min}}{T}}=e^{\frac{\mathrm{\Delta }U}{T}}`$ is also small. Note here, that the diffusion coefficient $`D`$ does not depend on the symmetry of the potential $`U`$ and is determined by the the narrow regions close to the top and bottom of potential. In the absence of AF the particle displays a diffusive motion and it is completely determined by $`D`$. The lowest band can be also mapped to a well known tight-binding model . The corresponding eigen function $`\psi _p(x)`$ has the form of a Bloch wave , namely $`\psi _p(x)=e^{ipx}r_p(x)`$, where the periodic function

$$r_p(x)\mathrm{exp}(\frac{U(x)}{2T})+pu_0(x)$$
(10)
in the limit of small $`p`$. Here, the first correction $`u_0(x)`$ in the region of small $`p`$ has been found in Refs. . Upper bands correspond to much smaller relaxation times $`\tau _n`$ but exactly these bands determine the directed transport.

In the limit of large time the main contribution to the probability distribution $`P(x,t)`$ results from a small momentum $`p`$. The presence of an odd component of $`Z_\lambda (p)p`$ leads to a nonzero value of $`v`$, and thus, by using (9), (10) and substituting (8) and (7) in (6) we arrive at the expression for $`v`$

$$v=\frac{\eta ^2}{\alpha ^2}\frac{2a}{3A}_0^a𝑑y_1_0^a𝑑y_2e^{\frac{3U(y_1)}{2T}}e^{\frac{3U(y_2)}{2T}}$$

$$\frac{d}{dy_1}Re\stackrel{~}{G}_{i\omega }(y_1,y_2),$$
(11)
where $`A=_0^a𝑑y_1_0^a𝑑y_2e^{\frac{U(y_1)}{T}}e^{\frac{U(y_2)}{T}}`$. Here, the Green function $`\stackrel{~}{G}_\lambda (x,x_0)`$ is a solution of the Eq. (5) with periodic boundary conditions.

The formula (11) is convenient for analysis of both general properties of Brownian particle motion and particular limits. Thus, a general feature of an overdamped Brownian particle transport is that the average velocity $`v`$ is determined by the properties of potential $`U(x)`$ on the whole range of its variation. Moreover, the directed transport of Brownian particle is an interplay of two effects: the fluctuation induced particle escape from the potential well (the coefficient $`1/A`$ in the Eq. (11)) and the AF induced relaxation processes inside the potential well (the Green function $`\stackrel{~}{G}_{i\omega }(y_1,y_2)`$ in the Eq. (11)). However, if the potential $`U(x)`$ displays a reflection symmetry, namely $`U(x)`$ is an even function of $`x`$, the average value of $`v`$ vanishes due to a fact that $`\frac{d}{dy_1}P_{i\omega }(y_1,y_2)`$ is an odd function of $`y_1`$ or $`y_2`$, and one of the integrals in the Eq. (11) is zero.

If the potential $`U(x)`$ has no reflection symmetry, the average velocity is not zero, and it increases as $`\eta ^2`$ in the limit of a small amplitude of AF $`\eta `$. In this limit by making use of the eigen functions $`\psi _n(x)`$ of the Eq. (5) we get:

$$<v>=\frac{\eta ^2}{\alpha ^2}\frac{2a}{3A}\underset{n}{}\frac{\tau _n}{1+(\omega \tau _n)^2}$$

$$_0^a𝑑y_1_0^a𝑑y_2e^{\frac{3U(y_1)}{2T}}e^{\frac{3U(y_2)}{2T}}\psi _n(y_2)\psi _n^{^{}}(y_1).$$
(12)
Thus, the average velocity $`v`$ is small in both limits of small and large effective temperature $`T`$:

$$v(T)\{\begin{array}{cc}\hfill e^{\frac{\mathrm{\Delta }U}{T}},& TU_{max}U_{min}\hfill \\ \hfill T^3f(\omega \tau _1),& TU_{max}U_{min}.\hfill \end{array}$$
(13)
The particular function $`f`$ depends on the ratio between the frequency $`\omega `$ and $`\tau _1^1`$ where $`\tau _1`$ is the maximum relaxation time of the particle in the potential $`U(x)`$. The Eq. (12) shows that the adiabatic regime is valid only in the limit of $`\omega 1/\tau _1`$, The relaxation time $`\tau _1`$ depends crucially on the effective temperature $`T`$. Thus, as the fluctuations are weak ($`T\mathrm{\Delta }U`$), $`\tau _1=\frac{\alpha }{\omega _0^2}`$, where $`\omega _0`$ is the frequency of small oscillations at the bottom of potential well . Moreover, $`\tau _1`$ can be even larger in the limit of small fluctuations as the potential $`U(x)`$ has a form of a double potential well . In the opposite regime of large fluctuations ($`T\mathrm{\Delta }U`$) the relaxation time depends on $`T`$ and decreases as

$$\tau _1=\frac{\alpha a^2}{4\pi ^2T}.$$
(14)
In the limit of a small frequency $`\omega `$ the main contribution to the directed transport results from the first ”relaxation time” band. It leads to a particular sign of $`v`$ that does not change with the temperature. However, as the frequency increases ($`\omega \tau _1^1`$) the upper ”relaxation time” bands dominate and the sign reversal can occur. As T increases the relaxation time $`\tau _1`$ decreases and the adiabatic limit is recovered.

This general scenario can be verified in the limit of large temperature $`T\mathrm{\Delta }U`$ for a particular model of low symmetry potential $`U(x)=U_0(\mathrm{cos}(\frac{2\pi x}{a})+\gamma \mathrm{sin}(\frac{4\pi x}{a}))`$, where the parameter $`\gamma `$ is of order one. As the temperature is large the eigen functions $`\psi _n(x)`$ of the Eq. (5) are plain waves, the main contribution to the $`v`$ results from the first and second relaxation time bands, and we obtain:

$$v(\omega )=\gamma \frac{\eta ^2}{\alpha ^2a}\frac{U_0^3}{T^3}\tau _1[\frac{1}{1+(\omega \tau _1)^2}\frac{4}{16+(\omega \tau _1)^2}],$$
(15)
where $`\tau _1`$ is determined by the Eq. (14). The dependence of $`v(\omega )`$ that displays a sign reversal, is presented in Fig. 2.

In the discussion presented to this point we assume that the amplitude of AF $`\eta `$ is small. As $`\eta `$ increases the peculiar oscillations appear in the dependence of $`v(\eta )`$. The period of these oscillations $`\mathrm{\Delta }\eta `$ can be estimated in the same limit of large temperature $`T`$. By making use of the Eq. (2) and a smallness of the potential $`U`$ we find:

$$\mathrm{\Delta }\eta (T)\{\begin{array}{cc}\hfill \alpha a\omega ,& \omega \tau _1^1\hfill \\ \hfill \alpha a\tau _1^1,& \omega \tau _1^1.\hfill \end{array}$$
(16)
In conclusion, we have presented an analysis of an overdamped Brownian particle motion in the presence of periodic potential and AF. We have shown that the directed transport in such a system is determined by two effects: the particle escape from a potential well and the relaxation processes inside the potential well. It is at variance with both the standard diffusive motion in the absence of AF and the nonlinear dynamic response of the system to a weak AF . We argue that the $`\omega `$ dependent reversal of the sign of $`v`$ found previously by numerical analysis is due to different strength and symmetry of relaxation processes. Thus, the observation of directed transport and especially its dependence on the frequency $`\omega `$ and amplitude $`\eta `$ of AF, allows to investigate the various relaxation processes of an overdamped Brownian particle. Moreover, the used method of taking into account the most singular diagrams where the momentum $`p`$ and parameter $`\lambda `$ are small, can be also useful in order to analyze the fluctuation induced transport in the presence of more complex (periodic or random) AF .

I thank S. Flach, O. Yevtushenko, Y. Zolotaryuk and A. V. Ustinov for useful discussions.


format: text

text:
# THE UNIQUE RAPID VARIABILITIES OF THE IRON K𝛼 LINE PROFILES IN NGC 4151

## 1 Introduction

Broad Fe K fluorescence line profiles, detected in many AGNs (Tanaka et al. 1995, Nandra et al. 1997, Wang et al. 1999a, Nandra et al. 1999, and references therein), are thought to arise from the accretion disk near the center massive black hole of AGNs. The line profile, thus, carries important information about the distribution and kinematics of optically material as well as about the space properties in the vicinity of the postulated super massive black hole. Doppler and gravitational shifts would imprint characteristic signatures on the line profile which map the geometric and dynamical distributions of matter surround the black hole.

Additional information concerning the geometry of X-ray source and matter in the active nucleus could in principle be derived by studying the rapid variability (on time scales of several 104 s or less) of the line profile, intensity and their relationship with the continuum variations. So far, rapid variability in Iron K line has been detected in NGC 7314 (Yaqoob et al. 1996), MCG –6-30-15 (Iwasawa et al. 1996, 1999), NGC 4051 (Wang et al. 1999b) and NGC 3516 (Nandra et al. 1999).

The observed Fe K$`\alpha `$ line rapid variabilities in different targets or during different observations of the same target are quite different. Yaqoob el al. (1996) presented the evidence for rapid variability of the Fe K line profile in the narrow-line Seyfert galaxy NGC 7314, which is consistent with a disk-line of constant equivalent width superposed on a constant flux narrow line (presumably from the torus). Iwasawa et al. (1996) discovered that, during the ASCA observation on MCG –6-30-15 in 1994, when the source is bright, the Fe K line is weak and dominated by the narrow core, whilst during a deep minimum, a huge red tail appears. The intensity of broad Fe K line correlates inversely with the continuum flux on the time scales of several 104 s. And Iwasawa et al. (1999) discovered that during the ASCA observation on the same target in 1997, the Fe K$`\alpha `$ line of a flare phase peaks around 5keV and most of its emission is shifted to below 6 keV with no component detected at 6.4keV. While for NGC 4051, during the ASCA observation in 1994, the equivalent width and the width of Fe K line correlate positively with the continuum flux (Wang et al. 1999b), which shows an opposite trend with the ASCA observation on MCG –6-30-15 in 1994. The rapid variabilities of Fe K$`\alpha `$ line profiles in NGC 3516 are also quite irregular (Nandra et al. 1999, Wang et al. 2000).

The nearby Seyfert 1.5 galaxy (z = 0.0033), NGC 4151, was first established as a bright X-ray source some thirty years ago (Gursky et al. 1971). Since that time, this archetypal low-luminosity active galaxy has invariably been considered as a prime target by all major X-ray astronomy missions. In 1995, NGC 4151 was observed by ASCA from May 10 to 12 (Leighly et al. 1997). The net exposure time of the observation is about 200 ks. Because of the enhanced sensitivity provided by ASCA and the much longer exposure time, the obtained Fe K$`\alpha `$ line profile has very high data quality (Wang et al. 1999a, here after W99). In this letter we report the discovery of unique and rapid variabilities of the iron line profiles in NGC 4151 during the ASCA observation in 1995.

## 2 Fitting the X-ray continuum

Only the SIS data are used in this paper to study the Fe K$`\alpha `$ line variabilities. Detailed description of the ASCA observation and data reduction can be found in W99. Figure 1a shows the background subtracted light curve (0.4 – 10.0 keV) with 256 s binning. Significant variability is absent on time scales of less than 1000s, but moderate variations can be seen on time scales of several $`10^4`$ s, exhibiting a continuous slow drifting behavior (e.g. Lawrence 1980). The light curve shows an intensity dip at 60 ks from the beginning, which lasts about 50 ks, and two neighboring flares at about 150 ks, for which each lasts about 30 ks. The intensity ratios of the flares to the dip are around 0.8. Apart from these, the count rate also increases slightly during the whole observation.

The time-averaged Fe K$`\alpha `$ line profile from this observation has been extracted by W99. The data quality of the line profile is one of the best that could be found in literatures (see Figure 1a in W99). In order to analysis the rapid variability properties of the Fe K$`\alpha `$ line profiles, we divide the whole observation into five time intervals (see Figure 1, I-1 to I-5). Following Weaver et al. (1994) and W99, we fit the underlying continuum of the each segment in the 1.0 – 4.0 and 8.0 – 10.0 keV bands (to exclude the broad iron line region) with a model which consists of a dual absorbed power law plus some fraction (the best-fit value for the time averaged spectrum is $``$ 5%) of the direct continuum scattered into our line of sight and absorbed only by the Galactic column of 2 $`\times `$ 1020 cm-2. We do not include a Compton reflection component in the spectral fits because of its small impact on the ASCA spectra. During the fitting, we fix the power law index and the flux of the scattered component at their best-fit values obtained from the time-averaged spectrum. Same to the time-averaged spectrum, the spectra of these five segments can also be well fitted by this model. We show the results in Table 1.

The power law indices for individual segments are all well consistent with the average value, i.e, there is no significant change in the continuum spectrum. However, the absorption columns density and the covering factor of the dual absorber showed significant variations (table 1), not correlated with variations of the continuum flux.

## 3 The Fe K$`\alpha `$ line profiles

Figure 2 shows the Fe K$`\alpha `$ line profiles for the whole observation and the five time intervals. The profiles are obtained from the ratio, the data divided by the best-fit continuum (folded through the detector response), multiplied by the power-law (original form), thus are detector effective area corrected, and also independent from the model used to fitting the lines. The time-averaged line shows a strong narrow peak around 6.4 keV and a huge red wing extending to $``$ 4.5 keV. Apart from a similar narrow core in the 6.4 keV, the line profile apparently differs by one interval to the other. The line profile of I-1 also shows a huge red wing, extending to about 4.0 keV with a possible red peak at around 4.9 keV; while for I-2, the red wing is much weaker. The line profile of I-5 is quite similar to the time averaged one. For I-3, besides the strong peak at 6.4 keV, there is also a seeming peak at around 5.4 keV. We noticed that between the two peaks there seems a valley in between which indicates that the line is double-peaked. The line profile of I-4 is similar to that of I-3, seemingly double-peaked, but much weaker. We conclude here that the Fe K$`\alpha `$ line profiles of NGC 4151 show clear and obvious variations during the ASCA observation.

Initially the line profiles are fitted with two gaussians and the results is given in Table 2. A narrow gaussian at around 6.4 keV and additional broad component can well fit the line profiles. and the results clearly present the variations of the line profiles and equivalent widths.

We also fit the Fe K line profiles with a model consist of a narrow component and a disk line (Fabian et al. 1989). The narrow core is likely emitted by material far away from the center black hole, such as torus, thus we would not expect substantial change in the line flux within the duration of the observation. During the fitting, we fix the flux of the narrow core for the five time intervals at the best-fit value from the time-averaged spectrum (which has a equivalent width of 70 eV according to the time-averaged spectrum). We also fix the disk-line energy at 6.4 keV in the source rest frame (W99), and fix the inclination angles of the disk at the best-fit value from the time-averaged spectrum. In fact, when we free the inclinations for the five time intervals, the best-fit values are all well consistent with that of the time-averaged spectrum (27 degree). As W99 did, the corresponding Fe K$`\beta `$ and Ni K$`\alpha `$ components (disk-line plus narrow core, see W99 for detail) are also included in the fit. The model can fit the line profiles quite well as the two Gaussians model do, and the results of the fit are shown in Table 3. We did not adopt the Model B in W99 because the scattered disk line component should be constant on such a short time scale, which is inconsistent with the rapid Fe K line variabilities we detected.

Nandra et al. (1999) claimed that they discovered a red shifted iron resonant absorption line at around 5.9 keV in the Fe K$`\alpha `$ profile of NGC 3516. We also add an additional absorption line component to the model, and find that including such a line can significantly improve the fit for I-4 ($`\delta \chi ^2`$ = 7). However, the absorption line is not required for other segments ($`\delta \chi ^2`$ $`<`$ 3.0).

## 4 Discussion

### 4.1 The seemingly double-peaked line profiles

Massive black holes are generally thought to exist at the centers of active galaxies (Rees 1984), but unambiguous identification of a black hole has been impeded by lack of evidence for the strong-field relativistic effects expected in the vicinity of black holes. Because of its unique anticipated profiles, the Fe K$`\alpha `$ fluorescence line has long been considered as an observational key to the study of the innermost regions of active galactic nuclei (AGN). Fabian et al. (1989) and Laor (1991) pointed out that the line profiles from the accretion disk should be broad, skewed and double-peaked (dependent on the inclination of the disk and the line-producing region in the disk), which are characteristics of the Doppler effects of the accretion disk and the strong gravitational field in the vicinity of the central black hole.

The Fe K$`\alpha `$ line profiles detected in many AGNs have been proved to be skewed and extremely broad (Tanaka et al. 1995, Nandra et al. 1997, Wang et al. 1999a, Nandra et al. 1999, and references therein), consistent with the disk line model, so give the ever strongest evidence for the presence of massive black holes and accretion disks in the center of AGNs.

However, some alternative models, such as Comptonization and jet, can also explain the broad and skewed Fe K$`\alpha `$ line profiles. Those models do not require strong gravitational field and so give no support to the presence of massive black holes (Fabian et al. 1995, Misra & Sutaria 1999). But different to the disk line model, those models predict a smooth line feature, but not a double-peaked one. Though some indirect evidences have shown that they are not satisfactory (Fabian et al. 1995, Misra 2000, Reynolds & Wilms 2000), as Misra & Kembhavi (1998) have pointed out, we infer that the detection of a double-peaked line profile is the key to confirm the disk line model for the observed Fe K$`\alpha `$ lines, and so on confirm the detection of extremely strong gravitational field.

In this paper we find that the Fe K$`\alpha `$ line profile of NGC 4151 is seemingly double-peaked during phase I-3 and I-4 in Fig1. Because an additional resonant absorption line at 5.9 keV has also been detected in the line profile in I-4, we only focus on the profile of I-3 in this section. We show the seemingly double-peaked line profile of I-3 in detail in Fig.3, and the best-fit two gaussians and disk-line models are also plotted.

The seemingly double-peaked line profile strongly suggest that the line arises from materials with disk-like geometry. Additional strong gravitational field is also needed to explain the red shift of the profile. Others models for the broad line profiles, such as Comptonization or jet, can not be responsible for a double-peaked line profile, and disk-like outflows or inflows with high velocity can not be responsible for the red shift of the line profile. The best disk line fits show that the line mainly arises from a narrow ring between 10.8 Rg and 18.6 Rg.

Currently, perhaps due to the limited energy resolution and data quality, obvious double-peaked Fe K$`\alpha `$ line profile has never been detected in any active galactic nuclei. We noticed that the line profile of MCG –6-30-15 in Fig. 1b of Wang et el. (1999a) looks double-peaked, but the plot was folded through the effective area of the detectors which rises sharply towards the low energy, so did not give the true profile of the line. While in plottings corrected for detectors’ effective area, we can not see obvious double peaks (see Fig. 7 of Iwasawa et al. 1996 and Fig. 6 of Fabian et al. 2000). We hope that the new generation X-ray satellites, such as XMM, can show us clearly double-peaked iron line profiles of AGNs.

### 4.2 The rapid variabilities of the Fe K$`\alpha `$ line profiles

Additional to the double-peaked line profiles, the rapid variabilities of the Fe K line profiles in NGC 4151 are also detected in this letter. We divided the whole ASCA observation of NGC 4151 in 1995 into five time intervals, and found that, despite the relatively small amplitude variations in the continuum flux, the Fe K$`\alpha `$ line flux and profiles show dramatic variations: during the first time interval I-1, the iron line profile is broad and skewed, showing a strong peak at 6.4 keV, and a huge red wing extending to $``$ 4.5 keV, similar to the time-averaged line profile; after several 104 s, the huge red wing becomes very weak (I-2); for I-3, except the strong peak at 6.4 keV, there is also a seeming peak at around 5.4 keV; the line profile of I-4 is similar to that of I-3, but much weaker; and finally, the line profile change back to single peak (I-5). The line equivalent widths of the five intervals are 375, 148, 329, 168 and 300 eV respectively, showing obvious variations.

The rapid variabilities of Fe K line profiles discovered previously have suggested that the X-ray continuum flux is generated by magnetic flares above the accretion disk, and at any given time there are only a few flares moving around above the disk (Iwasawa et al. 1996, Fabian 1997, Wang et al. 1999b, Iwasawa et al. 1999, Wang et al. 2000). Flares with different location and different motion can give different fluorescence line profiles.

Can the magnetic flare model also explain the observed rapid variabilities of Fe K$`\alpha `$ line in NGC 4151? For I-1 and I-5, the X-ray continuum (in the flare model context) could be dominated by a single or several strong flares located very close to the center black hole, while the inner disk radius from the disk line fits (7.4 and 8.3 Rg) gives strong supports to this model. For I-2, the X-ray continuum seems to be dominated by flares far away from the center black hole (Ri = 17 Rg). For I-3 and I-4, the X-ray continuum could be dominated by a single or several strong flares, which are $`r_g`$ above and coroating with the disk at around 14 Rg, with only the materials in the ring around 14 Rg being illuminated. The large dramatic decreasing in the line EWs from I-3 to I-4 might be due to different degree of the anisotropy in the X-ray continuum, i.e., larger fraction of the X-ray continuum received by the disk in I-3 than in I-4. However, it requires some fine tune in order to keep the observed continuum more or less constant. Alternatively, Misra (2000) suggested that the changes in the disk structure may explain the observed line profile variations, but, a disk undergoing structure changes is also hard to keep the same radiative spectrum. Finally, the decreasing of the fluorescent yielding may be due to the increasing ionization in the disk skin. Under flare context, large change in the ionization structure of disk skin is possible even for quite moderate variations in the observed continuum flux if the flares have different heights If this is case, we would expect that the flare is lower in the I-4 than in I-3. This should also reflect in the line profile. However, the uncertainty in the disk parameters do not allow us to verify this.

The observed properties strongly suggest that the rapid line variations may not be due the variations of X-ray continuum source. Instead, it may be due to changes in geometry or ionization level of the reflecting medium and at the same time, these changes should not affect the X-ray producing medium. If so, except for the variations of the X-ray continuum source, what can lead to the changes in geometry or ionization level of the reflecting medium?

An absorption line at 5.8 keV is marginally detected in I-4, but not significant in other intervals. Note a variable absorption line at the same energy has been discovered by Nandra et al. (1999) in NGC 3516. Nandra et al. (1999) interpreted this line as the redshifted Fe K resonant absorption line. The large redshift is either caused by the infalling absorbing material or by the strong gravity in the vicinity of the black hole. Given the fact that the line is highly variable and the large redshift, the agreement of the redshifts in this two objects perhaps is not purely coincidence, but have some physical origin. If so, the gravitational redshift explanation would be favored.

This work is supported by Chinese National Natural Science Foundation, Chinese Science and Technology ministry and Foundation of Ministry of Education.


text:
# Limits on the Curie temperature of (III,Mn)V ferromagnetic semiconductors

## Abstract

Mean-field-theory predicts that the Curie temperature $`T_c`$ of a (III,Mn)V ferromagnet will be proportional to the valence band density-of-states of its host (III,V) semiconductor, suggesting a route toward room-temperature ferromagnetism in this materials class. In this Letter, we use theoretical estimates of spin-wave energies and Monte-Carlo simulations to demonstrate that long-wavelength collective fluctuations, neglected by mean-field theory, will limit the critical temperature in large density-of-states materials. We discuss implications for high $`T_c`$ searches.

The recent discovery of ferromagnetism at relatively high temperatures ($`T_c>100\mathrm{K}`$) in III-V compound semiconductors containing Mn has generated intense interest, mainly because of the new technological roadways that would be opened by room temperature ferromagnetism in semiconductors with favorable materials properties. The search for systems in this materials class with higher critical temperatures is an important current activity that has been guided thus far by mean-field theoretical considerations. In this Letter we address the importance of collective magnetization fluctuations, neglected by mean-field theory, in limiting the critical-temperature and discuss the implications of these considerations for high-$`T_c`$ searches.

Our analysis is based on the kinetic-exchange model for interactions between the Mn spins and band electrons,

$$H=H_0+J_{pd}\underset{I}{}d^3r\stackrel{}{S}_I\stackrel{}{s}(\stackrel{}{r})\delta (\stackrel{}{r}\stackrel{}{R}_I),$$
(1)
where $`\stackrel{}{S}_I`$ describes a Mn spin with spin length $`S=5/2`$ at site $`\stackrel{}{R}_I`$, $`\stackrel{}{s}(\stackrel{}{r})`$ is the band-carrier density, and $`J_{pd}>0`$ represents the exchange integral. $`H_0`$ represents a simplified single-parabolic-band model for the host semiconductor valence bands.

The simplest treatment of this model is a mean field theory which takes the magnetizations of carriers and ion spins to be uniform in space and neglects correlations between them. A straightforward calculation yields for the critical temperature in mean field approximation

$$T_c^{MF}=\frac{\chi _P}{(g^{}\mu _B/2)^2}\frac{S(S+1)NJ_{pd}^2}{12},$$
(2)
where $`g^{}`$ is the g-factor of the carriers and $`\chi _P`$ is their Pauli susceptibility, which is proportional to the effective band mass. This observation has given rise to concrete predictions for critical temperatures for several host semiconductors based on their different band masses .

In mean-field theory, ferromagnetism occurs because the penalty in entropic free energy paid to polarize the Mn spins vanishes at $`T=0`$. Any coupling to a band-electron system with a finite magnetic susceptibility is sufficient to yield ferromagnetism. While this mean-field theory probably captures much of the physics of (III,Mn)V ferromagnetism, it has a qualitative deficiency which will have an important quantitative impact on $`T_c`$ predictions in circumstances we identify below. Mean-field theory fails to account for the small energy cost of magnetization configurations in which spin orientations vary slowly in space, reducing the average magnetization but maintaining local correlations between Mn and band-electrons spin orientations. In the following paragraphs we estimate the critical temperature for the case when these collective excitations dominate thermal magnetization suppression.

Isotropic ferromagnets have spin-wave Goldstone collective modes whose energies vanish at long wavelengths,

$$E(k)=Dk^2+\mathrm{},$$
(3)
where $`k`$ is the wavevector of the mode. Each spin-wave excitation reduces the total spin of the ferromagnetic state by 1. The coefficient $`D`$ is inversely proportional to the saturation magnetization and proportional to the exchange constant $`A`$ of classical micromagnetic theory, that parameterizes the free-energy cost of spatial variations in magnetization orientation. We have previously presented a theory of spin-wave excitations in (III,Mn)V ferromagnets . These collective excitations are not accounted for in the mean-field approximation. If the spin stiffness is small, they will dominate the suppression of the magnetization at all finite temperatures and limit the critical temperature. A rough upper bound  on the resulting critical temperature, $`T_c^{coll}`$, can be found by using the $`T=0`$ stiffness value and finding the temperature where the number of excited spin waves per unit volume equals the spin per volume of the ground state,

$$SN=\frac{1}{2\pi ^2}_0^{k_D}𝑑kk^2n[E(k)].$$
(4)
A Debye wavevector $`k_D=(6\pi ^2N)^{1/3}`$ cuts off the sum over wavevectors at the correct number of modes, and $`n(E)`$ is the Bose occupation number. We find that

$$k_BT_c^{coll}=\frac{2S+1}{6}Dk_D^2$$
(5)
for $`S5/2`$. To obtain this equation, we have assumed that the spin waves can be approximated as non-interacting Bose particles, replaced the dispersion by the long-wavelength limit Eq. (3), and noted that the critical temperature estimate is proportional to $`Dk_D^2`$, justifying the use of the classical expression for the mode occupation number $`n(k)k_BT/E(k)1/2`$. These considerations set an upper bound on the critical temperature which is proportional to the spin-stiffness, a bound not respected by the mean-field theory. A familar example of ferromagents in which long range order is suppressed by long-wavelength collective excitations is provided by the ferromagnetic transition metals Fe, Co, and Ni. In that case, an expression similar to Eq. (5) and proportional to the micromagnetic exchange constant, predicts critical temperatures with $`20\%`$ accuracy whereas mean-field-theory overestimates $`T_c`$ by factors of 5 to 10. In (III,Mn)V ferromagnets, we will see that both mean-field and collective regimes can occur depending on carrier density and host semiconductor band parameters.

Our theoretical results for (III,Mn)V ferromagnet collective modes lead to physically transparent results for the spin-stiffness in both strong and weak exchange coupling limits. The dimensionless parameter which characterizes the strength of the exchange coupling is the ratio $`\mathrm{\Delta }/ϵ_F`$, where $`ϵ_F`$ is the band-system Fermi energy and $`\mathrm{\Delta }=J_{pd}NS`$ is its mean-field spin splitting at $`T=0`$. For small $`\mathrm{\Delta }/ϵ_F`$, the RKKY regime, exchange coupling is a weak perturbation on the band system. In this regime we find that $`D=\delta /(12k_F^2)`$ where $`\delta =J_{pd}n\xi /2=(3/8S)(n/N)(\mathrm{\Delta }^2/ϵ_F)`$ is the energy cost of an uncorrelated spin reduction at a single Mn site. Note that in this regime $`\delta T_c^{MF}`$ and that

$$T_c^{coll,RKKY}=T_c^{MF}\frac{2S+1}{12(S+1)\sqrt[3]{2}}\left(\frac{N}{n}\right)^{2/3}.$$
(6)
In the weak coupling regime mean-field theory is reliable only for $`n/N1`$, as expected since in this case the RKKY interaction has a range which is long compared to the distance between Mn spins. In the large $`\mathrm{\Delta }/ϵ_F`$ regime exchange coupling completely polarizes the band-electron system. In this case (and for $`n2NS`$) we find that $`D=(n/2NS)(ϵ_F/k_F^2)`$. For a fully polarized band the energetic cost of varying the moment orientation direction is entirely due to band kinetic energy. We, thus, obtain as a third $`T_c`$ bound

$$T_c^{coll,s}=\frac{2S+1}{12S}ϵ_F\left(\frac{n}{N}\right)^{1/3}.$$
(7)
The different regimes deduced from these considerations are illustrated in Fig. 1.

To substantiate these qualitative considerations, we have performed Hybrid-Monte-Carlo  simulations, treating the Mn spins as discrete classical degrees of freedom, an approximation that is justified near the critical temperature. We allow for disorder by choosing the Mn positions randomly. Microscopic $`pd`$ exchange physics is modeled by allowing the interaction to have a finite range $`a_0`$ . We simulate this by replacing the delta function by a Gaussian distribution in Eq. (1).

An exhaustive description of our Monte Carlo approach, including a detailed account of all technical aspects such as thermalization procedures, finite-size effctes etc., will be given elsewhere . Here we shall, for brevity, concentrate on the results. In the following we consider the strong coupling regime where mean-field theory is not reliable and finite-size effects in our simulations are small. One important finding is that randomness in the Mn positions can enhance the spin stiffness (for large $`J_{pd}`$ by up to factor of two) in comparison to the estimates based on course-grained Mn spins with Debye cutoffs discussed above. The importance of disorder was also emphasized by Wan and Bhatt  for a model of interacting magnetic polarons where both the ion and the carrier spins are treated classically, while we account for the quantum-mechanical nature of the free carriers in our simulation.

Figure 2 shows a typical magnetization curve for manganese ions and carriers magnetizations in the collective regime. The maximum of the finite-size Mn susceptibility marks the ferromagnetic transition temperature. The lower inset compares the global and local free carrier polarizations, $`m_{loc}=\overline{|\stackrel{}{s}(\stackrel{}{r})|/n(\stackrel{}{r})}`$ . As anticipated by our qualitative discussion, the typical free carrier local band polarization remains large above the critical temperature; it would vanish in a mean-field picture. Ferromagnetism, and the technologically useful robust collective physics it gives rise to, disappears in this case only because of the loss of long-range spatial coherence.

The left panel of Fig. 3 shows the critical temperature as a function of the carrier effective mass. Mean-field theory predicts that $`T_c`$ will grow linearly with increasing mass. The Monte-Carlo results, however, are substantially lower and show a saturation of $`T_c`$ at carrier masses close to the bare electron mass. For even higher masses, we expect $`T_c`$ to decrease, reflecting the reduction of spin stiffness expected in this regime.

The right panel of Fig. 3 compares the mean-field prediction and the Monte Carlo results for $`T_c`$ as a function of the carrier density. For higher carrier densities we expect ferromagnetism to give way to spin-glass order.

We conclude from the present work that high critical temperatures cannot be achieved simply by narrowing the free carrier band or placing its Fermi energy at a density-of-states peak in order to enhance its Pauli magnetic susceptibility $`\chi _P`$. It will also be necessary to engineer a suppression of collective magnetization fluctuations.

We thank Glenn Martyna for pointing out to us the striking advantages of the Hybrid Monte Carlo algorithm, and T. Dietl, T. Jungwirth and H. Ohno for useful discussions. We acknowledge support from the Deutsche Forschungsgemeinschaft, the National Science Foundation, and the Indiana 21st Century Fund.


text:
# Electroweak phase transition beyond the Standard Model

## 1 Introduction

Even though a Higgs particle remains to be discovered, the success of precision computations gives strong support to the idea of an SU(2)$`\times `$U(1) electroweak theory with one or more Higgs doublets, some of which have, in any practical gauge, a non-vanishing expectation value at zero temperature. However, when the theory is placed in an environment with a sufficiently high temperature, the symmetry generically gets restored.  Thus we may imagine that there was an electroweak phase transition in the Early Universe.

An outstanding physics motivation for these considerations was perceived about fifteen years ago.  Indeed, an electroweak phase transition might in principle leave as a remnant the matter–antimatter (baryon) asymmetry presently observed. An essential ingredient in this statement is that there are anomalous baryon number violating interactions which are fast enough to be in thermal equilibrium at temperatures above the transition.  If the transition is sufficiently strongly of the first order, i.e. discontinuous, baryon number violating interactions are however switched off below the transition temperature, in the broken phase. Together with the CP violation known to exist in nature, as well as the non-equilibrium dynamics always present in first order transitions, these ingredients might suffice to produce an asymmetry with the correct order of magnitude (for a review, see ref. ).

In the present talk, we will not enter any details of the actual baryogenesis computations, but simply consider the existence and strength of a first order phase transition. As a reference value, we however recall that the constraint of sufficient baryon violation switch-off after the transition can be related to the Higgs expectation value over temperature $`v/T`$ in, say, the Landau gauge, or gauge invariantly to an appropriately regularized value of $`\varphi ^{}\varphi /T^2v^2/(2T^2)`$, where $`\varphi `$ is the Higgs doublet. In the case of several SU(2) Higgs doublets, the observable is correspondingly $`_i\varphi _i^{}\varphi _i/T^2`$. The numerical requirement, based on 1-loop saddle point computations (for a summary of their status and accuracy, see ) as well as on a non-perturbative evaluation,  is that the discontinuity in $`v/T`$ across the phase transition should exceed unity, $`\mathrm{\Delta }v/T>1.0`$.

## 2 General formulation of the problem

Let us start by stating more generally the types of quantities that we may want to determine. The most basic characteristic of a finite temperature system is its partition function $`Z(T,V)`$, or the free energy density $`f(T)`$, where

$$Z(T,V)=e^{\beta Vf(T)}=\mathrm{Tr}e^{\beta H}.$$
(1)
Here $`H`$ is the Hamiltonian, $`\beta =T^1`$, and $`V\mathrm{}`$ is the volume. The points $`T=T_c`$ where $`f(T)`$ is non-analytic are phase transition points, and if $`f^{}(T)|_{TT_c^\pm }`$ is discontinuous, the transition is of the first order. This implies also that the energy density, $`e(T)=T^2d[T^1f(T)]/dT`$, is discontinuous. In addition to $`f(T),e(T)`$, we also want to determine, as discussed above, the value of a (suitably regularized) $`\varphi ^{}\varphi (T)`$, which is very strongly correlated with the rate of baryon number violation in the broken phase.

In perturbation theory, one does not usually address directly $`f(T)`$ or $`\varphi ^{}\varphi `$, but rather the effective potential $`V(v,T)`$ for the length of the Higgs field, $`|\varphi |=v/\sqrt{2}`$. After gauge fixing, $`V(v,T)`$ can be defined for instance as the free energy density after the introduction of a constraint $`\delta (v/\sqrt{2}V^1d^3x|\varphi (x)|)`$ in the path integral. In the thermodynamic limit, the value of $`V(v,T)`$ at the global minimum gives the actual free energy density, $`f(T)`$.

We may now recall that for all such “static” observables, it is straightforward to write down an “imaginary time” Euclidian path integral expression. Thus, in principle we only need to compute simple well-defined functional integrals, without any ambiguities related e.g. to Minkowski space.

## 3 Why is it non-trivial?

One may ask, why not simply trust $`V(v,T)`$ computed from the functional integral using perturbation theory? Indeed, such computations have been carried out up to resummed 2-loop level for the Standard Model (SM) , as well as for the MSSM -. It turns out that at finite temperatures perturbative computations are not a priori reliable, however. We will meet circumstances below where even the order of the phase transition is wrong.

The reason for the breakdown is referred to as the infrared (IR) problem of finite temperature field theory.  It concerns “light” bosonic degrees of freedom, and can be understood as arising from expansion parameters of the type

$$g^2n_b(m)=\frac{g^2}{e^{m/T}1}\stackrel{mT}{}\frac{g^2T}{m},$$
(2)
where $`n_b(m)`$ is the Bose-Einstein distribution function, $`g^2`$ is a generic coupling constant, the SU(2) gauge coupling in particular, and $`m`$ some perturbative mass scale appearing in the propagators. Adding numerical factors from known contributions up to 3-loop order, the largest expansion parameter in the system can more concretely be estimated as

$$ϵ\frac{g^2T}{\pi m}.$$
(3)
Thus, relatively “light” ($`mg^2T`$) degrees of freedom interacting with the Higgses are a problem, and should be studied non-perturbatively. Such modes include, for instance, SU(2) gauge bosons near the symmetric phase.

In principle, a seemingly obvious way to study the system non-perturbatively is to carry out four-dimensional (4d) finite temperature lattice simulations. However, in the present context this turns out to be quite demanding. This is, first of all, because the coupling is weak, so that there are multiple length scales, $`(2\pi T)^1,(gT)^1,(g^2T)^1`$, difficult to fit simultaneously on a finite lattice. The second problem is that only the bosonic sector of the theory can be studied anyway, since chiral fermions cannot be efficiently handled in practice. Nevertheless, leaving out the fermions, large scale efforts have been undertaken, complete with an extrapolation to the continuum limit ( and references therein). The method we shall follow here is, however, different.

## 4 Dimensional reduction and 3d effective field theories

It turns out that all the problems mentioned above can be overcome in an economic and precise way by the method of finite temperature effective field theories (for reviews, see ). The basic idea of such an approach is to combine perturbation theory and simulations, in regions where they work best: one can integrate out massive modes perturbatively, which works well since the couplings are small, and then study the dynamics of the light modes non-perturbatively. In the first step, the original 4d theory reduces to a three-dimensional (3d) effective one. For studying the electroweak phase transition in a weakly coupled theory ($`m_H<250`$ GeV), this approach works with a practical accuracy at the percent level, both from the point of view of dimensional reduction,  as well as numerical simulations.

More concretely, the degrees of freedom integrated out  are all the heavy beyond-the-SM particles, as well as the non-zero Matsubara modes of the SM particles. This includes, in particular, all fermions. The effective masses of such modes are $`m>\pi T`$ so that there are no IR-problems but, according to Eq.(3), a small expansion parameter $`g^2/\pi ^2`$. The largest corrections are related to the top Yukawa coupling. After this step, the theory lives in 3d, since all non-zero Matsubara modes have been removed. A second step can also be taken:  the zero components of the gauge fields get radiatively a mass $`mgT`$. Thus there are again no IR-problems but, according to Eq. (3), a small expansion parameter $`g/\pi `$.

We will not discuss the precise details of these integrations in any more detail here. For previous reviews, we refer to , and for explicit generic rules for the integrations, to . For the logic behind the formal conjecture concerning the parametric accuracy obtained with the types of super-renormalizable effective theories discussed here (viz., $`\delta G/G𝒪(g^3)`$, where $`G`$ is a general non-vanishing bosonic Green’s function), we refer to .

## 5 Phase diagram for the Standard Model

After the heavy degrees of freedom have been integrated out perturbatively, one still has to deal with the non-perturbative bosonic scales $`m<g^2T`$. In the case of the Standard Model and many of its extensions,  the only infrared modes left are a single Higgs $`\varphi `$ and the spatial SU(2) and U(1) gauge fields, with field strength tensors $`F_{ij},B_{ij}`$. The Lagrangian is

$`_{3\mathrm{d}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{Tr}F_{ij}^2+{\displaystyle \frac{1}{4}}B_{ij}^2+(D_i\varphi )^{}D_i\varphi +m_3^2\varphi ^{}\varphi +\lambda _3(\varphi ^{}\varphi )^2.`$ (4)
All knowledge about the original modes with their physical zero temperature parameters, as well as about the temperature, is encoded in the expressions for the effective couplings $`m_3^2,\lambda _3,g_3^2,g_3^2`$. Here $`g_3,g_3^{}`$ denote the SU(2) and U(1) gauge couplings.

To get a feeling for the properties of the effective theory in Eq. (4), let us first apply 1-loop perturbation theory. Ignoring the tiny corrections from $`g_3^2/g_3^2`$ for the moment and tuning $`m_3^2`$ to be on the phase transition line, we find a first order transition, with the discontinuity

$$\mathrm{\Delta }\frac{v}{gT}=\frac{1}{8\pi }\frac{g_3^2}{\lambda _3}\left[1+𝒪(\frac{\lambda _3}{g_3^2})\right].$$
(5)
(The $`𝒪(\lambda _3/g_3^2)`$ corrections are in fact not analytic.) Using that in the SM,

$$\frac{\lambda _3}{g_3^2}\frac{1}{8}\frac{m_H^2}{m_W^2}+𝒪\left(\frac{g^2}{(4\pi )^2}\frac{m_{\text{top}}^4}{m_W^4}\right),$$
(6)
and inspecting, say, $`m_W(T)=gv/2`$ in the broken phase, we however see from Eq. (3) that for large (realistic) Higgs masses $`m_H>m_W`$ an IR sensitive computation of the type in Eq. (5) cannot be trusted at all.

Fortunately, we do not need to rely on perturbation theory. The theory in Eq. (4) is ideally suited for lattice simulations. Due to the simplicity of the action and its low dimensionality, extrapolations to the infinite volume and continuum limits can be systematically carried out .

The results of such lattice simulations from refs.  are summarised in Fig. 1. There is a line of first order phase transitions, which however ends at a critical point, after which there is only a crossover. The ending of the transition is in qualitative contrast with the perturbative prediction in Eq. (5).

The endpoint location has been determined with great precision, and corresponds in physical units to $`m_{H,c}=72.3(7)`$ GeV, $`T_c=109.2(8)`$ GeV . The errors here are statistical, remaining after continuum extrapolation. In case one would now not be a priori convinced about the accuracy of dimensional reduction at the percent level, we may note that a similar phase diagram exists also in the bosonic 4d SU(2)$`+`$Higgs theory, where also direct 4d simulations have been carried out.  While the statistical errors are slightly larger, the results agree completely with the results obtained after leaving out fermions from the dimensional reduction expressions, as well as changing $`g^2`$ to the larger value used in ref. . A detailed comparison of the whole phase diagram has been carried out in ref. , with full agreement for $`m_H>30`$ GeV. Thus, we can really rely on the result in Fig. 1.

Let us now discuss the physical conclusions drawn from Fig. 1. As pointed out in the Introduction, for baryogenesis we need a first order phase transition, and even a strong one, $`\mathrm{\Delta }v/T>1`$. We have learned that a first order transition only exists for $`m_H<72`$ GeV. Furthermore, it has in fact $`\mathrm{\Delta }v/T<1.0`$ down to $`m_H10`$ GeV!  Thus, we observe that experimentally allowed Higgs masses $`m_H>110`$ GeV  are very far from allowing for a first order electroweak phase transition, let alone for electroweak baryogenesis.

Finally, let us briefly mention a related physics topic, which is of some interest because its solution is not yet understood as well as the issues above. It is the question as to how the properties of the phase transition change if there are primordial magnetic fields present,  or, more generally, whether magnetic fields could in some way be related to the baryon asymmetry.

## 6 How to change the properties of the phase transition?

From the discussion in the previous section, it is clear that if we want to have an electroweak phase transition which is of the first order for realistic Higgs masses, and even strongly so, then a drastic change is needed with respect to the Standard Model.

What is it that should be changed? As the strength of the transition is determined by the scalar self-coupling, cf. Eq. (5), we apparently need some new degree of freedom which can decrease the effective $`\lambda _3`$ by $`𝒪(100\%)`$. To get such a large correction in the first place, we need a bosonic zero mode with an expansion parameter like in Eq. (3). A simple perturbative 1-loop computation shows that at least the sign of loop corrections of this type is the correct one, negative. But to have an effect of $`𝒪(100\%)`$, we need $`ϵ1`$, so that $`mg^2T/\pi `$. That is, the degree of freedom should itself be non-perturbative!

In principle, there are two kinds of possibilities for such degrees of freedom. Either it is a new type of a gauge field, maybe something like technicolour, or a scalar. In the latter, more natural case, there are thermal corrections in the effective mass parameter, appearing as $`m_3^2m_{\text{4d}}^2+g^2T^2`$. Thus, to get a total outcome of order $`(g^2T)^2`$, a cancellation must take place such that the mass parameter $`m_{\text{4d}}^2`$ is negative. At zero temperature, the physical mass is then roughly $`m_{\text{phys}}^2<m_{\text{4d}}^2+m_{\text{top}}^2<m_{\text{top}}^2`$. In order for such a relatively light degree of freedom not to have shown up so far in precision electroweak observables, it should be an SU(2) singlet.

It is interesting that scalar degrees of freedom satisfying these requirements can be found in realistic theories. Consider, for instance, the MSSM. Potential candidates are squarks and sleptons. Since they should couple relatively strongly to the Higgs, we should choose stops.

Now, the stops come left- and right-handed, $`\stackrel{~}{t}_L,\stackrel{~}{t}_R`$ (these states can also mix, but for lack of space we ignore this here). The requirement of having a small violation of the electroweak precision observables means that the weakly interacting left-handed one cannot be light, $`m_{\stackrel{~}{t}_L}m_{\text{top}}`$. This helps also in getting a large Higgs mass (see, e.g., ref. ),

$$m_H^2m_Z^2\mathrm{cos}^22\beta +\frac{3g^2}{8\pi ^2}\frac{m_{\text{top}}^4}{m_W^2}\mathrm{ln}\frac{m_{\stackrel{~}{t}_R}m_{\stackrel{~}{t}_L}}{m_{\text{top}}^2}.$$
(7)
On the contrary, the SU(2) singlet stop $`\stackrel{~}{t}_R`$ can be “light”, and serve as the desired new degree of freedom. (Then we of course lose half of the correction to $`m_H`$ in Eq. (7), but this is the price to pay.) It should not be so light that the stop direction gets broken before the electroweak phase transition, though, because then one cannot get back to the SM minimum afterwards.

Another example of a viable light scalar degree of freedom is the complex gauge singlet of the NMSSM, which can also get broken during the transition (for a recent study, see ).

## 7 Numerical results for the MSSM

Under the conditions described in the previous section, resummed 2-loop perturbation theory indicates that the electroweak phase transition can indeed be strong enough for baryogenesis. - To check the reliability of perturbation theory, a 3d effective theory was constructed in refs. . The Lagrangian in Eq. (4) is supplemented by the right-handed stop field $`U`$ and the SU(3) gauge fields, with the field strength tensor $`G_{ij}`$:

$$\delta _{\text{3d}}=\frac{1}{2}\mathrm{Tr}G_{ij}^2+(D_i^sU)^{}D_i^sU+m_{U3}^2U^{}U+\lambda _{U3}(U^{}U)^2+\gamma _3\varphi ^{}\varphi U^{}U.$$
(8)
Simulations with this action have been carried out in ref. , and with a more complicated version thereof, involving two Higgs doublets, in ref. .

The results are illustrated in Fig. 2 for a particular choice of parameters. We find the encouraging outcome that, in fact, for a light right-handed stop the 2-loop estimates are reliable and even somewhat conservative. This is in strong contrast to the case of the SM at realistic Higgs masses.

Next, we must ask more rigorously whether this parameter region is indeed in agreement with all experimental data. The most important constraint comes from the lower bound on the Higgs mass.  There is a parameter in the MSSM, $`m_A`$, which determines whether the Higgs sector resembles that in the Standard Model ($`m_A>120`$ GeV) or not ($`m_A<120`$ GeV). In the latter case, the experimental lower bound is relaxed,  and since the transition needs not always get significantly weaker (see Fig. 3, and ref.  for more details), this case is acceptable. The other possibility is having a large $`m_A`$, but then the left-handed stop should be quite heavy, $`m_{\stackrel{~}{t}_L}>2`$ TeV (see Fig. 3), in order to increase the Higgs mass.

## 8 CP violation in the background field configuration

So far, we have discussed the strength of the phase transition. Let us finally very briefly touch another issue relevant for baryogenesis, CP violation.

Even though the strength of the phase transition is determined by a single “dynamical” Higgs doublet, in the MSSM there are of course in principle two of them, of which the one discussed above, $`\varphi `$ in Eqs. (4), (8), was a particular linear combination (see, e.g., ref. ). However, the direction orthogonal to $`\varphi `$ could play a role for CP violation. Indeed, when there are two Higgses, $`H_1,\stackrel{~}{H}_2`$, then there is for instance a new low-dimensional gauge invariant CP-odd observable $`ImH_1^{}\stackrel{~}{H}_2`$. It could obtain a non-vanishing value, in principle even a large one, as is the case in spontaneous CP violation.  If so, then at the semi-classical level, fermion mass matrices may obtain non-trivial space-dependent phases, so that particles and anti-particles have different dispersion relations, which could be of relevance to baryon asymmetry generation.

However, in practice such interesting effects are quite small. The reason is that all the effects related to the orthogonal direction are suppressed at least by $`g^2T/[\pi (m_A^2+0.5T^2)^{\frac{1}{2}}]`$. Here $`m_A`$, experimentally $`>m_Z`$ measures the deviation from the single Higgs doublet dynamics at zero temperature, and the temperature correction is seen to increase the suppression further on.

This situation is illustrated in Fig. 4, which shows the profile of a stationary phase interface (for perturbative analyses of such objects, see ref.  and references therein). We observe that in field space, the transition takes place in one particular direction (see also ), and CP violation is small, as if the system effectively only had one dynamical CP even Higgs doublet.

## 9 Conclusions

We have reviewed the methods that can be employed in determining the thermodynamical properties of the phase transition related to electroweak gauge symmetry restoration, as well as the numerical results that are already available for the Standard Model and for its Minimal Supersymmetric extension. The general statement to be remembered, is that there is no phase transition at all in electroweak theories with realistic Higgs masses, unless there is also another light scalar degree of freedom available, which plays an essential role in the phase transition dynamics.

As a concrete phenomenological example, we have recalled that it is still possible to find parts of the MSSM parameter space, not excluded by experiment, where the electroweak phase transition is strongly enough of the first order to allow for baryon asymmetry generation. However, these parts are quite restricted: not only does one need one “light” stop, $`m_{\stackrel{~}{t}_R}<m_{\text{top}}`$, in order to get a strong transition, but also either another stop which is much heavier, $`m_{\stackrel{~}{t}_L}>10m_{\text{top}}`$, or a Higgs sector which has at least two independent light particles, $`m_A<120`$ GeV, in order not to violate experimental constraints.

## Acknowledgments

I thank D. Bödeker and K. Rummukainen for useful comments. This work was partly supported by the TMR network Finite Temperature Phase Transitions in Particle Physics, EU contract no. FMRX-CT97-0122, and by the RTN network Supersymmetry and the Early Universe, EU contract no. HPRN-CT-2000-00152.


text:
# Diabolical Points in Molecular Magnets with a Four-Fold Easy Axis

## I INTRODUCTION

The magnetic molecules Mn12 \[short for Mn12-acetate, or \[Mn12(CH3COO)16(H2O)4O12\] $``$2CH3COOH$``$4H2O\] and Fe8 \[short for Fe8O2(OH)12(tacn)6\]8+\] are among a few dozen that are currently being studied as extreme cases of superparamagnets . Both Mn12 and Fe8 have spin 10, and both of them display hysteresis at the molecular level , as do some of the others. In addition, however, Fe8 shows an effect that has not yet been seen in any of the other molecules: oscillation of the Landau-Zener-Stückelberg transition rate between low lying levels as a function of the applied static magnetic field . This oscillation is due to an oscillatory quenching of the underlying tunneling matrix element connecting the levels in question, and is an unambiguous signature of quantum tunneling. The occurrence and observability of tunneling in a spin of such a large magnitude is of much interest in itself, and the tunneling frequency, of the order of 100 Hz, is perhaps the lowest ever inferred or meaasured in any physical system. Further, the quenching effect can be regarded as due to the interference of different Feynman tunneling paths for the spin , and this is how it was first discovered . While massive particle tunneling in two or more spatial dimensions can also show such interference , the effect arises more directly in the spin problem since the kinetic term in the action has the mathematical structure of a Berry phase. This adds to the interest in the problem. Reciprocally, the experimental observations have motivated more careful investigations of spin-coherent-state path integrals which are more delicate than their massive particle counterparts .

The above intereference effect has also been sought in Mn12, but has not yet been seen. While the spin Hamiltonian for Fe8 has biaxial symmetry, Mn12 is tetragonal, hence the systematics of the effect are different, and it is interesting to calculate them. More specifically, with an external magnetic field $`𝐇`$, Mn12 is described by an anisotropy Hamiltonian

$$=AS_z^2BS_z^4+C(S_+^4+S_{}^4)g\mu _b𝐒𝐇,$$
(1)
where $`S=10`$, $`g2`$, and $`ABC>0`$. (The experimental values for $`A`$, $`B`$, and $`C`$ are 0.556, $`1.1\times 10^3`$, and $`3\times 10^5`$ K, respectively, so that $`\lambda _1=B/A=1.98\times 10^3`$, $`\lambda _2=C/A=5.4\times 10^5`$, and $`H_c=A/g\mu _B=0.414`$ T.) The easy axis is now $`z`$ with four fold symmetry, the hard axes are $`\pm x`$ and $`\pm y`$, and the medium axes are the lines $`y=\pm x`$ in the $`xy`$-plane. For values of $`|𝐇|`$ that are not too large, the spin coherent state  expectation value $`\widehat{𝐧}||\widehat{𝐧}`$, which may be regarded as a classical energy $`E(\theta ,\varphi )`$, may have two or more local minima, between which the spin can then tunnel. The issue is to calculate the tunnel splitting $`\mathrm{\Delta }`$, and especially the field values where this splitting is quenched, i.e., $`\mathrm{\Delta }=0`$. Since there may be other molecules with this symmetry, it is desirable to do this analytically, for all values of (or at least a wide range of) $`\lambda _1`$ and $`\lambda _2`$.

For the discussion that follows, it is useful to review some basic facts pertaining to degeneracy in quantum mechanics in the absence of symmetry. A point of degeneracy, or equivalently, a point where the splitting between two levels vanishes, is diabolical  in the terminology of Berry and Wilkinson , or a conical intersection in older terminology . As a rule, eigenvalues of a finite Hamiltonian are all simple, and for a general Hamiltonian, represented by a complex Hermitean matrix, we must be able to adjust at least three parameters in order to produce a degeneracy. A simple approximate argument is as follows . Let two states $`m`$ and $`m^{}`$ be approximately degenerate, and let the secular matrix between them be written as

$$\left(\begin{array}{cc}E_m& V_{mm^{}}\\ V_{m^{}m}& E_m^{}\end{array}\right),$$
(2)
with $`V_{m^{}m}=V_{mm^{}}^{}`$. The states will be truly degenerate only if the following two conditions are met:

$`B_{mm^{}}`$ $``$ $`E_mE_m^{}=0,`$ (3)
$`V_{mm^{}}`$ $`=`$ $`0.`$ (4)
It is convenient to refer to these as the no-bias and the no-mixing conditions, respectively . Since $`V_{mm^{}}`$ is in general complex, we have three real conditions, requiring three or more variable parameters for their satisfaction.

Precisely three parameters are available to us in Eq. (1) in the three components of $`𝐇`$. If the Hamiltonian matrix is real symmetric, the number of adjustable parameters required is lowered to two. In the present problem, this situation is realized when $`H_y=0`$, and so, as in the Fe8 case, we expect to find the degeneracies in the $`H_x`$-$`H_z`$ plane. Unless explicitly stated, we will henceforth take $`H_y=0`$, so that $`V_{mm^{}}=V_{m^{}m}`$.

Ignoring an additive constant, the energy surface in the vicinity of the degeneracy is given by

$$E=\pm (B_{mm^{}}^2+V_{mm^{}}^2)^{1/2}$$
(5)
which has the form of a double cone or a diabolo in $`H_x`$-$`H_z`$ space, which explains the term diabolical for these points.

In a previous paper , we have studied the Mn12 problem for $`𝐇\widehat{𝐱}`$. In this case, $`E(\theta ,\varphi )`$ is reflection symmetric in the equatorial plane, and for small enough $`H`$, has two degenerate minima, one in each hemisphere with $`S_z>0`$ and $`S_z<0`$. The problem is analogous to the tunneling of a massive particle in a symmetric double well. The approach used is a discrete phase integral (DPI) or Wentzel-Kramers-Brillouin (WKB) method, and it provides a good quantitative approximation for the tunnel splitting, $`\mathrm{\Delta }`$, as well as the diabolical fields where the splitting vanishes.

In this paper, we extend our studies to the case where the field has a nonzero $`z`$ component. The problem is now like a massive particle in an asymmetric double well. Our initial intent was to approach this case also using the DPI method, as it has been appplied successfully to the Fe8 problem. Indeed, the calculation reduces to no more than the evaluation of a handful of action integrals, and based on prior experience with Fe8 and other tunneling Hamiltonians, we are confident that it will be quantitatively accurate for Mn12 too. Unlike the Fe8 case, however, where the integrals can be found analytically, this turns out to be not so for Mn12, and so, even though the calculation is not as atomistic as a brute force diagonalization of the $`21\times 21`$ Hamiltonian matrix, it still does not yield final answers in an analytic form. To our pleasant surprise, however, we have found that a perturbative approach in the spirit of Ref.  not only yields quite good approximations for the locations of the diabolical points, it also reveals a pattern in the points which would be hard to discern from a numerical diagonalization by itself. Most importantly, it provides a scheme for indexing these points. It is the purpose of this paper to report this work.

The analysis is contained in Sec. II. We first develop the perturbation theory with $`B=0`$ (Sec. II A). This yields the diabolical points as the roots of polynomials in $`H_x`$ and $`C`$. For a given spin $`S`$, we have $`2S`$ polynomials. One of the unexpected bonuses is that a polynomial which applies to a given value of $`S`$ also applies to any other value of $`S`$. We find all these polynomials for $`S10`$. In subsection II B, we incorporate the effects of the $`B`$ term approximately, and compare our analytic results with those from explicit numerical diagonalization of the Hamiltonian with the parameters for Mn12 (see Table III). Our results are accurate to about 10% for Mn12, and we believe that they will also be useful for other systems with four-fold anisotropy. This is especially true for the low lying energy levels. In fact, for the higher pairs of levels, the diabolical points can be significantly moved or even eliminated altogether by still higher ansitropy terms in the Hamiltonian. In subsection II C we discuss some qualitative aspects of the degeneracies on the $`H_z`$ axis, and show that some of them behave as the merger of two or three diabolical points. A short summary (Sec. III) concludes the paper.

## II Perturbative Calculation

Let us first consider the problem when $`H_x=0`$, i.e., $`𝐇\widehat{𝐳}`$. The $`S_\pm ^4`$ terms couple states with $`m`$ differing by 4, and thus divide the Hamiltonian into four disjoint subspaces. Levels belonging to different subspaces can cross at fields which can be approximately obtained by neglecting the $`S_\pm ^4`$ terms. The crossing conditions are $`E_m=E_m^{}`$, where $`E_m=Am^2Bm^4g\mu _BH_zm`$. Improved formulas can be obtained by finding corrections to $`E_m`$ perturbatively in $`C`$. These intersections are easy to understand because a symmetry of the Hamiltonian (invariance under a rotation by 90 about $`\widehat{𝐳}`$) is easily recognized.

If $`H_x`$ and $`H_z`$ are both non-zero, there is no obvious symmetry. If $`H_x`$ and $`C`$ are both small, however, we may continue to label the states by the $`m`$ quantum numbers, and calculate the energies perturbatively in these two parameters. In terms of the secular matrix (2), the energies $`E_m`$, $`E_m^{}`$, and the bias $`B_{mm^{}}`$ are determined by the terms in $`S_z`$ in Eq. (1), and the mixing $`V_{mm^{}}`$ by the terms involving $`C`$ and $`H_x`$. The energies are trivial to find, so the problem is to find $`V_{mm^{}}`$.

It is convenient to divide the Hamiltonian by $`A`$ and to work with scaled quantities

$$\lambda _1=B/A,\lambda _2=C/A,h_{x,z}=H_{x,z}/SH_c,$$
(6)
with

$$H_cA/g\mu _B.$$
(7)
### A Simplified model: $`B=0`$.

To keep the problem tractable, let us set $`B=0`$ at this stage. Then, to zeroth order in both $`C`$ and $`H_x`$, $`E_m=Am^2g\mu _BH_zm`$. Hence, levels $`m`$ and $`m^{}`$ are degenerate when

$$h_z=\frac{m+m^{}}{S}.$$
(8)
It remains to find the off-diagonal element $`V_{mm^{}}`$. As we shall see, the choice $`B=0`$ simplifies the calculation greatly, for $`E_m`$ is then quadratic in $`m`$, and energy level differences are linear in $`m`$, and given by a fixed set of numbers whenever Eq. (8) holds.

To illustrate the calculation of $`V_{mm^{}}`$, we consider the case $`S=5`$. Suppose $`h_z=1/5`$, so that $`m=5`$ and $`m^{}=4`$ are degenerate (see Fig. 1). We will find $`V_{4,5}`$ to leading non-zero order in $`h_x`$ and $`C`$ as a double series in these variables. It is clear that a transition from $`m`$ to $`m^{}`$ can be made in three ways:

(1) act with $`h_xS_+`$ in 9th order.

(2) act with $`h_xS_+`$ in 5th order and $`CS_+^4`$in 1st order.

(3) act with $`h_xS_+`$ in 1st order and $`CS_+^4`$ in 2nd order.
We denote the corresponding contributions to $`V`$ by $`V^{(1)}`$, $`V^{(2)}`$, and $`V^{(3)}`$. Each of these involves a product of matrix elements and a product of energy denominators. For $`V^{(1)}`$, the former is

$$\frac{1}{2^9}5,4S_+^95,5(h_xS)^9=709\sqrt{10}(h_xS)^9W(h_xS)^9.$$
(9)
In fact the number $`W`$ will be common (up to some power of 2) to all three $`V^{(i)}`$. The energy denominators can be read off Fig. 1, and for $`V^{(1)}`$, these are

$$(1)^8(8\times 14\times 18\times 20)^2(1)^8K^2.$$
(10)
The factor $`(1)^8`$ appears here because all intermediate states are higher in energy than $`E_5`$ and $`E_4`$, and we have introduced the number $`K`$ because many of the energy denominator products for $`V^{(2)}`$ and $`V^{(3)}`$ contain the same factors. Putting together Eqs. (9) and (10),

$$V^{(1)}=\frac{W}{K^2}(h_xS)^9.$$
(11)
For $`V^{(2)}`$, the transition can occur in several ways, six to be precise, corresponding to where the $`S_+^4`$ term acts. Two of the ways to make the transition are $`5101234`$, and $`5401234`$. The product of matrix elements in each case is

$$\left(\frac{h_xS}{2}\right)^5\lambda _2<4S_+^95>=2^4W\lambda _2(h_xS)^5.$$
(12)
The energy denominators, however, depend on the transition path, and are listed in Table I. Adding together all the contributions, we get

$`V^{(2)}`$ $`=`$ $`2^4W{\displaystyle \frac{2}{K}}\left({\displaystyle \frac{1}{8}}+{\displaystyle \frac{1}{20}}+{\displaystyle \frac{20}{112}}\right)(h_xS)^4\lambda _2`$ (13)
$`=`$ $`{\displaystyle \frac{2^4\times 99}{140}}{\displaystyle \frac{W}{K}}(h_xS)^5\lambda _2.`$ (14)
Lastly, for $`V^{(3)}`$, there are three transition paths: (i) $`5134`$, (ii) $`5104`$, (iii) $`5404`$. The transition element product for all three is

$$2^8W(h_xS)\lambda _2^2,$$
(15)
and the energy denominator product is $`20\times 8`$, $`20\times 20`$, and $`8\times 20`$, respectively. Thus,

$$V^{(3)}=\frac{2^8\times 3}{200}W(h_xS)\lambda _2^2.$$
(16)
Adding together $`V^{(1)}`$, $`V^{(2)}`$, and $`V^{(3)}`$, we obtain the net $`V_{4,5}`$ (restoring the level quantum numbers). The condition for a diabolic point is that this quantity should vanish. In addition to $`h_x=0`$, this happens when

$$\xi ^2\frac{99K}{140}\xi +\frac{3K^2}{200}=0,$$
(17)
where

$$\xi =\frac{1}{\lambda _2}\left(\frac{h_xS}{2}\right)^4.$$
(18)
Solving Eq. (17), we obtain

$$h_x=\frac{2}{5}\left\{\frac{99}{280}\pm \left[\left(\frac{99}{280}\right)^2\frac{3}{200}\right]^{\frac{1}{2}}\right\}^{\frac{1}{4}}(K\lambda _2)^{\frac{1}{4}}.$$
(19)
With the scaled value $`\lambda _2=2.16\times 10^4`$ for $`S=5`$, this yields $`h_x=0.2643`$, and 0.6252. Direct numerical diagonalization yields $`h_x=0.2669`$, and 0.638.

Readers will undoubtedly have noted that apart from an overall factor of $`h_x`$ to some power, our perturbation method yields the off-diagonal element as a homogeneous polynomial in $`h_x^4`$ and $`\lambda _2`$. It is not difficult to see that this will be generally true, and also not difficult to justify. Let us first take the point that we have only included transition paths that go through the higher energy levels. Consider, e.g. the path for $`V_{4,5}`$ in the above calculation that goes from $`5`$ to $`+1`$ via six successive $`h_xS_+`$ terms, then to $`+5`$ via a $`CS_+^4`$ term, and finally to $`+4`$ via an $`h_xS_{}`$ term. This term is of order $`h_x^7\lambda _2`$, and should be regarded as a higher order correction to $`V^{(2)}`$. Secondly, it is positive and of the same sign as $`V^{(2)}`$, because it involves six negative and one positive energy denominator. This feature is also generally valid, and is important in light of the next point, which is that the sign of the terms in the polynomial for the off-diagonal element alternate when organized as a series in $`\lambda _2`$. Thus, for $`V_{4,5}`$, $`V^{(1)}`$ is negative, $`V^{(2)}`$ is positive, and $`V^{(3)}`$ is negative. This is a consequence of the fact that replacing four $`h_xS_+`$ terms by a single $`CS_+^4`$ term $`(a)`$ leaves the sign of the matrix element product unchanged, but $`(b)`$ replaces four negative energy denominators by a single negative one.

Let us call the polynomial that remains after we have cancelled off as many overall factors of $`h_xS`$ from $`V_{mm^{}}`$ as possible the underlying polynomial. It is clear that this polynomial is of degree

$$n_{mm^{}}=\left[\left|\frac{mm^{}}{4}\right|\right],$$
(20)
in $`h_x^4`$, where $`[x]`$ denotes the integer part of $`x`$, i.e., the largest integer less than or equal to $`x`$. The alternation of signs of successive powers of $`h_x^4`$ is a necessary (but not sufficient) condition for all $`n_{mm^{}}`$ roots to be positive . This means that, not including the points on the $`h_x`$ or $`h_z`$ axes, it is possible for states labelled by $`m`$ and $`m^{}`$ to intersect in a diabolical point up to $`n_{mm^{}}`$ times in the first quadrant of the $`h_x`$-$`h_z`$ plane. This appears to us to be a topological property of the Mn12 Hamiltonian, that is not altered by presence of higher order terms, as long as the symmetry is not changed. Of course, the number of diabolical points may be fewer, but we do not believe it can be greater, because if $`h_x`$ is sufficiently large, the term $`H_xS_x`$ dominates the energy in the equatorial plane, and the possibility of interfering trajectories is lost. We do not have a proof of these statements, which must be regarded as conjectures, but the similarity to Fe8, and all the empirical evidence we have gathered suggests that they are indeed true.

For $`S=5`$, we have found all the diabolical points using this perturbation approach, and also numerically. In all cases, the perturbative answers are nearly exact. The results are shown in Fig. 2.

At this point we wish to note a remarkable feature of this approximation, which may have been noticed by some readers. This is that the diabolical values of $`h_x`$ depend on $`m`$ and $`m^{}`$ only through the combination $`\mathrm{\Delta }m=m^{}m`$. This means that in Fig. 2, the theoretical points corresponding to the same value of $`\mathrm{\Delta }m=m^{}m`$ are vertically aligned. See, for example, the points corresponding to $`m=4`$, $`m^{}=4`$, and $`m=5`$, $`m^{}=3`$, or to $`m=5`$, $`m^{}=2`$, and $`m=4`$,$`m^{}=3`$. The reason is not hard to find. When we set $`B=0`$, the energy $`E_m`$ is a quadratic function of $`m`$, and when levels $`m`$ and $`m^{}>m`$ are degenerate,

$$E_{m+k}E_m=k(\mathrm{\Delta }mk).$$
(21)
Thus, the entire pattern of energy levels above the levels $`m`$ and $`m^{}`$ depends only on $`\mathrm{\Delta }m`$ (see Fig. 1), and since only these levels enter into our perturbation theory, the energy denominators are identical. The matrix elements are of course different, but since our transition paths involve no closed loops, they amount to a net factor of $`m^{}|S_+^{\mathrm{\Delta }m}|m`$ in each term, which therefore drops out of the underlying polynomial. In short, this polynomial depends only on $`\mathrm{\Delta }m`$.

Furthermore, we see from Eq. (21) that value of $`S`$ also does not enter in the energy denominators. This means that the polynomials found for $`S=5`$ are applicable to the $`S=10`$ problem for transitions with $`\mathrm{\Delta }m10`$, and makes it worthwhile to find all the remaining polynomials for $`S=10`$. The task is easily automated on a computer since the problem is essentially to enumerate the transition paths for a given order in $`h_x`$ and $`\lambda _2`$, and to add up the reciprocals of the corresponding energy denominators. It can be further simplified by noting that all denominators appear in the term proportional to $`h_x^{\mathrm{\Delta }m}`$, so that relative to this term, the coefficient for any other transition path appears as an energy numerator, consisting of all the missing denominators. The polynomials are given in Table II , along with the roots for $`h_x`$ for general $`S`$ and $`\lambda _2`$, as well as for the values applicable to Mn12.

One last general point worth noting is that for a diabolical point labelled by the pair $`(m,m^{})`$ with $`m^{}>m`$, the energy levels which are degenerate are numbers $`2S+(mm^{})+1`$ and $`2S+(mm^{})+2`$, where the ground state is given the number 1.

### B Inclusion of $`S_z^4`$ term

When we try and compare the results of the previous subsection with those obtained numerically from Eq. (1) with $`B0`$, we find that the systematics of the diabolical points are fully captured in that the analytic results provide a complete indexing scheme, but the quantitative disagreement with the Mn12 parameters is as bad as $`30\%`$ in some cases. We therefore seek some way to incorporate the $`B0`$ effects.

It is easy to let $`B0`$ in the no-bias condition. Equation (8) is modified to

$$h_z(m,m^{})=\frac{1}{S}(m+m^{})\left[1+\lambda _1(m^2+m^2)\right].$$
(22)
The no-mixing condition is clearly harder to evaluate. A simple minded approach is to shift the energy levels so as to retain the same relative spacings as when $`B=0`$, but allow the overall range to be modified.

With this in mind, let us first consider the energies when $`\lambda _1=0`$. When levels $`m`$ and $`m^{}`$ are degenerate, the level at the top of the barrier is given by the quantum number $`k=(m+m^{})/2`$, whenever this is an integer, or by the nearest integers if it is a half-integer. In the spirit of our approximations, keeping track of this distinction would be an overrefinement, so we will use the formula $`(m+m^{})/2`$ in both cases. The energy range is thus given by

$$\mathrm{\Delta }E^{(0)}=E_k^{(0)}E_m^{(0)}=\frac{1}{4}(mm^{})^2,$$
(23)
where the (0) superscript indicates that $`B=0`$. With $`B0`$, we get

$`\mathrm{\Delta }E^{(1)}`$ $`=`$ $`\mathrm{\Delta }E^{(0)}\lambda _1\left[\left({\displaystyle \frac{m+m^{}}{2}}\right)^4m^4\right]+\lambda _1(m+m^{})(m^2+m^2)\left[{\displaystyle \frac{m+m^{}}{2}}m\right]`$ (24)
$`=`$ $`\gamma _{mm^{}}\mathrm{\Delta }E^{(0)},`$ (25)
where

$$\gamma _{mm^{}}=1+\frac{\lambda _1}{4}(7m^2+10mm^{}+7m^2).$$
(26)
If we assume that the entire spectrum gets modified from its quadratic form by a uniform stretching factor $`\gamma _{mm^{}}`$, then the only change in our perturbation theory is that all energy denominators get multiplied by this factor. In the underlying polynomial, $`h_x`$ and $`\lambda _2`$ get replaced by $`h_x/\gamma _{mm^{}}`$ and $`\lambda _2/\gamma _{mm^{}}`$, and hence the no-mixing condition becomes

$$h_x(m,m^{})=\gamma _{mm^{}}^{3/4}r_\alpha (\mathrm{\Delta }m),\alpha =1,2,\mathrm{},n_{mm^{}},$$
(27)
where $`\{r_\alpha \}`$ are the original $`h_x`$ values obtained from the roots of the underlying polynomial $`P_{\mathrm{\Delta }m}`$.

The formulas (22) and (27) are compared with exact numerical results in Table III. The errors are now typically about 10%, and can be of either sign.

It is useful to briefly discuss our numerical procedure. For points lying on the $`H_z`$ or $`H_x`$ axis, the splitting is a function of one variable, and its zeros can be found by simple scanning. For the off-axis zeros, this is harder, and we resort to the Herzberg and Longuet-Higgins sign change theorem , which applied to the present problem states that, upon adiabatic traversal of a closed contour in the $`H_x`$-$`H_z`$ plane enclosing a single point of degeneracy of two states, the wavefunction of either of these two states returns to itself except for a change in sign. Conversely, there is no change in sign if the contour does not enclose a degeneracy. Hence, to find a diabolical point, we first find a sign-reversing rectangular contour by hit and trial. By bisecting this rectangle in the $`x`$ and $`z`$ directions alternately, and using the sign-change test at each bisection, we can corral the degeneracy ever more tightly, to the degree desired. We have found this procedure to be generally superior to a direct minimization of the energy difference for the reason that the diabolo in the vicinity of a degeneracy is highly asymmetrical in the $`x`$ and $`z`$ directions. Consider for example, the 4th and 5th energy levels from the bottom when $`h_z0.12`$$`0.13`$, corresponding approximately to the $`m`$ quantum numbers $`9`$ and $`8`$. Since these states are separated by a high barrier, the mixing element between them is best understood as arising from tunneling, and will therefore be proportional to an exponentially small Gamow factor $`e^\mathrm{\Gamma }`$, where $`\mathrm{\Gamma }`$ is the appropriate tunneling action. Thus the energy surface consists of a deep and narrow valley running nearly parallel to the $`h_x`$ axis, with a valley floor that goes to zero linearly at occasional points, and may rise and fall in between. Because of this shape, and because the surface is not analytic in the vicinity of the points being sought, standard methods for finding the minima of a function are often not well suited.

The above argument also enables to understand an observation made by Berry and Wilkinson  and Berry and Mondragon  in the study of two very different model problems, namely, that the energy cone near a diabolical point has very high eccentricity in terms of the physically natural parameters describing the system. In other words, the cross section of the energy surface is a very long and narrow ellipse in one direction. Let us see how this happens in the present problem. To save writing let us write just $`x`$ and $`z`$ for the deviations of $`h_x`$ and $`h_z`$ from a diabolical point at $`(h_{x0},h_{z0})`$. In the vicinity of this point, we have

$`B_{mm^{}}`$ $``$ $`az,`$ (28)
$`V_{mm^{}}`$ $``$ $`be^\mathrm{\Gamma }x,`$ (29)
where $`a`$ and $`b`$ are constants of order unity. Thus the energy surface is

$$E(a^2z^2+b^2e^{2\mathrm{\Gamma }}x^2)^{1/2}.$$
(30)
The cross section is an ellipse with major axis parallel to $`x`$ and eccentricity $`e^\mathrm{\Gamma }1`$. This scenario is expected to be quite general. The no-bias condition defines a line in parameter space. The gradient of the bias normal to this line is generally expected to be of order unity in the natural physical variables. The mixing element also varies on the same order unity scale in the parameter space, but because it arises from tunneling, its absolute scale is very small. The result is an energy surface of high eccentricity of the type just described.

### C Merged diabolical points

Our discussion of diabolical points needs some elaboration for certain degeneracies lying on the $`h_z`$ axis. For the pair $`(m,m^{})=(10,9)`$, e.g., (more generally any pair with $`\mathrm{\Delta }m=4n+3`$), we peeled off a factor of $`h_x^3`$ from the mixing element. Writing $`x=h_x`$, and $`z=h_z`$ as in Eq. (30), and $`z_0=h_{z0}`$ for the point of degeneracy, the bias and mixing are given by

$`B_{mm^{}}`$ $``$ $`zz_0,`$ (31)
$`V_{mm^{}}`$ $``$ $`x^3+O(x^7),`$ (32)
ignoring multiplicative constants. Correspondingly, the energy surface is $`[(zz_0)^2+x^6]^{1/2}`$, whose cross section is no longer an ellipse. There is also no reason for the simple sign-change result to hold a priori.

These conclusions are based on perturbation theory, however. More generally, we can only argue on grounds of symmetry that, under $`xx`$, $`B_{mm^{}}B_{mm^{}}`$, and $`V_{mm^{}}\pm V_{mm^{}}`$. Instead of Eqs. (31) and (32), we should therefore expect the general expansion to take the form

$`B_{mm^{}}`$ $``$ $`zz_0+ax^2+O((zz_0)^2,x^2(zz_0)^2,x^4),`$ (33)
$`V_{mm^{}}`$ $``$ $`x^3bx+cx(zz_0)+O(x^3(zz_0),x^5),`$ (34)
where $`a`$, $`b`$, and $`c`$ are constants, all of which we expect to be very small on account of the quantitative accuracy of the perturbative approach. Ignoring the higher order terms, the bias vanishes on the parabola $`z=z_0ax^2`$. On this parabola, the mixing is given by

$$(1ac)x^3bx,$$
(35)
which vanishes at $`x=0`$, and $`x=\pm [b/(1ac)]^{1/2}\pm x_1`$, assuming that $`b/(1ac)>0`$. Thus instead of a single degeneracy at $`(0,z_0)`$, we have three closely spaced degeneracies at

$$(0,z_0),(\pm x_1,z_0ax_1^2).$$
(36)
The energy surface in the immediate vicinity of each one of these points is now diabolical in the ordinary sense. A small circuit of each of these points separately will therefore lead to a sign reversal, as will a larger circuit enclosing all three of them. If we ignore the splitting, we may regard the original degeneracy on the $`h_z`$ axis as a triply merged diabolical point.

It is useful to think of the coefficients $`a`$, $`b`$, and $`c`$ as depending on parameters in the Hamiltonian other than the components of $`𝐇`$, e.g., $`\lambda _1`$ and $`\lambda _2`$. It may be that as these parameters are varied, the quantity $`b/(1ac)`$ becomes negative, so that the roots $`\pm x_1`$ cease to be real. We can think of the two off-axis diabolical points as having annihilated each other, leaving behind only one true diabolical point on the axis. Unless one is very close to this point, however, the energy surface may still resemble that of a triply merged point.

For the parameters appropriate to Mn12, we find that the points are located at

$$(h_x,h_z)=(0.0,0.135836224),(\pm 0.01855,0.135832551).$$
(37)
These numbers are obtained by using the same sign-reversal theorem as previously described. Because the energy difference depends so sensitively on $`h_z`$, however, we have confirmed them in another way. For any given value of $`h_x`$, we first find the minimum of the relevant energy gap $`\mathrm{\Delta }`$ with respect to $`h_z`$. In essence, we find the value of the gap at the bottom of the parabolic trench where the bias vanishes. A plot of this gap versus $`h_x`$ should be given by the absolute value of the expression (35). As shown in Fig. 3, this is indeed so, and the split off point is again found to be at $`h_x=0.01855`$.

In exactly the same way, we may also consider doubly merged points, corresponding to the tunneling of states with $`\mathrm{\Delta }m=4n+2`$, e.g., $`(m,m^{})=(10,8)`$, when $`h_x0`$. We can continue to expand the bias as in Eq. (33), but the leading term in $`V_{mm^{}}`$ is now proportional to $`x^2`$, so that instead of Eq. (34), we have

$$V_{mm^{}}x^2+O((zz_0),x^2(zz_0),x^4).$$
(38)
It is then obvious that the double zero of $`V_{mm^{}}`$ at $`x=0`$ can not be split. This conclusion can also be reached in another way. Symmetry would require that, if they split, the points be located at $`(\pm x_0,z_0)`$, with $`x_00`$. The selection rule argument given at the very beginning of Sec. I shows, however, that this is impossible, as there must be a crossing of any two levels $`m`$ and $`m^{}`$ with $`\mathrm{\Delta }m4n`$ as $`H_z`$ is varied with $`H_x=0`$. In Fig. 3, we also show the $`(m,m^{})=(10,8)`$ gap at the bottom of the no-bias trench. It is apparent that now the diabolical points remain unsplit at $`h_x=0`$, and the curve is extremely well fit by a parabola, as required by Eq. (38).

Finally, the points on the $`h_z`$ axis, corresponding to $`\mathrm{\Delta }m=4n+1`$ are singly diabolical to begin with, so the issue of splitting does not arise.

## III Summary

We have studied the diabolical points of a spin Hamiltonian describing molecules such as Mn12 that have an easy axis of four-fold symmetry, for arbitrarily directed magnetic fields. A perturbation theory in the parameters $`H_x`$ and $`C`$ is found to give a very good qualitative and even quantitative understanding. Our central results are the formulas (22) and (27), which along with the results in Table II give the full set of diabolical points for any molecule with $`S10`$.

###### Acknowledgements.

We are indebted to Wolfgang Wernsdorfer for very useful correspondence. AG’s research is supported by the NSF via grant number DMR-9616749.