text:
# Quantum Tunneling of Magnetization in a Large Molecular Nanomagnet – Approaching the Mesoscale

## Abstract

A Mn30 molecular cluster is established to be the largest single-molecule magnet (SMM) discovered to date. Magnetization versus field measurements show coercive fields of about 0.5 T at low temperatures. Magnetization decay experiments reveal an Arrhenius behavior and temperature-independent relaxation below 0.2 K diagnostic of quantum tunneling of magnetization through the anisotropy barrier. The quantum hole digging method is used to establish resonant quantum tunneling. These results demonstrate that large molecular nanomagnets, having a volume of 15 nm3, with dimensions approaching the mesoscale can still exhibit the quantum physics of the microscale.

The study of the interface between classical and quantum physics has always been a fascinating area, but its importance has nevertheless grown dramatically with the current explosive thrusts in nanoscience. Taking devices to the limit of miniaturization (the mesoscale and beyond) where quantum effects become important makes it essential to understand the interplay between the classical properties of the macroscale and the quantum properties of the microscale. This is particularly true in nanomagnetism, where many potential applications require monodisperse, magnetic nanoparticles. One source of such species are single-molecule magnets (SMMs) Christou00 ; Sessoli93b ; Sessoli93 ; Aubin96 ; Boskovic02 , individual molecules that function as single-domain magnetic particles. Below their blocking temperature, they exhibit magnetization hysteresis, the classical macroscale property of a magnet, as well as quantum tunneling of magnetization (QTM) Friedman96 ; Thomas96 ; Sangregorio97 ; Aubin98 ; Aubin98b ; Boskovic01  and quantum phase interference WW\_Science99 ; Garg93 , the properties of a microscale entity. QTM is advantageous for some potential applications of SMMs, e.g. in providing the quantum superposition of states for quantum computing Leuenberger01 , but it is a disadvantage in others such as information storage where it would lead to loss of preferential spin alignment. Large SMMs approaching the mesoscale have long been the target of synthesis, and an important question for these (and other mesoscale magnetic particles) is whether they might still unequivocally demonstrate quantum properties, as theoretically predicted Chudnovsky88 ; Barbara90 .

In this letter, we show that a Mn30 molecular cluster is by far the largest SMM to date, and show that it unambiguously still exhibits QTM. This establishes that the quantum physics of the microscale can also be observed in large molecular nanomagnets.

The Mn30 molecular cluster has the formula \[Mn30O24(OH)8(O2CCH2CMe3)32(H2O)2(MeNO2)4\]. xCH2Cl2.yH2O and crystallizes in the monoclinic space group C2/c with each Mn30 on a crystallographic C2 rotation axis Soler01 . The molecules are thus all parallel with their molecular C2 axes ($`z`$ axes) aligned along the crystal $`b`$ axis. Each Mn30 comprises a Mn/O magnetic core (Fig. 1) of volume $``$4 nm3 enveloped in a hydrocarbon shell $``$0.5 nm thick composed of carboxylate CH2CMe3 groups. This essentially magnetically isolates neighboring Mn/O cores. However, the large size of Mn30 (volume $``$15 nm3, molar mass $``$6335 Daltons), akin to a small protein, causes inefficient packing in the crystal, with voids between molecules that contain disordered solvent molecules, primarily CH2Cl2. In addition, some CMe3 groups are statically disordered about two sites.

Each Mn30 contains one Mn4+, twenty-six Mn3+ and three Mn2+ ions, and exchange coupling between the many spins of the Mn/O core leads to Mn30 having a small net uncompensated spin in the ground state of S = 5. The latter was determined from temperature and field magnetization measurements on powdered, microcrystalline samples in the 1.8-4.0 K and 0.1-0.4 T ranges. The data were fit by diagonalization of the spin Hamiltonian matrix assuming only the ground state is populated and incorporating axial anisotropy $`(DS_z^2)`$ and Zeeman terms, and employing a full powder average Yoo01 ; the Mn30 is thus modelled as a giant spin with Ising-like anisotropy. The corresponding Hamiltonian is given by equation (1):

$$=DS_z^2+_{\mathrm{trans}}+g\mu _\mathrm{B}\mu _0\stackrel{}{S}\stackrel{}{H}$$
(1)
where $`D`$ is the axial anisotropy constant, $`S_z`$ is the $`z`$-component of the spin operator $`\stackrel{}{S}`$, $`_{\mathrm{trans}}`$ is the transverse anisotropy, containing $`S_x`$ and $`S_y`$ operators, $`g`$ is the electronic g-factor, $`\mu _\mathrm{B}`$ is the Bohr magneton, $`\mu _0`$ is the vacuum permeability, and $`\stackrel{}{H}`$ is the applied field. The last term in eq. 1 is the Zeeman energy associated with an applied field. The fit parameters were $`S=5`$, $`D`$ = -0.73 K, and $`g=2.00`$, with the transverse anisotropy $`_{\mathrm{trans}}`$ not included for simplicity (in fact, inclusion of a small transverse anisotropy term involving the transverse anisotropy constant, $`E`$, had no significant effect on the fit). These values suggest an upper limit to the potential energy barrier $`U`$ to magnetization relaxation (reversal) of $`U=S^2|D|=18`$ K, a value potentially large enough to make Mn30 a SMM at sufficiently low temperatures. This was therefore explored by hysteresis loop, relaxation, and hole digging measurements.

Studies were performed on single crystals using an array of micro-SQUIDs WW\_ACP\_01 . Hysteresis in magnetization versus magnetic field scans was observed, establishing Mn30 as a new SMM, the largest to date. Fig. 2 shows typical hysteresis loops, with the field applied approximately along the easy axis of magnetization. The blocking temperature is about 1.2 K, and the coercivity increases with decreasing temperature and increasing sweep rate, as expected for the superparamagnet-like behavior of a SMM. The loops do not display step-like features characteristic of resonant QTM between the energy states of the Mn30 molecule, as observed for several other SMMs Friedman96 ; Thomas96 ; Sangregorio97 ; Aubin98 ; Aubin98b ; Boskovic01 ; WW\_Science99 . However, it is possible that steps are present but smeared out by broadening effects due to a distribution of magnetization relaxation barriers (i.e. $`D`$ values), consistent with the distribution of Mn30 environments resulting from the disordered CMe3 groups and solvent molecules observed in the crystal structure. It is thus not possible to determine directly from the hysteresis loops whether resonant QTM is occurring in this large Mn30 molecule. Nevertheless, the first indication of quantum tunneling is the temperature independence of the hysteresis loops below about 0.3 K (Fig.2). Detailed measurements of the hysteresis loops show however that the loops are still time-dependent below 0.3 K. This is presented in the inset of Fig. 3 showing the temperature dependence of the coercive field at different field sweep rates. The strong time-dependence in the temperature-independent regime is another important indication of QTM.

QTM can also be confirmed with magnetization decay studies. The magnetization was first saturated at 5 K with a large applied field, the temperature decreased to a chosen value, and then the field removed and the magnetization decay monitored with time. This provided magnetization relaxation rates $`\mathrm{\Gamma }`$ at different temperatures allowing us to make an Arrhenius plot (Fig. 3) that is based on the Arrhenius law:

$$\mathrm{\Gamma }=\mathrm{\Gamma }_0exp(U_{\mathrm{eff}}/k_\mathrm{B}T)$$
(2)
where $`\mathrm{\Gamma }_0`$ is the pre-exponential factor, $`U_{\mathrm{eff}}`$ is the mean effective barrier to relaxation, and $`k_\mathrm{B}`$ is the Boltzmann constant. The slope in the thermally activated region above 0.5 K gave $`U_{\mathrm{eff}}/k_\mathrm{B}`$ = 15 K and $`\mathrm{\Gamma }_0=1.5\times 10^8`$ s-1. The mean barrier $`U_{\mathrm{eff}}`$ is smaller than the calculated $`U=S^2|D|=18`$ K, as expected for QTM between higher energy levels $`M_s`$ of the $`S=5`$ manifold. More importantly, the relaxation time becomes temperature-independent below 0.2 K at about 10-6 s-1, diagnostic of ground state QTM between the lowest energy $`m=\pm 5`$ levels. The crossover temperature between thermally activated relaxation and ground state tunneling is between 0.2 and 0.3 K. The ground state tunneling at $`<`$ 0.2 K rationalizes the appearance at these temperatures of a QTM step at H = 0 in the loops of Fig. 2. Since the tunneling is now only between the lowest $`m`$ levels, the distribution of Mn30 environments no longer has such a large broadening effect on the step at zero field, and a better defined step is seen.

Additional confirmation of QTM in Mn30 was obtained from the ‘quantum hole-digging’ method Prokofev98 ; WW\_PRL99 ; WW\_PRL00 ; Fernandez01 ; Tupitsyn03 . This is a relatively new method that can, among other things, establish whether resonant tunneling occurs even when steps are absent in hysteresis loops due to a distribution of energy barriers. The method is based on the simple idea that after a rapid field change, the resulting magnetization relaxation at short time periods is directly related to the number of molecules in resonance at the applied field; Prokof’ev and Stamp proposed Prokofev98  that this short time relaxation should follow a $`\sqrt{t}`$ ($`t`$ = time) relaxation law. Thus, the magnetization of the Mn30 molecules in the crystal was first saturated with a large negative field, and then a ‘digging field’ $`H_{\mathrm{dig}}`$ of 0.275 T was applied at 0.04 K for a chosen ‘digging time’ $`t_{\mathrm{dig}}`$. If QTM can occur in Mn30, then that fraction (and only that fraction) of the molecules that is in resonance at $`H_{\mathrm{dig}}`$ can undergo magnetization tunneling. After $`t_{\mathrm{dig}}`$, a field $`H_{\mathrm{probe}}`$ is applied and the magnetization relaxation rate is measured for short time periods; from this is calculated the short-time relaxation rate $`\mathrm{\Gamma }_{\mathrm{sqrt}}`$, which is related to the number of Mn30 molecules still available for QTM WW\_ACP\_01 . The entire procedure is then repeated at other $`H_{\mathrm{probe}}`$ fields. The resulting plot of $`\mathrm{\Gamma }_{\mathrm{sqrt}}`$ versus $`H_{\mathrm{probe}}`$ reflects the distribution of spins still available for tunneling after $`t_{\mathrm{dig}}`$, and it will display a ‘hole’ if resonant QTM did indeed occur during application of $`H_{\mathrm{dig}}`$ for time $`t_{\mathrm{dig}}`$. The obtained plot for Mn30 (Fig. 4) does show a hole in this distribution, corresponding to a depletion of those spins that were in resonance and able to tunnel during $`t_{\mathrm{dig}}`$ at $`H_{\mathrm{dig}}`$. The occurrence of resonant QTM in Mn30 is thus confirmed, since a classical system or a quantum system with an continuum of levels will not show a ‘quantum hole’ but rather a depletion of all spins with low barriers (Fig. 4).

The above results unambiguously demonstrate that Mn30 undergoes tunneling, involving the coherent reversal of about 600 interacting electron spins within the Mn/O magnetic core. Of the three diagnostic tests for resonant QTM, Mn30 clearly demonstrates two of them, temperature-independent relaxation and quantum hole digging. The third is steps in hysteresis loops, but these are broadened beyond resolution by the distribution of barriers ($`D`$ values) resulting from a distribution of Mn30 environments. Mn30 thus represents the largest SMM by far to unequivocally demonstrate QTM, establishing that quantum effects can be clearly observed, and studied, even in very large ’magnetic particles’ with dimensions approaching the mesoscale. This is made possible in Mn30 by its monodisperse nature and highly ordered arrangement within the crystal (notwithstanding the solvent and CMe3 disorder, which in an absolute sense is a small perturbation), which prevent complications from distributions of particle size, shape, surface roughness, and spin. The latter have severely hampered previous attempts to demonstrate QTM in large magnetic particles, both for classical magnetic materials such as Co metal and molecules such as the ferritin protein. The results have consequently usually been negative, unclear or controversial Awschalom92 ; Gider95 ; Tejada96 ; Mamiya02 ; WW\_PRL97\_BaFeO . In contrast, the present work demonstrates that crystalline assemblies of monodisperse molecular nanomagnets allow clear observation of quantum properties under straightforward conditions even for large systems approaching the mesoscale, the hydrocarbon shell of organic groups in Mn30 ensuring the Mn/O cores are exactly identical in every molecule.

From a broader point of view, the present work confirms that large SMMs can be prepared and studied as crystalline assemblies of monodisperse particles. Mn30 is much bigger than other currently known SMMs such as Mn4, Mn12, etc, and thus represents proof-of-feasibility of extending to larger species the many advantages of SMMs over classical magnetic particles. These include monodispersity, crystalline order, room-temperature solution synthesis, true solubility (rather than colloidal suspension) in common organic solvents, and an insulating shell of organic groups around the magnetic core, which additionally can be varied at will using standard synthetic chemistry methods. Such factors are clear benefits for the many proposed applications of SMMs, for example as qubits in quantum computing Leuenberger01 , where the larger is an individual SMM, the easier it should be to manipulate and measure its quantum state. The synthetic variation of the hydrocarbon shell can also provide a facile way to increase (or decrease, if desired WW\_Nature02 ) the isolation of a single SMM from its neighbors, and perhaps thus also help decrease decoherence problems.

We thank the USA National Science Foundation for support.


format: text

text:
# Charge symmetry breaking in 𝒏𝒑→𝒅𝝅^𝟎

## I The Experiment

The experiment was performed at TRIUMF with a 279.5 MeV neutron beam, a liquid hydrogen target, and the SASP magnetic spectrometer Walden and et al. (1999) positioned at $`0^{}`$. With these near threshold kinematics and the large acceptance of SASP, the full deuteron distribution from $`npd\pi ^0`$ was detected in one setting of the spectrometer thereby eliminating many systematic uncertainties. These deuterons form a distinct kinematic locus in momentum vs laboratory scattering angle, which is shown in fig. 1 for the collected data.

The TRIUMF CHARGEX facility Helmer and et al. (1987) produced the neutron beam by passing a high intensity proton beam through a thin $`{}_{}{}^{7}\mathrm{Li}`$ target. A sweeping magnet deflected the primary proton beam into a well- shielded dump. The liquid hydrogen target ($`\mathrm{LH}_2`$) was centered 92 cm downstream from the $`{}_{}{}^{7}\mathrm{Li}`$ target and was contained within a flat cylindrical volume, 10 cm in diameter with a nominal thickness of 2 cm. Two sets of veto counters (FEV1, FEV2) and a trigger counter set (FET) were each composed of a pair of plastic scintillators positioned above one another. This allowed more stable operation in the high (few MHz) particle rate environment. The thick veto scintillators were upstream of the $`\mathrm{LH}_2`$ and shadowed it. The FET counters were positioned immediately downstream of the $`\mathrm{LH}_2`$.

Three multi-wire proportional chambers, positioned upstream of the SASP entrance (FECs i.e. Front-End Chambers), provided tracking information for charged particles. Each FEC consisted of a pair of orthogonal wire planes. The first and last FECs, separated by 33 cm, were mounted to measure vertical and horizontal coordinates. The third FEC was positioned midway between the other two and rotated $`40^{}`$ with respect to them for efficiency measurements and to aid in multi-hit track reconstruction. Particle tracking near the SASP focal plane was provided by two vertical-drift chambers (VDCs). Three sets of scintillators, downstream from the VDCs, provided timing and particle ID information as well as sufficient redundancy to determine the efficiencies of all focal plane area detectors.

Measurements of $`np`$ elastic scattering with incident neutron beams that filled the same target space and produced protons that spanned the momentum distribution of the $`npd\pi ^0`$  reaction provided a stringent test of the description of the spectrometer acceptance. Further details on the apparatus and other technical aspects of the measurement are found in reference Hutcheon and et al. (2001).

## II Extraction of $`𝐀_{\mathrm{𝐟𝐛}}`$

Close to threshold, the $`npd\pi ^0`$  cross-section in the center-of-mass frame is given by

$$\frac{d\sigma }{d\mathrm{\Omega }}(\theta )=A_0+A_1P_1(\mathrm{cos}\theta )+A_2P_2(\mathrm{cos}\theta ),$$
(2)
where $`P_1`$ and $`P_2`$ are Legendre polynomials. The $`A_0`$ and $`A_2`$ coefficients were previously measured Hutcheon and et al. (1991) at a number of energies within 10 MeV above threshold. The presence of charge symmetry breaking is reflected in the $`A_1`$ term as it is odd in $`\mathrm{cos}\theta `$. In this standard parametrization, the angle integrated form of $`A_{\mathrm{fb}}`$ is given by $`A_{\mathrm{fb}}=\frac{1}{2}A_1/A_0`$.

For a given beam energy, $`\mathrm{cos}\theta `$ varies linearly with the longitudinal component of deuteron momentum in the laboratory reference frame. Ideally, the $`\mathrm{cos}\theta `$ distribution would be found by a suitable, simple projection of the data of fig. 1. However, the measured deuteron locus is distorted by energy loss, multiple scattering, energy spread of the beam, and spectrometer acceptance making a direct extraction of $`A_{\mathrm{fb}}`$ impossible. Instead, the data were binned according to laboratory momentum and angle (as in fig. 1) and compared to a model which represented the background due to $`\mathrm{C}(\mathrm{n},\mathrm{d})`$ reactions as a low-order polynomial and generated the locus of $`\mathrm{H}(\mathrm{n},\mathrm{d})\pi ^0`$ events by Monte Carlo simulation of the beam, target, reaction cross section, spectrometer and detectors.

The simulation was based on GEANT3. It began with a proton beam incident on the $`{}_{}{}^{7}\mathrm{Li}`$ target and included energy loss by the proton beam as well as the angular and energy distribution of neutrons from the $`{}_{}{}^{7}\mathrm{Li}(\mathrm{p},\mathrm{n})`$ reaction. Production of deuterons according to the distribution of equation 2 was allowed in the $`\mathrm{LH}_2`$ target and other hydrogenous material such as scintillators and their wrapping. Standard GEANT tracking options were adopted for deuteron energy loss and multiple scattering but the reaction losses, which amount to 1–2% and are momentum dependent, were parametrized from data on deuteron elastic and reaction cross-sections from hydrogen and carbon Auce and et al. (1987). Tracking through the SASP dipole used a field map obtained at 875 Amp and scaled up to the operating current of 905 Amp. Data were acquired in 10 different periods spanning two years and the simulation accounted for measured detector efficiencies, scintillator thresholds, missing FEC wires, and known changes in target thickness in a manner consistent with the actual running periods.

To reduce the possibility of psychological bias in matching simulation to data, a blind analysis technique was used which incorporated a hidden offset to the $`A_1/A_0`$ asymmetry parameter of the $`npd\pi ^0`$  generator. The collaborators developing the simulation and extracting the observable did not know the value of the offset until all consistency checks had been satisfied.

## III Systematic Effects

The acceptance of SASP is a function of the initial target position and direction of the deuteron as well as its momentum. Non–uniformities in the momentum acceptance of SASP would systematically produce a false asymmetry and had to be limited. High-statistics data from $`np`$ elastic scattering were collected and compared to model simulations to determine a fiducial volume of uniform acceptance. For these calibration measurements, the SASP magnets were set to their values for the $`npd\pi ^0`$  running, but the primary beam energy was adjusted so that the elastically scattered protons had a momentum deviation $`\delta =(pp_0)/p_0=4,0`$ or $`+4\%`$ compared to the central momentum of the deuterons of interest. Projections of the $`np`$ elastic data direction for position slices were formed, and the ratios of yields at $``$4% vs +4% were formed for both data and simulation; see fig. 2. The analysis software acceptance cuts in position and direction were then limited to the regions common to both data and simulation which were uniform in momentum to the statistical precision of the data.

Simulation vs simulation comparisons were carried out to determine how strongly experimental parameters and other effects were correlated with $`A_1/A_0`$. For example, momentum dependent deuteron reaction losses and detection efficiencies are obvious mechanisms which can mimic the effect of a non–zero $`A_1/A_0`$. Combining each correlation with the independently-determined uncertainty of its parameter gave the systematic contributions shown in Table 1. However, for the $`\mathrm{LH}_2`$ target thickness, the proton beam energy ($`T_{\mathrm{beam}}`$) and the central momentum of SASP ($`p_0`$) the independent information was not a sufficient constraint. Therefore, these three parameters, along with $`A_1/A_0`$, were treated as free parameters and their values extracted from fitting the data. To this end, simulations were made and $`\chi ^2`$ calculated for 81 points in a four-dimensional space, in which each of the four free parameters was stepped above and below a nominal value. $`\chi ^2`$ minimization techniques Press et al.  were then used to obtain the values of the parameters at the global $`\chi ^2`$ minimum, while the local curvature of the $`\chi ^2`$ surface gave their errors and mutual correlations.

## IV Results and Discussion

As a test of the model, the $`\chi ^2`$ calculations and fitting were repeated on subsets of the data and simulated data, selected according to whether the reaction occurred in the Top or Bottom part of the $`\mathrm{LH}_2`$ target. A second test divided events into those originating in the Left or Right part of the target. The best fit values and errors of $`A_1/A_0`$ (after removal of the offset) and the other three parameters are presented in Table 2.

The root mean square (rms) systematic error for the full acceptance and the four subspaces is $`2.7\%`$ with the standard binning scheme of 50 bins in $`\delta `$ and 20 bins in $`\theta `$, indicating a substantial discrepancy between data and the simulation. Pixel by pixel examination of the contribution to $`\chi ^2`$ revealed a systematic difference in the profile of the locus along lines of steepest ascent. The sign of the differences tended to be positive at the peak and negative at both the “inner” and “outer” margins of the locus, possibly due to inadequate treatment of deuteron scattering in the simulation.

A change in $`A_1/A_0`$ will not change the ratio of counts in peak vs margins of the locus because it multiplies $`\mathrm{cos}(\theta _{cm})`$. In contrast, the $`\mathrm{LH}_2`$ thickness, $`p_0`$, and $`\mathrm{T}_{\mathrm{beam}}`$ all shift or broaden the locus and thus are sensitive to the ratio of locus counts at the peak vs margins. It is reasonable to expect further rebinning to remove sensitivity to unimportant details of the simulation without losing sensitivity to $`A_1/A_0`$. We repeated the $`\chi ^2`$ grid search using 20 bins in $`\delta `$ and 10 bins in $`\theta `$, and again with 10 bins in $`\delta `$ and 5 bins in $`\theta `$. As expected, the fractional error dropped to 2.1% and 1.4%, respectively, with $`A_1/A_0`$ remaining consistent within errors. A more sophisticated binning scheme which treated the locus as a set of “elliptical” and “radial” bins on top of rectangular background bins produced an rms error of 2.1% for 2500 background bins and 36 locus bins, and an rms error of 1.2% for 100 background bins and 6 locus bins. In all fits and binning schemes the best fit values of the asymmetry in the acceptance subspaces agreed within errors with the value for the full acceptance, which is $`(34.4\pm 16)\times 10^4`$, implying $`A_{\mathrm{fb}}=[17.2\pm 8(\mathrm{stat})\pm 5.5(\mathrm{sys})]\times 10^4`$.

Theoretical predictions of $`A_{\mathrm{fb}}`$ have been made by Niskanen Niskanen (1999) using a meson-exchange coupled-channel model which showed that the major contribution by far is due to $`\pi \eta `$ (and $`\pi \eta ^{}`$) mixing in both the exchange and produced (outgoing) meson. At our energy the prediction is $`A_{\mathrm{fb}}=28\times 10^4`$, when accepted values are used for the $`\eta NN`$ coupling constant and the $`\pi ^0\eta `$ mixing matrix element ($`g_{\eta NN}^2/4\pi =3.68`$ from meson exchange $`NN`$ potential models Dumbrajs and et al. (1983) and $`\pi ^0||\eta =0.0059`$ $`\mathrm{GeV}^2`$ from analysis of $`\eta `$ decay data Coon and et al. (1986)). More recently, $`A_{\mathrm{fb}}`$ was revisited van Kolck et al. (2000) within the framework of chiral effective field theory where the issue of charge symmetery breaking in the rescattering amplitude of the exchanged pion was addressed. The resulting additional contribution to $`A_{\mathrm{fb}}`$ is then expressed in terms of two parameters $`\delta m_N`$ and $`\overline{\delta }m_N`$ representing contributions from the $`ud`$ quark mass difference and from electromagnetic effects within the nucleon, respectively. Specifically, at our energy, $`A_{\mathrm{fb}}`$ is expressed as

$$A_{\mathrm{fb}}=0.28\%\times \left[\left(\frac{g_{\eta NN}}{\sqrt{4\pi (3.68)}}\right)\left(\frac{\pi ^0||\eta }{0.0059\mathrm{GeV}^2}\right)\frac{0.87}{\mathrm{MeV}}(\delta m_N\frac{\overline{\delta }m_N}{2})\right]$$
(3)
where the first term arises from $`\pi \eta `$ mixing and the second from $`\pi ^0N`$ scattering. With the introduction of the new term, $`A_{\mathrm{fb}}`$ changes sign and becomes positive with an estimated upper value around $`+69\times 10^4`$, when large but reasonable values of $`\delta m_N`$ and $`\overline{\delta }m_N`$ are used van Kolck et al. (2000). Our positive experimental result strongly suggests, therefore, that such isospin violating $`\pi ^0N`$ interactions as outlined in reference van Kolck et al. (2000) are indeed significant.

The parameters $`\delta m_N`$ and $`\overline{\delta }m_N`$ are also constrained by the proton-neutron mass difference as

$$\mathrm{\Delta }_N=m_nm_p=\delta m_N+\overline{\delta }m_N=1.29\mathrm{MeV}.$$
(4)
When our $`A_{\mathrm{fb}}`$ result is combined with equations 3 and 4, and the values given above for the $`\eta NN`$ coupling constant and $`\pi \eta `$ mixing matrix element, we find that $`\delta m_N=1.66\pm 0.27`$ MeV and $`\overline{\delta }m_N=0.36\pm 0.27`$ MeV, assuming no theoretical uncertainties. We emphasize, however, that this last exercise is only meant to illustrate the significance and potential important implications of our $`A_{\mathrm{fb}}`$ result. Further theoretical studies are currently underway Fonseca et al.  to accommodate simultaneously the new CSB result of our study and that of a recent cross-section measurement of the isospin forbidden reaction $`dd\alpha \pi ^0`$ Stephenson .

###### Acknowledgements.

This work was supported by grants from The Natural Sciences and Engineering Research Council of Canada, The National Science Foundation, and The Ohio Supercomputer Center. TRIUMF is operated under a grant from the National Research Council of Canada.


text:
# Coulomb’s law modification in nonlinear and in noncommutative electrodynamics

## I INTRODUCTION

Renewed interest in non-linear electrodynamics (Born-Infeld theory) has originated in recent string theory investigations. This is primarily because the low energy dynamics of D-branes have been described by a nonlinear Born-Infeld type action. It is worth recalling at this stage that Born and Infeld  suggested to modify Maxwell’s electromagnetism so as to get rid of the divergencies of the theory such as the infinite self-energy of a point charge. The resulting theory is a nonlinear gauge theory endowed with interesting features, like finite electron self-energy and a regular point charge electric field at the origin. However, due to the nonlinearity, the corresponding field equations are very difficult to solve. In addition to the string interest, the Born-Infeld theory has also attracted considerable attention from different viewpoints. For example, in connection to duality symmetry, also in magnetic monopoles studies, in generation of multipole moments for charged particles, and possible experimental determination of the parameter that measures the nonlinearity of the theory. The advent of noncommutative field theories also called attention to nonlinear theories. In fact, recently there have been indications that the $`\theta `$-expanded version of the noncommutative U(1) gauge theory is equivalent to the expansion of the Born-Infeld action up to order $`F^3`$, which is something that we intend to check in the present work, for the specific case of the static potential between charges.

On the other hand, one may improve our understanding of gauge theories through a proper study of the concepts of screening and confinement. In this respect the interaction energy of an infinitely heavy quark-antiquark pair is a key tool which plays an important role in the understanding of quark confinement. Moreover, the static potential is an essential concept both in electrodynamics and in gravitation. The importance of the static potential is further shown, for example, in the description of non-relativistic bound systems like quarkonia, as well as in the definition of the lattice coupling. In actual calculations it is obtained most directly when a correct separation between gauge-invariant and gauge dependent degrees of freedom is made. Previously, we proposed a general framework for studying the confining and screening nature of the potential in gauge theories in terms of gauge-invariant but path-dependent field variables . An important feature of this methodology is that it provides a physically-based alternative to the usual Wilson loop approach. In this paper, we will examine another aspect of nonlinear theories, namely, the lowest-order modification of the static potential due to the presence of Born-Infeld type terms, and in noncommutative electrodynamics. We address this issue along the lines of reference .

In Sec. II we will calculate the lowest-order correction to the Coulomb energy of a fermion-antifermion system, for both Born-Infeld electrodynamics and for the $`\theta `$-expanded version of noncommutative U(1) gauge theory. Our calculations show that the static interaction of fermions in Abelian gauge theories is determined by the geometrical condition of gauge invariance.

## II Interaction Energy

### A Two-dimensional Born-Infeld Electrodynamics

As already stated, our principal purpose is to calculate explicitly the interaction energy between static pointlike sources for Born-Infeld electrodynamics. To this end we will calculate the expectation value of the energy operator $`H`$ in the physical state $`|\mathrm{\Phi }`$, which we will denote by $`H_\mathrm{\Phi }`$. However, before going into the four-dimensional Born-Infeld electrodynamics, we would like to first consider the two-dimensional case. This would not only provide the theoretical setup for our subsequent work, but also fix the notation. The starting point is the two-dimensional space-time Lagrangian:

$$=\beta ^2\left\{1\sqrt{1+\frac{1}{2\beta ^2}F_{\mu \nu }F^{\mu \nu }}\right\}A_0J^0,$$
(1)
where $`J^0`$ is the external current. The parameter $`\beta `$ measures the nonlinearity of the theory and in the limit $`\beta \mathrm{}`$ the Lagrangian (1) reduces to the Maxwell theory. In order to handle the square root in the Lagrangian (1) we introduce an auxiliary field $`v`$, such that its equation of motion gives back the original theory. This allows us to write the Lagrangian as

$$=\beta ^2\left\{1\frac{v}{2}\left(1+\frac{1}{2\beta ^2}F_{\mu \nu }F^{\mu \nu }\right)\frac{1}{2v}\right\}A_0J^0.$$
(2)
Once this is done, the canonical quantization of this theory from the Hamiltonian analysis point of view is straightforward and follows closely that of references . The canonical momenta read $`\mathrm{\Pi }^\mu =vF^{0\mu }`$, and one immediately identifies the two primary constraints $`\mathrm{\Pi }^0=0`$ and $`p\frac{}{v}=0`$. The canonical Hamiltonian is thus

$$H_C=𝑑x\left\{\beta ^2+\mathrm{\Pi }_1^1A^0+\frac{1}{2v}\left(\beta ^2\mathrm{\Pi }_1\mathrm{\Pi }^1\right)+\frac{\beta ^2}{2}v+A_0J^0\right\}.$$
(3)
The consistency condition $`\dot{\mathrm{\Pi }}_0=0`$ leads to the secondary constraint $`\mathrm{\Gamma }_1\left(x\right)_1\mathrm{\Pi }^1J^0=0`$. The consistency condition for the $`p`$ constraint yields no further constraints and just determines the field $`v`$,

$$v=\sqrt{1\frac{1}{\beta ^2}\mathrm{\Pi }_1\mathrm{\Pi }^1},$$
(4)
which will be used to eliminate $`v`$. The extended Hamiltonian that generates translations in time then reads $`H=H_C+𝑑x\left(c_0(x)\mathrm{\Pi }_0(x)+c_1(x)\mathrm{\Gamma }_1(x)\right)`$, where $`c_0(x)`$ and $`c_1(x)`$ are the Lagrange multipliers. Since $`\mathrm{\Pi }_0=0`$ for all time and $`\dot{A}_0\left(x\right)=[A_0\left(x\right),H]=c_0\left(x\right)`$, which is completely arbitrary, we discard $`A_0\left(x\right)`$ and $`\mathrm{\Pi }_0\left(x\right)`$ because they add nothing to the description of the system. Then, the Hamiltonian takes the form

$$H=𝑑x\left\{\beta ^2\left(\sqrt{1\frac{1}{\beta ^2}\mathrm{\Pi }_1\mathrm{\Pi }^1}1\right)c^{}\left(x\right)\left(_1\mathrm{\Pi }^1J^0\right)\right\},$$
(5)
where $`c^{}\left(x\right)=c_1\left(x\right)A_0\left(x\right)`$.

The quantization of the theory requires the removal of non-physical variables, which is done by imposing a gauge condition such that the full set of constraints becomes second class. A convenient choice is found to be

$$\mathrm{\Gamma }_2\left(x\right)\underset{C_{\xi x}}{}𝑑z^\nu A_\nu \left(z\right)\underset{0}{\overset{1}{}}𝑑\lambda x^1A_1\left(\lambda x\right)=0,$$
(6)
where $`\lambda `$ $`(0\lambda 1)`$ is the parameter describing the spacelike straight path $`x^1=\xi ^1+\lambda \left(x\xi \right)^1`$, and $`\xi `$ is a fixed point (reference point). There is no essential loss of generality if we restrict our considerations to $`\xi ^1=0`$. With this choice the nontrivial Dirac bracket is given by

$$\{A_1\left(x\right),\mathrm{\Pi }^1\left(y\right)\}^{}=\delta ^{\left(1\right)}\left(xy\right)_1^x\underset{0}{\overset{1}{}}𝑑\lambda x^1\delta ^{\left(1\right)}\left(\lambda xy\right).$$
(7)
Having outlined the necessary aspects of quantization, we now turn to the problem of obtaining the interaction energy between pointlike sources in Born-Infeld theory, where a fermion is localized at the origin $`\mathrm{𝟎}`$ and an antifermion at $`𝐲`$. As we have already indicated, we will calculate the expectation value of the energy operator $`H`$ in the physical state $`|\mathrm{\Phi }`$. From (5) we then get for the expectation value

$$H_\mathrm{\Phi }=\mathrm{\Phi }\left|𝑑x\left\{\beta ^2\left(\sqrt{1\frac{1}{\beta ^2}\mathrm{\Pi }_1\mathrm{\Pi }^1}1\right)\right\}\right|\mathrm{\Phi }.$$
(8)
As remarked by Dirac, the physical state can be written as

$$|\mathrm{\Phi }|\overline{\mathrm{\Psi }}\left(𝐲\right)\mathrm{\Psi }\left(\mathrm{𝟎}\right)=\overline{\psi }\left(𝐲\right)\mathrm{exp}\left(ie\underset{\mathrm{𝟎}}{\overset{𝐲}{}}𝑑z^iA_i\left(z\right)\right)\psi \left(\mathrm{𝟎}\right)|0,$$
(9)
where $`|0`$ is the physical vacuum. As before, the line integral appearing in the above expression is along a spacelike path starting at $`\mathrm{𝟎}`$ and ending at $`𝐲`$, on a fixed time slice. It is worth noting here that the strings between fermions have been introduced in order to have a gauge-invariant function $`|\mathrm{\Phi }`$. In other terms, each of these states represents a fermion-antifermion pair surrounded by a cloud of gauge fields sufficient to maintain gauge invariance.

Since we are interested in estimating the lowest-order correction to the Coulomb energy, we will retain only the leading quadratic term in the expression (8). Thus the expectation value simplifies to

$$H_\mathrm{\Phi }=\mathrm{\Phi }\left|𝑑x\left\{\frac{1}{2}\left(\mathrm{\Pi }_1\right)^2\frac{1}{8\beta ^2}\left(\mathrm{\Pi }_1\right)^4\right\}\right|\mathrm{\Phi }.$$
(10)
It is easy to see that the first term inside the curly bracket comes from the usual Maxwell theory while the second one is a correction which comes from the Born-Infeld modification. From our above Hamiltonian analysis we observe that

$$\mathrm{\Pi }_1\left(x\right)|\overline{\mathrm{\Psi }}\left(y\right)\mathrm{\Psi }\left(0\right)=\overline{\mathrm{\Psi }}\left(y\right)\mathrm{\Psi }\left(0\right)\mathrm{\Pi }_1\left(x\right)|0e_0^y𝑑z_1\delta ^{\left(1\right)}\left(z_1x\right)|\mathrm{\Phi }.$$
(11)
Substituting this back into (10), we obtain

$$H_\mathrm{\Phi }=H_0+\frac{e^2}{2}𝑑x\left(_0^y𝑑z_1\delta ^{\left(1\right)}\left(z_1x\right)\right)^2\frac{e^4}{8\beta ^2}𝑑x\left(_0^y𝑑z_1\delta ^{\left(1\right)}\left(z_1x\right)\right)^4,$$
(12)
where $`H_0=0\left|H\right|0`$. We further note that

$$\frac{e^2}{2}𝑑x\left(_0^y𝑑z\delta ^1\left(z_1x\right)\right)^2=\frac{e^2}{2}L,$$
(13)
with $`|y|L`$. By employing Eq. (13) we can reduce Eq. (12) to

$$V=\frac{e^2}{2}\left(1\frac{e^2}{4\beta ^2}\right)L.$$
(14)
Hence we see that the static interaction between fermions arises only because of the requirement that the $`|\overline{\mathrm{\Psi }}\mathrm{\Psi }`$ states be gauge invariant. The above result reveals that the effect of adding the Born-Infeld term is to decrease the energy. Nevertheless, the confining nature of the potential is preserved.

Eq. (14) exhibits the same formal structure as the one obtained for the massive Schwinger model. In fact, the same calculation for the massive Schwinger model  gives

$$V=\frac{q^2}{2}\left(1+\frac{e^2}{4\pi ^2m\mathrm{\Sigma }}\right)^1L,$$
(15)
where $`\mathrm{\Sigma }=\left(\frac{e}{2\pi ^{\frac{3}{2}}}\right)\mathrm{exp}\left(\gamma _E\right)`$ with $`\gamma _E`$ the Euler-Mascheroni constant, $`m`$ and $`e`$ are the mass and charge of the dynamical fermions. Here $`q`$ refers to the probe charges. Considering the limit $`me`$ and $`qe`$, we get

$$V=\frac{e^2}{2}\left(1\frac{e^2}{4\pi ^2m\mathrm{\Sigma }}\right)L.$$
(16)
Therefore, for this special case the massive Schwinger model simulates the features of the Born-Infeld theory. This analysis suggests the interesting possibility of identifying the parameter $`\beta `$ with the mass of the dynamical fermions. This, however, is a separate question and which we do not intend to address here.

Before concluding this subsection we discuss an alternative derivation of the result (14), which highlights certain distinctive features of our methodology. We start by considering

$$Ve\left(𝒜_0\left(\mathrm{𝟎}\right)𝒜_0\left(𝐲\right)\right),$$
(17)
where the physical scalar potential is given by

$$𝒜_0(x^0,x^1)=_0^1𝑑\lambda x^1E_1\left(\lambda x^1\right).$$
(18)
This follows from the vector gauge-invariant field expression

$$𝒜_\mu \left(x\right)A_\mu \left(x\right)+_\mu \left(_\xi ^x𝑑z^\mu A_\mu \left(z\right)\right),$$
(19)
where, as in Eq.(6), the line integral is along a spacelike path from the point $`\xi `$ to $`x`$, on a fixed time slice. The gauge-invariant variables (19) commute with the sole first constraint (Gauss’ law), confirming that these fields are physical variables . Note that Gauss’ law for the present theory reads

$$_1\frac{E^1}{\sqrt{1\frac{1}{\beta ^2}\left(E^1\right)^2}}=J^0,$$
(20)
where $`E^1`$ is the one-dimensional electric field. For $`J^0(t,x)=e\delta ^{\left(1\right)}\left(x\right)`$, the electric field is given by

$$E^1=\frac{e}{2}\frac{1}{\sqrt{1+\frac{e^2}{4\beta ^2}}}\widehat{x}^1,$$
(21)
here $`\widehat{x}^1`$ is an unit vector ($`\widehat{x}^1=\frac{x^1}{|x^1|}`$). Finally, making use of (21) and (18) in (17), we find

$$V=\frac{e^2}{2}\left(1\frac{e^2}{4\beta ^2}\right)L,$$
(22)
where $`|y|L`$.

### B Four-dimensional Born-Infeld Electrodynamics

We now turn our attention to the calculation of the interaction energy between static pointlike sources for the four-dimensional Born-Infeld electrodynamics. In such a case the Lagrangian reads

$$=\beta ^2\left\{1\sqrt{1+\frac{1}{2\beta ^2}F^{\mu \nu }F_{\mu \nu }\frac{1}{16\beta ^4}\left(^{\mu \nu }F_{\mu \nu }\right)^2}\right\}A_0J^0.$$
(23)
where $`\frac{1}{4}F_{\mu \nu }F^{\mu \nu }=\frac{1}{2}\left(𝐄^2𝐁^2\right)`$ , $`F_{\mu \nu }^{\mu \nu }=4\left(𝐄𝐁\right)`$ and $`J^0`$ is the external current.

Before we proceed to work out explicitly the energy, we shall begin by summarizing the Hamiltonian analysis of the theory (23). Once again, we will introduce an auxiliary field $`v`$ to handle the square root in the Lagrangian (23). Expressed in terms of this field, the Lagrangian (23) takes the form

$$=\beta ^2\left\{1\frac{v}{2}\left(1+\frac{1}{2\beta ^2}F_{\mu \nu }F^{\mu \nu }\frac{1}{16\beta ^4}\left(^{\mu \nu }F_{\mu \nu }\right)^2\right)\frac{1}{2v}\right\}A^0J_0.$$
(24)
With this in hand, the canonical momenta are $`\mathrm{\Pi }^\mu =v\left(F^{0\mu }\frac{1}{4\beta ^2}F_{\alpha \beta }^{\alpha \beta }^{0\mu }\right)`$, and one immediately identifies the two primary constraints $`\mathrm{\Pi }^0=0`$ and $`p\frac{}{v}=0`$. The canonical Hamiltonian of the model can be worked out as usual and is given by the expression

$$H_C=d^3x\left\{\beta ^2+\mathrm{\Pi }_i^iA^0+\frac{1}{2v}\left(𝚷^2+\beta ^2\right)+\frac{v}{2}\left(𝐁^2+\beta ^2\right)\frac{1}{2v\beta ^2}\frac{\left(𝚷𝐁\right)^2}{\left(1+\frac{𝐁^2}{\beta ^2}\right)}+A_0J^0\right\}.$$
(25)
Requiring the primary constraint $`\mathrm{\Pi }_0`$ to be preserved in time yields the secondary constraint (Gauss’ law) $`\mathrm{\Gamma }_1\left(x\right)_i\mathrm{\Pi }^iJ^0=0`$. Similarly for the constraint $`p`$, we get the auxiliary field $`v`$ as

$$v=\frac{1}{\beta ^2\left(1+\frac{1}{\beta ^2}𝐁^2\right)}\sqrt{\beta ^2\left(𝚷^2+\beta ^2\right)\left(1+\frac{1}{\beta ^2}𝐁^2\right)\left(𝚷𝐁\right)^2}.$$
(26)
The extended Hamiltonian that generates translations in time then reads $`H=H_C+d^3x\left(c_0(x)\mathrm{\Pi }_0(x)+c_1(x)\mathrm{\Gamma }_1(x)\right)`$, where $`c_0(x)`$ and $`c_1(x)`$ are Lagrange multipliers. As before, neither $`A_0\left(x\right)`$ nor $`\mathrm{\Pi }_0\left(x\right)`$ are of interest in describing the system and may be discarded from the theory. Thus we are left with the following expression for the Hamiltonian

$$H=d^3x\left\{\sqrt{\beta ^2\left(𝚷^2+\beta ^2\right)\left(1+\frac{𝐁^2}{\beta ^2}\right)\left(𝚷𝐁\right)^2}\beta ^2+c^{}(x)\left(_i\mathrm{\Pi }^iJ^0\right)\right\},$$
(27)
where $`c^{}(x)=c_1(x)A_0(x)`$.

Since our main motivation is to compute the static potential for the Born-Infeld theory, we will adopt the same gauge-fixing condition that was used in the last subsection, that is,

$$\mathrm{\Gamma }_2\left(x\right)\underset{C_{\xi x}}{}𝑑z^\nu A_\nu \left(z\right)\underset{0}{\overset{1}{}}𝑑\lambda x^iA_i\left(\lambda x\right)=0.$$
(28)
Here again $`\lambda `$ $`(0\lambda 1)`$ is the parameter describing the spacelike straight path $`x^i=\xi ^i+\lambda \left(x\xi \right)^i`$ with $`i=1,2,3`$, and $`\xi `$ is a fixed (reference) point. There is no essential loss of generality if we restrict our considerations to $`\xi ^i=0`$. In this way, the fundamental Dirac bracket can be rewritten as

$$\{A_i\left(x\right),\mathrm{\Pi }^j\left(y\right)\}^{}=\delta _i^j\delta ^{\left(3\right)}\left(xy\right)_i^x\underset{0}{\overset{1}{}}𝑑\lambda x^i\delta ^{\left(3\right)}\left(\lambda xy\right).$$
(29)
We are now in a position to compute the potential energy for static charges in this theory. To do this, we will use the gauge-invariant scalar potential which is now given by

$$𝒜_0(t,𝐫)=_0^1𝑑\lambda r^iE_i(t,\lambda 𝐫).$$
(30)
It follows from the above discussion that Gauss’ law takes the form

$$_i\frac{E_i}{\sqrt{1\frac{𝐄^2}{\beta ^2}}}=J^0.$$
(31)
For $`J^0(t,𝐫)=e\delta ^{\left(3\right)}\left(𝐫\right)`$, the electric field follows as

$$E^i\left(𝐫\right)=\frac{e}{4\pi }\frac{1}{\sqrt{|𝐫|^4+\rho _0^2}}\widehat{r}^i,$$
(32)
where $`\rho _0\frac{e}{4\pi \beta }`$ and $`\widehat{r}^i=\frac{𝐫^i}{|𝐫|}`$. From this expression it should be clear that the electric field of a pointlike charge is regular at the origin, in contrast to the usual Maxwell theory. As a consequence, equation (30) becomes

$$𝒜_0(t,𝐫)=\frac{e}{4\pi }_0^1𝑑\lambda \frac{r}{\sqrt{\left(\lambda r\right)^4+\rho _0^2}}.$$
(33)
Again, as in the two-dimensional case, it is sufficient to retain the leading quadratic term in Eq. (33). Thus we obtain

$$𝒜_0(t,𝐫)=\frac{e}{4\pi r}_0^1𝑑\lambda \left\{\frac{1}{\lambda ^2}+\frac{1}{2}\frac{a^4}{\lambda ^6}\right\},$$
(34)
where $`a^4\frac{\rho _0^2}{r^4}=\frac{e^2}{16\beta ^2\pi ^2r^4}`$.

In terms of $`𝒜_0(t,𝐫)`$, the potential for a pair of static pointlike opposite charges located at $`\mathrm{𝟎}`$ and $`𝐋`$, that is, $`J^0(t,𝐫)=e\left\{\delta ^{\left(3\right)}\left(𝐫\right)\delta ^{\left(3\right)}\left(𝐫𝐋\right)\right\}`$, is given by

$$Ve\left(𝒜_0\left(\mathrm{𝟎}\right)𝒜_0\left(𝐋\right)\right)=\frac{e^2}{4\pi }\frac{1}{L}\left(1+\frac{e^2}{160\pi ^2\beta ^2}\frac{1}{L^4}\right),$$
(35)
with $`L=|𝐋|`$.

This result shows the usual Coulomb potential with a long-range correction due to the second term of the form

$$|\mathrm{\Delta }V|=\frac{e^4}{640\pi ^3\beta ^2}\frac{1}{L^5}.$$
(36)
To $`𝒪\left(\frac{1}{\beta ^2}\right)`$ Born-Infeld electrodynamics displays a marked qualitative departure from the usual Maxwell theory. In Fig.$`1`$ we show the effect of the $`\beta `$ -correction, for the case $`\beta =10`$. Note that we have plotted the logarithm of $`V(L)`$ as a function of $`L`$. It is important to notice that the presence of the second term on the right-hand side of Eq. (35), which dominates for small $`L`$ values, causes $`V`$ to decrease. Thus, from a physical point of view, this discussion allows us to say that the effect of adding the Born-Infeld term is to generate more stable bound states of charged particles. Let us also mention here that if we had considered the value of $`\beta `$ as predicted in , the correction to the Coulomb potential would have been negligible. Hence it becomes important to obtain an independent lower bound for $`\beta `$ from our long-range correction, which we hope to report elsewhere.

### C Non-commutative electrodynamics

We now want to extend what we have done to non-commutative electrodynamics to leading order in $`\theta `$. As before, we will concentrate on the effect of including the $`\theta `$ term in the static potential. In such a case the Lagrangian  reads

$$=\frac{1}{4}F_{\mu \nu }^2+\frac{1}{8}\theta ^{\alpha \beta }F_{\alpha \beta }F_{\mu \nu }^2\frac{1}{2}\theta ^{\alpha \beta }F_{\mu \alpha }F_{\nu \beta }F^{\mu \nu }A_0J^0,$$
(37)
( $`\mu ,\nu ,\alpha ,\beta =0,1,2,3`$ ). Here, the field-strength tensor is expressed in terms of potentials in the usual way, $`F_{\mu \nu }=_\mu A_\nu _\nu A_\mu `$, and $`J^0`$ is an external current. As is well known $`\theta ^{\alpha \beta }`$ is a real constant antisymmetric tensor, and from now on we take $`\theta ^{0\alpha }=0`$ and $`\theta ^{ij}=\epsilon ^{ijk}\theta ^k`$.

The above Lagrangian will be the starting point of the Dirac constrained analysis . The canonical momenta following from Eq. (37) are $`\mathrm{\Pi }^\mu =\left(1\frac{1}{2}\theta ^{\alpha \beta }F_{\alpha \beta }\right)F^{\mu 0}\theta ^{\mu \beta }F_{\nu \beta }F^{0\nu }\theta _{\alpha \beta }F^{0\alpha }F^{\mu \beta }`$, which results in the usual primary constraint $`\mathrm{\Pi }^0=0`$ and $`\mathrm{\Pi }^i=\left(1\frac{1}{2}\theta ^{kl}F_{kl}\right)F^{i0}\theta ^{il}F_{jl}F^{0j}\theta _{kl}F^{0k}F^{il}`$ ( $`i,j,k,l=1,2,3`$ ). Defining the electric and magnetic fields by $`E^i=F^{i0}`$ and $`B^i=\frac{1}{2}\epsilon ^{ijk}F_{jk}`$, respectively, the canonical Hamiltonian assumes the form

$$H_C=d^3x\left\{\frac{1}{2}\left(𝐄^2+𝐁^2\right)\left(1+𝜽𝐁\right)\left(𝜽𝐄\right)\left(𝐄𝐁\right)A_0\left(_i\mathrm{\Pi }^iJ^0\right)\right\}.$$
(38)
Time conservation of the primary constraint leads to the secondary constraint $`\mathrm{\Gamma }_1(x)_i\mathrm{\Pi }^iJ^0=0`$, and the time stability of the secondary constraint does not induce more constraints, which are first class. It should be noted that the constrained structure for the gauge field remains identical to the Born-Infeld theory. Thus, the quantization can be done in a similar manner to that in the previous subsection. In view of this situation, we pass now to the calculation of the interaction energy.

Following our earlier procedure, we will compute the expectation value of the noncommutative electrodynamics Hamiltonian in the physical state $`|\mathrm{\Phi }`$ (Eq. (9)). That is,

$$H_\mathrm{\Phi }=\mathrm{\Phi }\left|d^3x\left\{\frac{1}{2}\left(𝐄^2+𝐁^2\right)\left(1+𝜽𝐁\right)\left(𝜽𝐁\right)\left(𝐄𝐁\right)\right\}\right|\mathrm{\Phi }.$$
(39)
From our above Hamiltonian analysis, Eq. (39) can be simplified

$$H_\mathrm{\Phi }=\mathrm{\Phi }\left|d^3x\left\{\frac{1}{2}𝐄^2\left(1+𝜽𝐁\right)\left(𝜽𝐁\right)\left(𝐄𝐁\right)\right\}\right|\mathrm{\Phi }.$$
(40)
According to the definition of the canonical momenta $`\mathrm{\Pi }^i`$, we may also write

$$E^i=\left(1+𝜽𝐁\right)\mathrm{\Pi }^i\left(𝜽𝚷\right)B^i\theta ^i\left(𝚷𝐁\right),$$
(41)
to lowest order in $`\theta `$. Using Eq. (41) we can rewrite Eq. (40) in the following way

$$H_\mathrm{\Phi }=\mathrm{\Phi }\left|d^3x\left\{\frac{1}{2}𝚷^2+\frac{3}{2}\left(𝜽𝐁\right)𝚷^23\left(𝜽𝚷\right)\left(𝚷𝐁\right)\right\}\right|\mathrm{\Phi }.$$
(42)
Taking into account the preceding Hamiltonian structure, we first note that

$$\mathrm{\Pi }_i\left(𝐱\right)=\overline{\psi }\left(𝐲\right)\psi \left(\mathrm{𝟎}\right)\mathrm{\Pi }_i\left(𝐱\right)|0+e_\mathrm{𝟎}^𝐲𝑑z_i\delta ^{\left(3\right)}\left(𝐱𝐳\right)|\mathrm{\Phi }.$$
(43)
Using this in (42) we then evaluate the expectation value in the presence of the static charges

$$H_\mathrm{\Phi }=H_0+V_1+V_2^\theta +V_3^\theta ,$$
(44)
where $`H_0=0\left|H\right|0`$. The $`V_1,V_2^\theta ,V_3^\theta `$ terms are given by

$$V_1=\frac{e^2}{2}d^3x\left(_\mathrm{𝟎}^𝐲𝑑z_i\delta ^{\left(3\right)}\left(𝐳𝐱\right)\right)^2,$$
(45)
$$V_1^\theta =\frac{3}{2}e^2d^3x\left(𝜽𝐁\left(𝐱\right)\right)\left(_\mathrm{𝟎}^𝐲𝑑z_i\delta ^{\left(3\right)}\left(𝐳𝐱\right)\right)^2,$$
(46)
$$V_2^\theta =3e^2d^3x\left(\theta _iB_j\left(𝐱\right)\right)_\mathrm{𝟎}^𝐲𝑑z_i\delta ^{\left(3\right)}\left(𝐳𝐱\right)_\mathrm{𝟎}^𝐲𝑑z_j^{}\delta ^{\left(3\right)}\left(𝐳^{}𝐱\right).$$
(47)
The integrals over $`z_i`$ and $`z_j^{}`$ are zero except on the contour of integration.

Here we make the following observations. First, we note that the term (45) may look peculiar, but it is nothing but the familiar Coulomb interaction plus a self-energy term. In effect, as was explained in , by using spherical coordinates the integral $`_\mathrm{𝟎}^𝐲𝑑z_i\delta ^{\left(3\right)}\left(𝐳𝐱\right)`$ can also be written as

$$_\mathrm{𝟎}^𝐲𝑑z_i\delta ^3\left(𝐱𝐳\right)=\frac{𝐲_i}{|𝐲|}\frac{1}{|𝐲𝐱|^2}\underset{l,m}{}Y_{lm}^{}(\theta ^{},\varphi ^{})Y_{lm}(\theta ,\varphi ).$$
(48)
By means of (48) and using usual properties for the spherical harmonics, the term (45) reduces to the Coulomb energy after subtracting the self-energy term. On the other hand, it should be noted that in order to evaluate Eqs. (46) and (47) we need to know the magnetic field $`𝐁(𝐱)`$. However, since we are dealing with an external constant field $`\theta `$, we restrict ourselves to constant magnetic fields, $`𝐁(𝐱)=𝐁(\mathrm{𝟎})`$. Then, with this assumption and using (48), the interquark potential at lowest order in $`\theta `$ becomes

$$V=\frac{e^2}{4\pi }\frac{1}{L}\left[1+3\left(𝜽𝐁\left(\mathrm{𝟎}\right)\right)6\left(𝜽\widehat{𝐫}\right)\left(𝐁\left(\mathrm{𝟎}\right)\widehat{𝐫}\right)\right],$$
(49)
where $`L𝐲`$ and $`\widehat{𝐫}_i=\frac{𝐲_i}{|𝐲|}`$.

Accordingly, to lowest order in $`\theta `$, the nature of the static potential remains unchanged. Nevertheless, the introduction of the noncommutative parameter induces a renormalization of the charge, which is absent in the corresponding ordinary spacetime. In such a case, Eq. (49) may be rewritten as

$$V=\frac{e_{ren}^2}{4\pi }\frac{1}{L}.$$
(50)
To conclude, the expressions for the corrections to the static energy obtained from both the Born-Infeld and noncommutative electrodynamics (to order $`e^2`$) are quite different. This means that the two theories are not equivalent. Born-Infeld electrodynamics has a rich structure which is reflected in a long-range correction to the Coulomb potential, which is not present in its noncommutative counterpart. On the other hand, the present investigation reveals the general applicability of our methodology. It seems a challenging work to extend to higher orders the above analysis. We expect to report on progress along these lines soon.

## III ACKNOWLEDGMENTS

Work supported in part by Fondecyt (Chile) grant 1030355.


text:
# Biased Metropolis Sampling for Rugged Free Energy Landscapes

## 1 Introduction

Reliable simulations of biomolecules are one of nowadays grand challenges of computational science. In particular the problem of protein folding thermodynamics starting purely from an aminino-acid sequence has received major attention. Until recently the prevailing view has been that it is elusive to search for the native states with present day simulational techniques due to limitations of time scale and force field accuracy. Now a barrier appears to be broken, as it was reported Pande BBA5 (2002b) that large scale distributed computing allows to achieve folding of the 23-residue design mini-protein BBA5, which is relatively insensitive to inaccuracies of the force field. The relaxation dynamics of the computer simulation is found to be in good agreement with the experimentally observed folding times and equilibrium constants.

Molecular dynamics (MD) technics, for a review see Frenkel and Smit (1996a), tend to be the method of first choice for simulations of biomolecules. One of the attractive features of MD is that it allows to follow the physical time evolution of the system under investigation. Nevertheless, there has also been activity based on the Metropolis method Metropolis (1953b), which allows only for limited dynamical insights and has its strength in the generation of configurations which are in thermodynamical equilibrium. A major advantage of Metropolis simulations of biomolecules is that they allows for updates which are large moves when one has to follow dynamical trajectories. Such updates may help to overcome the kinetic trapping problem, which is due to a large number of local minima in the free energy space of typical biomolecules. Already Metropolis simulations of the canonical Gibbs-Boltzmann ensemble may jump certain barriers. Using generalized ensembles, which enlarge the Gibbs-Boltzmann ensemble and/or replace the canonical weights by other weighting factors, further progress has been made. For reviews see Okamoto review (2001b); Berg review (2002b).

To my knowledge umbrella sampling Torrie and Valleau (1977b) was the first generalized ensemble method of the literature. Its potential for applications to a wide range of interesting physical problems remained for a long time dormant. Apparently, one reason was that the computational scientists in various areas shied away from performing simulations with an a-priori unknown weighting factor. As Li and Scheraga put it Li and Scheraga (1988b): “The difficulty of finding such a weighting factor has prevented wide applications of the “umbrella sampling” method to many physical systems.” This changed with the introduction of the multicanonical approach Berg and Neuhaus (1992b) to complex systems Berg and Celik (1992b); Hansmann and Okamoto (1993b); Hao and Scheraga (1994b). Besides having the luck that the controversy surrounding its first application got quickly resolved by analytical results Borgs and Janke (1992b), the multicanonical approach addressed aggressively the problem of finding reliable weights. Nowadays a starter kit of Fortran programs is available on the Web Berg Fortran (2003b).

The replica exchange or parallel tempering (PT) method Geyer PT (1996b); Hukusima and Nemeto PT (1996b) uses a generalized ensemble which enlarges the canonical configuration space by allowing for the exchange of temperatures within a set of canonical ensembles. In this and other generalized ensembles the Metropolis dynamics at low temperatures can be accelerated by excursions to disordered configurations at higher temperatures. Note that, in contrast to multicanonical simulations, the probabilities of canonically rare configurations are not enhanced by PT. For the purpose of simulating biomolecules PT was, e.g., studied in Ref.Hansmann PT (1997b); Okamoto PT (1999b); Pande PT (2003b).

In this article we discuss a biased updating scheme Berg (2003b), which enhances already the dynamics of canonical Metropolis simulations and is easily integrated into generalized ensemble simulations too. The latter point is illustrated for parallel tempering. The biased updating scheme is inspired by the funnel picture of protein folding Funnel (1987b). At relatively high temperatures probability densities of the dynamical variables are estimated and used at lower temperatures to enhance the chance of update proposals to lie in the statistically relevant regions of the configuration space. In some sense this is nothing else but an elaboration of the original importance sampling concept of Metropolis et al.Metropolis (1953b).

This article is organized as follows. All-atom protein models, the funnel picture and the rugged Metropolis updating are introduced in the next section. In the subsequent section numerical results are presented for the brain peptide Met-Enkephalin. A short summary and conclusions are given in the final section.

## 2 All-Atom Protein Models

Proteins are linear polymers of the 20 naturally occurring amino acids. Small proteins are called peptides. The problem of protein folding is to predict (at room temperature and in solvents) the 3d conformations (native structures) from the sequence of amino acids. A conformational energy function models the interactions between the atoms (units kcal/mol):

$$E_{\mathrm{tot}}=E_{\mathrm{es}}+E_{vdW}+E_{\mathrm{hb}}+E_{\mathrm{tors}}+E_{\mathrm{sol}}.$$
(1)
The contributions to the energy function are (charges are in units of the elementary charge):

$`\mathrm{Electrostatic}E_{\mathrm{es}}`$ $`=`$ $`{\displaystyle \underset{ij}{}}{\displaystyle \frac{332q_iq_j}{ϵr_{ij}}},`$ (2)
$`\mathrm{van}\mathrm{der}\mathrm{Waals}E_{\mathrm{vdW}}`$ $`=`$ $`{\displaystyle \underset{ij}{}}\left({\displaystyle \frac{A_{ij}}{r_{ij}^{12}}}{\displaystyle \frac{B_{ij}}{r_{ij}^6}}\right),`$ (3)
$`\mathrm{hydrogen}\mathrm{bond}E_{\mathrm{hb}}`$ $`=`$ $`{\displaystyle \underset{ij}{}}\left({\displaystyle \frac{C_{ij}}{r_{ij}^{12}}}{\displaystyle \frac{D_{ij}}{r_{ij}^{10}}}\right),`$ (4)
$`\mathrm{torsion}E_{\mathrm{tors}}`$ $`=`$ $`{\displaystyle \underset{l}{}}U_l\left[1\pm \mathrm{cos}(n_l\alpha _l)\right].`$ (5)
The solvent accessible surface method allows for an approximative inclusion of solvent interactions:

$$E_{\mathrm{sol}}=\underset{i}{}\sigma _iA_i,$$
(6)
where $`\sigma _i`$ is the solvation parameter for atom $`i`$ and $`A_i`$ the conformation dependent solvent accessible surface area.

The $`r_{ij}`$ are the distances between the atoms and the $`\alpha _l`$ are the torsion angles for the chemical bonds. The parameters $`q_i,A_{ij},B_{ij},C_{ij},D_{ij},U_l`$ and $`n_l`$ are determined from crystal structures of amino acids and a number of thus obtained force fields are given in the literature. Bond length and angles fluctuate little and are normally set constant. The important degrees of freedom are the dihedral angles $`\varphi ,\psi ,\omega `$ and $`\chi `$, see Fig.1.

The funnel picture of Bryngelson and Wolyness Funnel (1987b) gives qualitative insight into the process of protein folding, see Fig.2, where a schematic sketch of the free energy versus a suitable reaction coordinate is given. However, there is no generical definition of a good reaction coordinate, as the funnel lives in the high-dimensional configuration space. Here we give a parameter free funnel description Berg (2003b) from higher to lower temperatures, which suggests a method for designing the a-priori Metropolis weights for simulations of biomolecules.

To be definite, we use the all-atom energy function ECEPP/2 (1987b) ECEPP/2 (Empirical Conformational Energy Program for Peptides). Our dynamical variables $`v_i`$ are then the dihedral angles, each chosen to be in the range $`\pi v_i<\pi `$, so that the volume of the configuration space is $`K=(2\pi )^n`$. Let us define the support of a pd of the dihedral angles. The support of a pd is the region of configuration space where the protein wants to be. Mathematically, we define $`K^p`$ to be the smallest sub-volume of the configuration space for which

$$p=_{K^p}\underset{i=1}{\overset{n}{}}dv_i\rho (v_1,\mathrm{},v_n;T)$$
(7)
holds. Here $`0<p<1`$ is a probability, which ought to be chosen close to one, e.g., $`p=0.95`$. The free energy landscape at temperature $`T`$ is called rugged, if the support of the pd consists of many disconnected parts (this depends of course a bit on the adapted values for $`p`$ and “many”). That a protein folds at room temperature, say $`300K`$, into a unique native structure $`v_1^0,\mathrm{},v_n^0`$ means that its pd $`\rho (v_1,\mathrm{},v_n;300K)`$ describes small fluctuation around this structure. We are now ready to formulate the funnel picture in terms of pds

$$\rho _r(v_1,\mathrm{},v_n)=\rho (v_1,\mathrm{},v_n;T_r),r=1,\mathrm{},s,$$
(8)
which are ordered by the temperatures $`T_r`$, namely

$$T_1>T_2>\mathrm{}>T_f.$$
(9)
The sequence (8) constitutes a protein funnel when, for a reasonable choice of the probability $`p`$ and the temperatures (9), the following holds:

1. The pds are rugged.
2. The support of a pd at lower temperature is contained in the support of a pd at higher temperature

$$K_1^pK_2^p\mathrm{}K_f^p,$$
(10)
e.g. for $`p=0.95`$, $`T_1=400K`$ and $`T_f=300K`$.
3. With decreasing temperatures $`T_r`$ the support $`K_r^p`$ shrinks towards small fluctuations around the native structure.

Properties 2 and 3 are fulfilled for many systems of statistical physics, when some groundstate stands in for the native structure. The remarkable point is that they may still hold for complex systems with a rugged free energy landscape, i.e., with property 1 added. In such systems one finds typically local free energy minima, which are of negligible statistical importance at low temperatures, while populated at higher temperatures. In simulations at low temperatures the problem of the canonical ensemble approach is that the updating tends to get stuck in those local minima. This prevents convergence towards the native structure on realistic simulation time scales. On the other hand, the simulations move quite freely at higher temperatures, where the native structure is of negligible statistical weight. Nevertheless, the support of a protein pd may already be severely restricted, as we shall illustrate. The idea is to use a relatively easily calculable pd at a higher temperature to improve the performance of the simulation at a lower temperature.

The Metropolis importance sampling would be perfected, if we could propose new configurations $`\{v_i^{}\}`$ with their canonical pd $`\rho (v_1^{},\mathrm{},v_n^{};T)`$. Due to the funnel property 2 we expect that an estimate $`\overline{\rho }(v_1,\mathrm{},v_n;T^{})`$ from some sufficiently close-by higher temperature $`T^{}>T`$ will feed useful information into the simulation at temperature $`T`$. The potential for computational gains is large because of the funnel property 3. The suggested scheme for the Metropolis updating at temperature $`T_r`$ is to propose new configurations $`\{v_i^{}\}`$ with the pd $`\overline{\rho }_{r1}(v_1^{},\mathrm{},v_n^{})`$ and to accept them with the probability

$$P_a=\mathrm{min}[1,\mathrm{exp}\left(\frac{E^{}E}{kT_r}\right)\frac{\overline{\rho }_{r1}(v_1,\mathrm{},v_n)}{\overline{\rho }_{r1}(v_1^{},\mathrm{},v_n^{})}].$$
(11)
This equation biases the a-priori probability of each dihedral angle with an estimate of its pd from a higher temperature. In previous literature Bruce (1985b); Binder et al. (1986b) such a biased updating has been used for the $`\varphi ^4`$ theory, where it is efficient to propose $`\varphi (i)`$ at each lattice size $`i`$ with its single-site probability.

For our temperatures $`T_r`$ the ordering (9) is assumed. With the definition $`\overline{\rho }_0(v_1,\mathrm{},v_n)=(2\pi )^n`$ the simulation at the highest temperature, $`T_1`$, is performed with the usual Metropolis algorithm. We have thus a recursive scheme, called rugged Metropolis (RM) in the following. When $`\overline{\rho }_{r1}(v_1,\mathrm{},v_n)`$ is always a useful approximation of $`\rho _r(v_1,\mathrm{},v_n)`$, the scheme zooms in on the native structure, because the pd at $`T_f`$ governs its fluctuations.

To get things working, we need to construct an estimator $`\overline{\rho }(v_1,\mathrm{},v_n;T_r)`$ from the numerical data of the RM simulation at temperature $`T_r`$. Although this is neither simple nor straightforward, a variety of approaches offer themselves to define and refine the desired estimators. In the following we work with the approximation

$$\overline{\rho }(v_1,\mathrm{},v_n;T_r)=\underset{i=1}{\overset{n}{}}\overline{\rho }_i^1(v_i;T_r)$$
(12)
where the $`\overline{\rho }_i^1(v_i;T_r)`$ are estimators of reduced one-variable pds defined by

$$\rho _i^1(v_i;T)=_\pi ^{+\pi }\underset{ji}{}dv_j\rho (v_1,\mathrm{},v_n;T).$$
(13)
The implementation of the resulting algorithm, called RM1, is straightforward, as estimators of the one-variable reduced pds are easily obtained from the time series of a simulation. The CPU time consumption of RM1 is practically identical with the one of the conventional Metropolis algorithm.

## 3 Numerical Results

To illustrate the developed ideas, we rely on the brain peptide Met-Enkephalin, which is a numerically well-studied Li and Scheraga (1987b); Okamoto et al. (1992b); Hansmann and Okamoto (1993b); Meirovitch et al. (1994b); Hansmann et al. (1999b). Met-Enkephalin is determined by the amino-acid sequence Tyr–Gly–Gly–Phe–Met or, in the short notation, Y–G–G–F–M, where

$`\mathrm{Tyr}(\mathrm{Y})`$ $``$ $`\mathrm{Tyrosine}`$
$`\mathrm{Gly}(\mathrm{G})`$ $``$ $`\mathrm{Glycine}`$
$`\mathrm{Phe}(\mathrm{F})`$ $``$ $`\mathrm{Phenylalanine}`$
$`\mathrm{Met}(\mathrm{M})`$ $``$ $`\mathrm{Methionine}`$
Our simulations are performed with a variant of SMMP SMMP (2001b) (Simple Molecular Dynamics for Protein) using fully variable $`\omega `$ torsion angles. For the data analysis we keep a times series by writing out configurations every 32 sweeps.

Fig.3 show the probability densities of the dihedral angle $`v_{13}`$ (for the notation see Table 1) at 400 K and 300 K. This angle is chosen because it illustrates the possibility of large moves in the Metropolis updating, which jump barriers of a molecular dynamics simulation. Namely, a single update may take us directly from each of the three populated regions to each other, whereas one encounters barriers when one has to move in small increments $`\mathrm{}v_{13}`$.

To evaluate the relative performance of different algorithms, we measure the integrated autocorrelation times $`\tau _{\mathrm{int}}`$ for the energy and each dihedral angle. The integrated autocorrelation times are directly proportional to the computer run times needed to achieve the same statistical accuracy. For an observable $`f`$ the autocorrelations are

$$C(t)=f_0f_tf^2$$
(14)
where $`t`$ labels the computer time. Defining $`c(t)=C(t)/C(0)`$, the time-dependent integrated autocorrelation time is given by

$$\tau _{\mathrm{int}}(t)=1+2\underset{t^{}=1}{\overset{t}{}}c(t^{}).$$
(15)
Formally the integrated autocorrelation time $`\tau _{\mathrm{int}}`$ is defined by $`\tau _{\mathrm{int}}=lim_t\mathrm{}\tau _{\mathrm{int}}(t)`$. Numerically this limit cannot be reached as the noise of the estimator increases faster than the signal. Nevertheless, one can calculate reliable estimates by reaching a window of $`t`$ values for which $`\tau _{\mathrm{int}}(t)`$ becomes flat, while its error bars are still reasonably small.

Columns 5 and 6 of Table 1 collect our estimates of the integrated autocorrelation times from conventional, canonical Metropolis simulations at 400 K and 300 K. We find a remarkable slowing down. For the energy, which is characteristic for the over-all behavior, $`\tau _{\mathrm{int}}`$ increases by a factor of about ten. For certain angles, e.g. $`v_{10}`$, $`\tau _{\mathrm{int}}`$ increases even by factors larger than twenty. At 400 K conventional Metropolis simulations allow to calculate observables with relative ease, while this is no longer the case at 300 K. Our aim is to use information from the 400 K simulation to improve on the performance at 300 K.

One finds that the energy histograms of the 400 K and the 300 K simulation have a considerable overlap, see Fig.4. Therefore, the two temperatures are well-suited to be combined into a PT Geyer PT (1996b); Hukusima and Nemeto PT (1996b) simulation. We abstain here from introducing additional temperatures, because our aim is to get a clear understanding of the improvement of the 300 K simulation due to PT input from 400 K.

In Fig.5 integrated autocorrelations times for the energy variable are compared at 300 K. The order of the curves agrees with the order in the figure legend. From up to down: Conventional canonical Metropolis simulation, RM1 improved canonical Metropolis simulation with input from 400 K, PT simulation coupled to 400 K and the RM1 improved PT simulation. The variable $`t`$ of equation (15) is given in units of 32 sweeps due to the way our data are recorded. In the range shown, i.e. up to $`t=200\times 32`$ sweeps, a window exist for the PT and the RM1 improved PT simulations, which allows to estimate the integrated autocorrelation times of these simulations. For the other two simulations one has to go to even larger $`t`$ values, but it remain possible to estimate $`\tau _{\mathrm{int}}`$. For the energy and all 24 dihedral angles the $`\tau _{\mathrm{int}}`$ estimates are collected in Table 1.

For the energy we find a decrease from $`\tau _{\mathrm{int}}50`$ to $`\tau _{\mathrm{int}}25`$ due to the RM1 improvement of the canonical Metropolis simulation, i.e. approximately a factor of 2. The PT improvement of this simulation is even larger, namely from $`\tau _{\mathrm{int}}50`$ to $`\tau _{\mathrm{int}}20`$, i.e. approximately a factor of 2.5. In both cases the CPU time spent at 400 K is not part of the equation. This is justified as one should anyway understand the system first at temperatures where one does not suffer from long autocorrelation times. For PT the factor 2.5 is the improvement in real time when two identical PC nodes with some parallel software like MPI (Message Passing Interface) are available. Most remarkably, the two improvements multiply: For the RM1 improved PT simulation the autocorrelation time is down to $`\tau _{\mathrm{int}}10`$, only about one more factor of two away from the $`\tau _{\mathrm{int}}`$ value at 400 K.

Inspection of the $`\tau _{\mathrm{int}}`$ values for the 24 dihedral angles shows large differences. For the conventional canonical simulation the range is from111 A value $`\tau _{\mathrm{int}}=1`$ means that our resolution of measurements every 32 sweeps is too crude to show the autocorrelations, which one may still expect on the scale of a few sweeps. $`\tau _{\mathrm{int}}1`$ for $`v_{21}`$ to $`\tau _{\mathrm{int}}200`$ for $`v_{16}`$. For the PT-RM1 simulation it becomes reduced to a range from $`\tau _{\mathrm{int}}1`$ to $`\tau _{\mathrm{int}}20`$. This suggests that it is not efficient to simulate by sweeps, where each dihedral angle is updated once per sweep. Instead, an algorithm (systematic or random) where the number of updates per angle is proportional to its integrated autocorrelation time is expected to be more efficient. Obviously, this can be implemented without major changes of the existing code, but tests of this idea have not yet been completed. After all these improvements the 300 K simulation is expected to be no more autocorrelated than the conventional, canonical simulation at 400 K.

All-atom Metropolis simulations of small biomolecules deserve a place on their own in the arsenal of computational biophysics. A good understanding of smaller pieces of a larger protein is one of the ingredients, which can pave the way towards the understanding of large systems like proteins. Major algorithmic problems remain to be solved before Metropolis simulations of all-atom models may become directly suitable for applications like the folding of small to medium sized proteins.

One algorithmic problem is that of correlated moves of two or more dihedral angles. Such moves promise to overcome (jump) essential free energy barriers in the case of larger molecules. Within conventional canonical Metropolis simulations the acceptance rates for simultaneous moves of two and more angles are prohibitively small. The RM concept promises major inroads, but details have not yet been tested out.

An even more challenging problem is the inclusion of solvent effects. Here the large moves of Metropolis updates create the cavity problem. We need a cavity in the surrounding water to accommodate the move and updates of the dihedral angles do not create one. This appears to be one of the major reasons why large scale simulations of proteins with a surrounding solvent are to an overwhelming extent done within the MD framework, where the cavity problem is avoided due to joint, small moves of all degrees of freedoms. Within the Metropolis approach effective solvent models may provide a solution. They are, e.g., defined by a solvent contribution like $`E_{\mathrm{sol}}`$ in Eq.(6). Unfortunately, no reliable determination of the parameters entering this equation exists presently in the literature, as was demonstrated by recent simulations Hansmann solvent (2003b); Berg and Hsu (2003b). However, it appears to be within the reach of Metropolis simulations to lead the way towards determining reliable parameterizations.

## 4 Summary and Conclusions

* The RM1 approximation to the rugged Metropolis (RM) method Berg (2003b) leads already to considerable improvements over conventional Metropolis simulations of Met-Enkephalin at 300 K. As RM1 is easily implemented and needs no additional computer time, it should be used whenever Metropolis simulations of suitable systems are done.
* For larger systems one-variable moves alone will not work due to correlations between the dihedral angles. The RM approach promises sufficiently large acceptance rates for multi-variable moves.
* Ultimately, each biomolecule of interest may need its own, specifically designed, Metropolis algorithm. This task includes to determine reliable parameters for a solvent model like the one of Eq.(6).

I am indebted to Yuko Okamoto for many useful discussion and to Robert Swendsen for kindly informing me during this meeting about a dynamically optimized Monte Carlo method Swendsen et al. (1992b), which is tailored for the simulation of biomolecules. Further, I would like to thank James Gubernatis and the other organizers of the Los Alamos Metropolis workshop for their kind hospitality. This work was partially supported by the U.S. Department of Energy under contract No. DE-FG02-97ER41022.


text:
# Surface photometry of LSB edge-on galaxies

## 1 Introduction

One of the main features of the low surface brightness galaxies (LSB) which make them different from the ”usual” high surface brightness (HSB) galaxies is that they are considered to be dark matter dominated objects.

We conducted the study of a sample of several edge-on LSB and HSB galaxies in order to compare their photometric parameters (including the stellar disk thickness) and the relative mass of their dark halos.

## 2 Sample of galaxies and observations

Our sample consists of several objects taken from the Revised Catalog of Flat Galaxies (, RFGC hereafter). All those galaxies included to the catalog are highly inclined objects. We selected four objects of faintest surface brightness class as candidates to the LSB galaxies (according to ), and seven objects with intermediate and high surface brightness class as ”reference” HSB objects.

The surface photometry observations were made with the 6-m telescope at the Special Astrophysical Observatory (Russia). The large aperture of the telescope provided us with high angular resolution in images (0”.2/pix) together with a good sensitivity to the very faint regions of the galaxies. The chosen galaxies and the set of calibrating frames were taked during one observing run in December 2000 using R Johnson-Cousins photometric band.

The combined and calibrated R images were used to obtain the radial scale length of a galactic disk ($`h`$), the disk vertical scale heights ($`z_0`$, assuming $`sech^2`$ law), the ”face-on” central surface brightness ($`\mu _0`$) of the disk (see  for the details), and the bulge to disk luminosity ratio $`Lb/Ld`$. The final value of $`\mu _0`$ was corrected for the extinction in our Galaxy (according to the LEDA). The fitting profiles have been convolved with the atmospheric smearing function.

The distribution of $`\mu _0`$ values indicates the presence of two subsamples: those with $`\mu _0`$ values greater than 23.5 $`mag/arcsec^2`$ which we defined as LSB and those with higher surface brightness designated HSB in this work. Note that all galaxies of the faintest surface brightness class (according to the RFGC) were included to our LSB objects.

Although our sample enables us to compare of the photometric parameters of LSB and HSB disks, the sample is very limited. We also incorporated one more sample of edge-on galaxies whose photometric parameters have been published by Barteldrees etc. . The authors used the similar passband and the same fitting functions to extract the photometric parameters.

## 3 LSB versus HSB: the vertical scale height of galactic disk <br>as a new feature to compare.

As it was shown in  and , the galaxies of lower surface brightness tend to show the higher $`z_0/h`$ ratio. However, mostly HSB galaxies were considered in those previous studies. Here we show the dependence of $`z_0/h`$ on the central surface brightness $`\mu _0`$ for our sample. Our galaxies are denoted by the squares in Fig.1. The open squares are for the HSB subsample whereas the filled ones designate the LSB galaxies. The galaxies taken from  are plotted as crosses. Futhermore, the near-infrared (K band) sample of edge-on galaxies from  is available for the comparison (the 2MASS sample hereafter). The systematic difference between the R and K images in $`z_0/h`$ (1.4 times, see  for details), in the surface brightness (2.4 $`mag`$, see ), and the internal extinction (1.2 $`mag`$) were taken into account.

Fig.1 shows all three samples together. The 2MASS sample is denoted by the small filled triangles. A trend in Fig.1 is seen well, in average of 2 mag difference in $`\mu _0`$ leads to 1.5 change in the ratio of the scales. Note that there is no correlation seen when $`h`$ or $`z_0`$ were drawn against $`\mu _0`$ separately.

Following , we calculated the ratio of total mass $`M_t`$ to B-band disk luminosity $`L_B`$: $`M_t=0.5G4hV_{m}^{}{}_{}{}^{2}/L_B`$. Here $`V_m`$ is the maximum of rotation curve (taken from the LEDA). The value of $`L_B`$ is inferred from the absolute B-magnitude of a galaxy taken from the LEDA and corrected for the internal absorption. In Fig. 3a one can see the values of $`M_t/L_B`$ plotted versus the ratio $`z_0/h`$. The notation in the figure is the same as in Fig.1. The three curves in Fig.3a are the same as in Fig.2 which were recalculated for the values of the mass to luminosity ratio $`M/L`$ of 1, 2 and 4. As it is seen in Fig.3a, most of the galaxies have a reasonable value of $`M/L`$ between 1 and 4. The mass to luminosity ratio is higher systematically for our LSB galaxies in comparison with HSB ones.

Since the luminosity ratio $`Lb/Ld`$ reflects the mass ratio, the mass of the dark halo $`M_h`$ can be estimated from the dependence in Fig. 3. We drew the ratio of dark to luminous matter $`M_h/(M_b+M_d)`$ for our galaxies ishown n Fig.3b. Here we designate the bulge mass as $`M_b`$. We show the ratio of ”dark” to ”luminous” mass $`M_h/(M_b+M_d)`$ for our galaxies. Here we designate the bulde mass as $`M_b`$. We assume that the bulge and the disk have roughly the same $`M/L`$. Surprisingly, there is no systematic difference between the values of $`M_h/(M_b+M_d)`$ for the galaxies of different central surface brightness.

We also enable to compare the mass of the spherical subsystem $`M_s`$ for our galaxies. In Fig. 3c we present how the ratio $`M_s/M_d`$ depends on the disk central surface brightness. We kept the same notation as in Fig. 1 and Fig. 3a. The figure indicates that the LSB galaxies do not have the most massive dark matter halos, but have the most massive spherical subsystem. That supports the result by Graham () that not all LSB galaxies are the dark-matter dominated objects.

## 4 Conclusions

1) We present the results of the photometric observations made for the sample of edge-on galaxies. Our sample contains four LSB galaxies as well as seven HSB ones. The photometric disk scales (vertical and radial), the disk central surface brightness and the bulge to disk luminosity ratio were determined.

2) Stellar disks of our LSB galaxies are thinner than HSB ones. There is a good correlation between their central surface brightness and their vertical to radial scales ratio.

3) Our LSB galaxies tend to harbor the massive spherical subsystems as well as to havehigher values of the mass-to-luminosity ratio in their disks when compared to the HSB objects. Nevertheless, the dark halo is not strictly the most massive subsystem in our LSB galaxies. The LSB galaxies appear to be the spherical dominated systems, but not the ”dark matter dominated” ones.

Acknowledgments D.B. is supported by NASA/JPL through grant NRA-99-04-OSS-058. The project was partially supported by Russian Foundation for Basic Research via grant 01-02-17597. We have made use of the LEDA database.
REFERENCES
 Karachentsev, I., Karachentseva, V., Kudrya, Y., et al.// 1999, Bull. Of Special Astrophys. Obs., 1999, 47, 5, (RFGC)  Bizyaev, D.// 2000, Astron. Lett. 26, 219.  van der Kruit, P., Searle, L. // 1981, A&A, 95, 105 and 116  Barteldrees, A., Dettmar, R.-J. // 1994, A&AS, 103, 475  Bizyaev, D., Mitronova, S. // 2002, A&A, 389, 795  Zasov, A., Bizyaev, D., Makarov, D. etc. // 2002, Astron.Lett. 28, 527  de Jong, R. // 1996, A&A, 313, 45  Zasov, A., Makarov, D., Mikhailova, E. // 1991, Astron.Lett., 17, 374  Mikhailova, E., Khoperskov, A., Sharpak, S. // 2001, in Stellar Dynamics: From Classic to Modern , Ed. by L.Ossipkov and I.Nikiforov, p. 147  Graham, A. // 2002, MNRAS, 334, 721

## FIGURE CAPTIONS

Fig.1. The vertical to radial scale ratio $`z_0/h`$ is shown against the disk central surface density $`\mu _0`$. The open squares are for the HSB subsample, the filled ones designate our LSB galaxies. The galaxies taken from  are shown by crosses, the 2MASS sample is designated by the small filled triangles.

Fig.2. The ratio od scales $`z_0/h`$ is plotted versus the relative mass of the spherical galactic subsystem $`M_s/M_d`$ according to the numerical simulations described in .

Fig.3. a) The scales ratio $`z_0/h`$ and the total mass to disk luminosity rati $`M_t/L_B`$. The same notation applies here as in Fig.1. The three curves show the model results (see Fig.2) recalculated for $`M/L`$ values of 1 (solid curve), 2 (dash-dotted) and 4 (dashed). The value of M/L is systematically higher for our LSB galaxies.

b) The ratio of dark to luminous mass $`M_h/(M_b+Md)`$ is plotted against the disk central surface brightness $`\mu _0`$ with the same notation as in Fig. 3a. There is no systematic difference in the relative dark matter mass between the LSB and HSB galaxies.

c) The ratio of the spherical to disk mass $`M_s/M_d`$ is shown in dependence on the disk central surface brightness. The LSB galaxies have the more massive spherical subsystem when compared to HSB ones.