text:
# Neural networks and photometric redshifts

## 1 Introduction

Ongoing and planned digital surveys such as the Sloan Digital Sky Survey and the implementation of an International Virtual Observatory, will enormously increase the quality and the amount of data available to the astronomical community. The scientific exploitation of these data, while promising new answers to old questions, is also stimulating the implementation of new tools capable to deal in an effective way with unprecedently large volumes of data. One field which will be deeply affected by these new data sets is that of the large scale of the universe. Planned or ongoing surveys such as the Sloan Digital Sky Survey (SDSS,;York et al. 2000), the VIRMOS-VLT Survey (Le Fèvre et al. 2000), the VST Survey (G. Busarello, private communication) will provide a huge amount of high accuracy spectroscopic and photometric data which will enable observational cosmologists to map with unprecedented accuracy and detail the properties and the structure of the Universe.

In this respect it needs to be stressed that, in spite of the recent and well known advances in multiobject spectroscopy, photometric redshifts (cf. Baum 1962, Pushell et al. 1982) derived from multicolor photometry have been for a long time and still are the only tool which may be effectively used to evaluate the distances of large number of galaxies. Most photometric redshifts methods rely on a $`\chi ^2`$ fitting of a library of template Spectral Energy Distributions (hereafter SEDs) to the observed data points, and differ mainly in how the SEDs are derived and on how they are fitted to the data.

SEDs may either be derived from population synthesis models (cf. Bruzual and Charlot 1993) or be spectra of real objects selected in order to ensure a sufficient coverage of morphological types and/or luminosity classes. As effectively stressed by Koo (1999), both approaches (synthetic and empirical) have their pro’s and con’s.

Synthetic spectra, for instance, sample an ’a priori’ defined grid of mixtures of stellar populations and may either include unrealistic combinations of parameters or exclude some unknown cases, while empirical templates are usually derived from nearby and bright galaxies and may therefore be not representative of the spectral properties of galaxies falling in other redshift ranges. The various methods are extensively compared and discussed in several papers (cf. Koo 1999; Fernandez-Soto et al. 2001) and have been applied to many different data sets such as the Hubble Deep Field (cf. Massarotti et al. 2001a, 2001b) .

Another approach, which is in the same line of the one discussed in this paper, can be applied only to what we shall call ’mixed surveys’, id est datasets where accurate and multiband photometric data for a large number of objects are supplemented by spectroscopic redshifts for a small but statistically significant subsample of the same objects. In this case, the spectroscopic data can be used to constrain the fit of a polynomial function mapping the photometric data (cf. Connolly et al. 1995, Wang et al. 1998, Brunner et al. 2000). It needs to be stressed that, at difference with the SED fitting methods, this interpolative approach cannot be effectively applied to objects fainter than the spectroscopic limit since, in absence of an ’a priori’ knowledge, impossible extrapolations would be required. It could be argued that the needed ’a priori’ knowledge could be extracted from population synthesis models, but it is apparent that, in this case, the uncertainties of the two methods would add up and SED’s fitting methods would - in any case - be more accurate and preferable.

Interpolative methods, however, offer the great advantage that they are trained on the real Universe and do not require strong assumptions on the physics of the formation and evolution of stellar populations. Neural Networks (hereafter NNs) are known to be excellent tools for interpolating data and for extracting patterns and trends (cf. the standard textbook by Bishop 1995) and in this paper, we shall discuss the application of a set of neural tools to the determination of photometric redshifts in large ”mixed surveys” (cf. Giordano 2001). In section 2 we introduce the basic concepts of Neural Networks paying special attention to the Multi Layer Perceptron and the Self Organising Maps; in Section 3 we discuss a first application to the SDSS Early Data Release data (Stoughton et al. 2001) and in Section 4 we show how SOM can be used to evaluate the degree of contamination of the final redshift catalogues. In Section 5, finally we shall draw our conclusions and discuss some possible applications of the neural tools. In a forthcoming paper we shall compare the results of different photometric redshifts methods applied to the same SDSS-EDR data and will discuss the properties of the photometric redshift catalogue derived with the method described in this paper (Longo et al. 2002, in preparation).

## 2 Neural Networks

NNs, over the years, have proven to be a very powerful tool capable to extract reliable information and patterns from large amounts of data even in the absence of models describing the data (cf. Bishop 1995) and are finding a wide range of applications also in the astronomical community: catalogue extraction (Andreon et al. 2001), star/galaxy classification (Bertin and Arnout, 1996, Andreon et al. 2001), galaxy morphology (Storrie-Lombardi et al. 1992; Lahav et al. 1996), classification of stellar spectra (Bailer and Jones 1998, Allende Prieto et al. 2000, Weaver 2000), data quality and data mining (Tagliaferri et al. 2002).

A NN is usually structured into an input layer of neurons, one or more ”hidden” layers and one output layer. Neurons belonging to adiacent layers are usually fully connected and the various types and architectures of the NNs are identified by the different topologies adopted for the connections and by the choice of the activation function (details can be found in the standard book by Bishop 1995). From the operational point of view, NNs can be divided into two main types, supervised and unsupervised systems, accordingly to the type of learning.

In the first case, NNs learn how to recognise a specific pattern or characteristic on a set (hereafter ”training set”) of labeled data containing for each input vector also the desired output (hereafter ”target”). In unsupervised systems, instead, NNs look for the statistical similarity of the imput data. In this work both types of NNs are used to perform different tasks.

The AstroMining software (Longo et al. 2001) is a package written in the MatLab © environment to perform a large number of data mining and knowledge discovery tasks, both supervised and unsupervised, in large multiparametric astronomical datasets. The package relies also on the Matlab © ”Neural Network”, the ”SOM” (Vesanto 1997) and the ”Netlab” (Nabney and Bishop 1998) toolboxes.

AstroMining accepts as input any ASCII table containing a header describing the contents of each colum and then a set of parameters. Via interactive interfaces, it is possible to perform a large number of operations: i) manipulation of the input data sets; ii) selection of relevant parameters; iii) selection of the type of neural architecture; iv) selection of the training validation and test sets construction procedure; v) etc. The package is completed by a large set of visualization and statistical tools which allow to estimate the reliability of the results and the performances of the network. The user friendliness and the generality of the package allow both a wide range of applications and the easy execution of experiments (more details on other aspects of the AstroMining tool which are not relevant to the present work may be found in Tagliaferri et al. 2002).

Let us focus now on some fundamental aspects connected with the use of supervised NNs. In order to perform correctly, almost all supervised NNs need to be trained, validated and tested on three independent datasets. In order to achieve good generalization performances, the training set needs to be representative of the typical data which will be passed to the network in the application phase. The validation data set (which is often and erroneously ignored in many NN applications) is a second dataset disjoined from the training set but having the same statistical properties. The role played by the validation set is subtle but crucial: by using the training set alone, in fact, the NNs (and the MLP in particular) may easily run into overfitting errors, thus loosing all generalization properties. It needs to be stressed in fact, that while the error computed on the training set may decrease asintotically, the capability of the network to reproduce patterns not encountered during the training phase may decrease (this is known as ’over-fitting’ condition), in other words, the NN will learn how to reproduce the patterns in the training set but will produce completely wrong results when applied to other data sets. The validation set prevents this from happening via a so called ’regularization technique’: during the training phase, at regular intervals, the training is interrupted (and the weights of the neurons are frozen), then the net is run on the validation set in order to compute the error with respect to the desired output; the training is stopped when the error computed on the validation set shows a significant increasing trend and, finally, the NN corresponding to the minimum error is selected.

After this phase the final performances of the resulting weight configuration are tested on a third data set, the so called ’test set’, which is, once more, completely disjoined from the previous two. Another regularization method not requiring a validation set, namely the Bayesian Learning Approach, will be described in the next paragraph.

### 2.1 The Multi Layer Perceptron - MLP

Due to its interpolation and capabilities, the Multi Layer Perceptron (MLP) is one of the most widely used neural architectures. As most other networks, the MLP is structured into an input layer, one or more hidden layers and one output layer (see figure 1). We implemented an MLP with one hidden layer and $`n`$ imput neurons, where $`n`$ is the number of parameters selected by the user as input in each experiment.

As just mentioned, it is also possible to train NN’s in a Bayesian framework, which allows to find the more efficient among a population of NN’s differing in the hyperparameters controlling the learning of the network (cf. Bishop 1995), in the number of hidden nodes, etc. The most important hyperparameters being the so called $`\alpha `$ and $`\beta `$. The parameter $`\alpha `$ is related to the weights of the network: a larger value for a component of $`\alpha `$ implying a less meaningful corresponding weight, thus allowing to estimate the relative importance of the different inputs (Automatic Relevance Determination; Bishop 1995) and, therefore, the selection of the input parameters which are more relevant to a given task. The parameter $`\beta `$ is instead related to the variance of the noise: a smaller value corresponding to a larger value of the noise and therefore to a lower reliability of the network. The Bayesian method allows the values of the regularization coefficients to be selected using only the training set, without the need for a validation set.

The implementation of a Bayesian framework requires several steps: initialization of weights and hyperparameters; training the network via a non linear optimization algorithm in order to minimize the total error function. Every few cycles of the algorithm, the hyperparameters are re-estimated and eventually the cycles are reiterated.

### 2.2 The Self Organizing Maps - SOM

The Self-Organizing Map (SOM), developed by Kohonen (1995), is one of the most used NN model. The SOM algorithm is based on unsupervised competitive learning, id est the training is entirely data-driven and all neurons of the map compete with each other producing only one winning neuron for each input vector. This property turns SOM into an ideal tool for KDD and expecially for its exploratory phase: data mining (Vesanto 1997). Among various other advantages, SOM allow an approximation of the probability density function of the training data, the derivation of prototype vectors best describing the data, and a highly visualized and user friendly approach to the investigation of the data.

A SOM is composed by neurons located on a regular, usually 1 or 2-dimensional grid. Each neuron $`i`$ of the SOM is represented by an $`n`$-dimensional weight or reference

$`m_i=[m_{i_1},m_{i_2},\mathrm{},m_{i_n}]^T`$

where $`n`$ is the dimension of the input vectors. Higher dimensional grids are not generally used since their visualization is much more problematic. Usually the map topology is a rectangle but also toroidal topologies have been used successfully. The neurons of the map are connected to adjacent neurons by a neighborhood relation dictating the structure of the map. In the 2-dimensional case the neurons of the map can be arranged either on a rectangular or on a hexagonal lattice. The number of neurons determines the granularity of the resulting mapping, which affects the accuracy and the generalization capability of the SOM.

The use of SOM for data mining requires several steps : construction, normalization and initialization of the Data Set, (unsupervised) training, visualization of the resulting map, and, finally, analysis of the results. The first two steps depend on the individual data set to be processed and the normalization is made in order to achieve $`mean=0`$ and $`variance=1`$.

In the basic SOM algorithm, the topological relations and the number of neurons are fixed from the beginning. The number of neurons should usually be selected by trial and error, with the neighborhood size controlling the smoothness and generalization of the mapping. Before training, in the course of the initialization phase, initial values are given to the weight vectors. The SOM is robust regarding the initialization, but if properly accomplished it allows the algorithm to converge faster to a good solution. Typically, any of the following initialization procedures may be used:

* random initialization, with the weight vectors initialized with small random values;
* sample initialization, where the weight vectors are initialized with random samples drawn from the input data set;
* linear initialization, where the weight vectors are initialized in an orderly fashion along the linear subspace spanned by the two principal eigenvectors of the input data set. These eigenvectors can be calculated using Gram-Schmidt procedure (Kohonen 1995).

During the training phase, one sample vector $`𝐱`$ from the input data set is randomly chosen and a similarity measure is calculated between it and all the weight vectors of the map. The Best-Matching Unit (BMU), denoted as $`c`$, is the unit with weight vector having the greatest similarity with the input sample $`𝐱`$. The similarity is usually defined by means of a distance measure, typically an Euclidean distance. Formally the BMU is defined as the neuron for which:

$`𝐱𝐦_c=\mathrm{min}_i𝐱𝐦_i`$

where $`.`$ denotes the distance measure. After finding the BMU, the weight vectors of the SOM are updated. The weight vectors of the BMU and its topological neighbors are moved in the direction of the input vector, in the input space. The SOM update rule for the weight vector of the unit $`i`$ is:

$`𝐦_i\left(t+1\right)=𝐦_i\left(t\right)+h_{ci}\left(t\right)\left[𝐱\left(t\right)𝐦_i\left(t\right)\right]`$

Where $`t`$ denotes the time, $`𝐱(t)`$ is the input vector and $`h_{ci}(t)`$ denotes the neighborhood kernel around the winner unit. The neighborhood kernel is a non-increasing function of time and of the distance of unit $`i`$ from the winner unit $`c`$. It defines the region of influence that the input sample has on the SOM. The kernel is composed by two parts: the neighborhood function $`h(d,t)`$ and the learning rate function $`\alpha (t)`$:

$`h_{ci}\left(t\right)=h(𝐫_c𝐫_i,t)\alpha (t)`$

where $`r_i`$ is the location of unit $`i`$ on the map grid. The neighborhood function used in our experiments is the Gaussian neighborhood function:

$`exp\left(\frac{𝐫_c𝐫_i^2}{2\sigma ^2\left(t\right)}\right)`$

The learning rate $`\alpha (t)`$ is a decreasing function of time. Two commonly used forms are a linear function and a function inversely proportional to time:

$`\alpha (t)=\frac{A}{t+B}`$

where $`A`$ and $`B`$ are some suitably selected constants. The training is usually performed into two phases. In the first phase, relatively large initial $`\alpha `$ value and neighborhood radius are used. In the second phase both the $`\alpha `$ value and the neighborhood are small from the beginning. This procedure corresponds to first tuning the SOM approximately to the same space as the input data and then fine-tuning the map. The SOM toolbox (Vesanto 1997) includes the tools for the visualization and analysis of SOM and, since the weight vectors are ordered on the grid, the visualization of the U matrix turns out to be expecially useful in the data understanding/survey phase. The U matrix visualizes the clustering structures of the SOM as distances (in the assumed metric) between neighboring map units, thus high values of the U-matrix indicate a cluster border, uniform areas of low values indicate clusters themselves.

In Figure 2, we show both the U matrix for the whole data set, and the structure of the individual components.

Another advantage of SOM is that it is relatively easy to label individual data, id est to identify which neuron is activated by a given input vector. The utility of these properties of the SOM will become clear in the next paragraphs.

## 3 Application to the SDSS-EDR data

A preliminary data release (Early Data Release or EDR) of the SDSS was made available to the public in 2001 (Stoughton et al. 2001). This data sets provides photometric, astrometric and morphological data for an estimated 16 millions of objects in two fields: an Equatorial $`2^{circ}`$ wide strip of constant declination centered around $`\delta `$=0 and a rectangular patch overlapping with the SIRTF First Look Survey.

The EDR provides also spectroscopic redshifts for a little more than $`50.000`$ galaxies distributed over a large redshift range and is therefore representative of the type of the data which will be produced by the next generation of large scale surveys. In order to build the training, validation and test sets, we first extracted from the SDSS-EDR a set of parameters ($`u`$, $`g`$, $`r`$, $`i`$, $`z`$, both total and petrosian magnitudes, petrosian radii, $`50\%`$ and $`90\%`$ petrosian flux levels, surface brightness and extinction coefficients; Stoughton et al. 2001) for all galaxies in the spectroscopic sample.

In this data set, redshifts are distributed in a very dishomogeneous way over the range $`07.0`$ (Figure 3; $`93\%`$ of the objects have $`z<0.7`$).

Table 1: Cuts applied to magnitudes and colors.    $`15<u<21.5`$ $`1.5<ug<5`$  $`15.0<g<21.7`$ $`1<gr<3`$  $`14.0<r<21.7`$ $`2<ri<2`$  $`14.0<i<20.8`$ $`2<iz<2`$  $`14.0<z<20.0`$

It needs to be stressed that the highly dishomogeneous distribution of the objects in the redshift space implies that the density of the training points dramatically decreases for increasing redshifts, and that: i) unless special care is paid to the construction of the training set, all networks will tend to perform much better in the range where the density of the training points is higher; ii) the application to the photometric data set will be strongly contaminated by the spurious determinations.

### 3.1 The construction of the training set

In order to achieve an optimal traning of the NNs, two different approaches to the construction of the training, validation and test sets were implemented: the uniform sampling and the clustered sampling (via K-means and/or SOM).

In both cases the training set data are first ordered by increasing redshift, then, in the case of uniform sampling, after fixing the number of training objects (which needs in any case to be smaller than 1/3 of the total sample) objects are extracted following a decimation procedure. This approach however, is undermined by the fact that the input parameter space is not necessarily uniformously sampled, thus causing a loss in the generalization capabilities of the network.

In the clustered sampling method, objects in each redshift bin are first passed to a SOM or a K-means algorithm which performs an unsupervised clustering in the parameter space looking for the most significant statistical similarities in the data. Then, in each bin and for each cluster, objects are extracted in order to have an uniform sampling of the parameter space. This second procedure, while being slower than the uniform sampling allows a complete and statistically homogeneous coverage of the parameter space.

### 3.2 The photometric redshift evaluation

In order to evaluate the performances of the software as close as possible to the detection limit of the data, we did not introduce strong cuts on the limiting magnitudes. Hence the filters listed in Table 1, were applied to the magnitudes and to the colors. The latter were introduced in order to remove a few spurious objects present in the original data set.

The experiments were performed using the NNs in the Matlab and Netlab Toolboxes, with and without the Bayesian framework. All NNs had only one hidden layer and the experiments were performed varying the number of the input parameters and of the hidden units. Extensive experiments lead us to conclude that the Bayesian framework provides better generalization capabilities with a lower risk of overfitting, and that an optimal compromise between speed and accuracy is achieved with a maximum of 22 hidden neurons and 10 Bayesian cycles.

In Table 2, we summarize some of the results obtained from the experiments and, in Figure 4, we compare the spectroscopic redshifts versus the photometric redshifts derived for the test set objects in the best experiment.

### 3.3 Contamination of the catalogues

In practical applications, one of the most important problems to solve is the evaluation of the contamination of the final photometric redshift catalogues or, in other words, the evaluation of the number of objects which are erroneously attributed a $`z_{phot}`$ significantly (accordingly to some arbitrarily defined treshold) different from the unknown $`z_{spec}`$. This problem is usually approached by means of extensive simulations. The problem of contamination is even more relevant in the case of NNs based methods, since NNs are necessarily trained only in a limited range of redshifts and, when applied to the real data, they will produce misleading results for most (if not all) objects which ”in the real word” have redshifts falling outside the training range. This behaviour of the NNs is once more due to the fact that while being good interpolation tools, they have very little, if any, extrapolation capabilities. Furthermore, in mixed surveys, the selection criteria for the spectroscopic sample tend to favour the brightest (and, on average, the closer) galaxies with respect to the fainter and more distant ones and, therefore, the amount of contamination encountered, for instance, in the test set sets only a lower limit to the percentage of spurious redshifts in the final catalogue.

To be more specific: in the SDSS-EDR spectroscopic sample, over a total of 54,008 objects having $`z>0`$, only $`88\%`$, $`91\%`$ and $`93\%`$ have redshift $`z`$ lower than, respectively than $`0.3`$, $`0.5`$ and $`0.7`$. To train the network on objects falling in the above ranges implies , respectively, a minimum fraction of $`12\%`$, $`9\%`$ and $`7\%`$ of objects in the photometric data set having wrong estimates of the photometric redshift. On the other hand, as we have shown, the higher is the cut in redshifts, the lower is the accuracy and a compromise between these two factors needs to be found on objective grounds.

An accurate estimate of the contamination may be obtained using unsupervised SOM clustering techniques over the training set.

In Figure 5 we show the position of the BMU as a function of the redshift bin. Each exagon represents a neuron and it is clearly visible that low redshift ($`z<0.5`$) tend to activate neurons in the lower right part of the map, intermediate redshift ones ($`0.5<z<0.7`$) neurons in the lower left part and, finally, objects with redshift higher than $`0.7`$ activate only the neurons in the upper left corner. The labeling of the neurons (shown in the upper left map) was done using the training and validation data sets in order to avoid overfitting, while the confidence regions were evaluated on the test set.

Therefore, test set data may be used to map the neurons in the equivalent of confidence regions and to evaluate the degree of contamination to be expected in any given redshift bin. Conversely, when the network is applied to real data, the same confidence regions may be used to evaluate whether a photometric redshift correspondent to a given input vector may be trusted upon or not.

The above derived topology of the network is also crucial since it allows to derive the amount of contamination. In order to understand how this may be achieved, let us take the NN whose topological properties are shown in Figure 5, and consider the case of objects which are attributed a redshifts $`z_{phot}`$ $`<0.5`$. This prediction has a high degree of reliability only if the input vector activates a node in the central or right portions of the map. Vector producing a redshift $`z_{phot}`$ $`<0.5`$ but activating a node falling in the upper left corner of the map are likely to be misclassified. In Figure 6, we plot the photometric versus spectroscopic redhift for all test set objects having $`z_{phot}`$ $`<0.5`$ and activating nodes in the correct region of the map.

As it can be seen, out of 9270 objects with $`z_{phot}`$ $`<0.5`$, only 39 (id est, 0.4% of the sample) have discordant spectroscopic redshift. A confusion matrix helps in better quantifying the quality of the results. In Table 3, we give the confusion (or, in this case, ’contamination’) matrix obtained dividing the data in three classes accordingly to their spectroscopic redshifts, namely class I: $`0<z<0.3`$, class II: $`0.3<z<0.5`$, class III: $`z>0.5`$. The elements on the diagonal are the correct classification rates, while the other elements give the fraction of objects belonging to a given class which have been erroneously classified into an other class.

As it can be seen, in the redshift range $`(0,0.3)`$, 95.4% of the objects are correctly identified and only $`4.6\%`$ is attributed a wrong redshift estimate. In total, $`94.2\%`$ are correctly classified. By taking into account only the redshift range $`0<z<0.5`$, this percentage becomes $`97.3\%`$. From the confusion matrix, we can therefore derive a completeness of $`97.8\%`$ and a contamination of about $`0.5\%`$.

A simple step filter applied to the accepted BMU’s allows therefore to optimise the filter performances. For instance, it allows to choose whether to minimize the number of misclassified objects (thus reducing the completeness) or to minimize the classification error in a given redshift bin more than in another one.

Another possible use of the topological properties of the SOM will be discussed in a forthcoming paper and concerns the use of BMU to choose for each given input vector the optimal NN.

## 4 Summary and conclusions

The application of NNs to mixed data, id est spectroscopic and photometric surveys, allows to derive photometric redshifts over a wide range of redshifts with an accuracy equal if not better to that of more traditional techniques.

The method makes use of three different neural tools: i) an unsupervised SOM used to cluster the data in the training, validation and test set in order to ensure a complete coverage of the input parameter space; ii) a MLP in Bayesian framework used to estimate the photometric redshifts; iii) a supervised SOM used to derive the completeness and the contamination of the final catalogues. On the SDSS-EDR, the best result (interq. error = 0.020) was obtained by a MLP with 1 hidden layer of 22 neurons, after 5 Bayesian cycles. Once they are trained, NNs are extremely effective in terms of computational costs (the 16 million objects in the SDSS-EDR are processed in less than 50 min on a laptop; Longo et al. 2002, in preparation).

The method fully exploits the wealth of data provided by the new digital surveys since it allows to take into account not only the fluxes, but also the morphological and photometric parameters.

The proposed method will be particularly effective in mixed surveys, id est, in surveys were a large amount of multiband photometric data is complemented by a small subset of objects having also spectroscopic redshifts. It needs also to be stressed that the foreseen implementation of the Virtual Observatory will provide an ideal framework to NN based data mining tools: the availability of continuously updated data sets from where to extract reliable and extensive training sets will allow a widespread use of supervised NNs for the fast and accurate derivation of secondary parameters such as the photometric redshifts.


format: text

text:
# Dielectric behaviors of corrugated membranes

## I Introduction

It is known that lipid bilayers (abound in living cells membranes) exhibit a ripple phase of the bilayer-water interface, in a narrow temperature range . The ripple phase is characterized by permanent wave-like deformations of the interfaces, the origin of which is still a widely debated topic. However, it was observed that a giant dielectric dispersion can occur in the radiofrequency range, signified by a large dielectric increment in that frequency range . In practice, these membrane systems are under intensive research because they are responsible for delivery and retention of drugs .

We aim to obtain some results on the dielectric dispersion spectrum of corrugated membranes, by using our recently developed Green’s function formalism of periodic interfaces . In that formalism, we obtained an analytic expression for the Greenian that includes the effect of periodicity. In this work, we extend the Green’s function formalism to compute the local field distribution for a lipid bilayer membrane of arbitrary shape, separating two media of different dielectric constants. We will calculate the effective dielectric constant of the membrane subject to an ac applied field.

## II Green’s Function Formalism

The Green’s function formalism has been published recently . Here we re-iterate the formalism to establish notation. We will apply the formalism to a single interface and then extend to a bilayer membrane. The electrostatic potential satisfies the Laplace’s equation:

$`[ϵ(𝐫)\varphi (𝐫)]=4\pi \rho (𝐫),`$ (1)
with standard boundary conditions on the interface, where $`\rho (𝐫)`$ is the free charge density, $`ϵ(𝐫)`$ equals $`ϵ_2`$ in the host and $`ϵ_1`$ in the embedding medium.

Let $`V_1`$ and $`V_2`$ be the volume of the embedding and host medium, separated by an interface $`S`$. Denoting $`\theta (𝐫)=1`$ if $`𝐫V_1`$ and 0 otherwise, leads to an integral equation :

$`[1u\theta (𝐫)]\varphi (𝐫)=\varphi _0(𝐫)+{\displaystyle \frac{u}{4\pi }}{\displaystyle _S}𝑑S^{}\left[\widehat{𝐧}^{}^{}G(𝐫,𝐫^{})\right]\varphi (𝐫^{}),`$ (2)
where $`u=1ϵ_1/ϵ_2`$, $`\widehat{𝐧}^{}`$ is unit normal to $`S`$, $`G(𝐫,𝐫^{})=1/|𝐫𝐫^{}|`$ and $`\varphi _0`$ is the solution of $`^2\varphi _0=4\pi \rho /ϵ_2`$.

Accordingly, our approach aims to solve a surface integral equation for the potential at the expense of a two-step solution :

1. step 1: determine $`\varphi (𝐫)`$ for all r $`S`$ by solving Eq.(2), and then
2. step 2: obtain $`\varphi (𝐫)`$ for all r by using Eq.(2) and the results of step 1.

In step 1, we encounter a singularity when the integration variable $`𝐫^{}`$ approaches the point of observation r. To circumvent the problem, we take an infinitesimal volume around r and perform the surface integral analytically, we find

$`\left(1{\displaystyle \frac{u}{2}}\right)\varphi (𝐫)=\varphi _0(𝐫)+{\displaystyle \frac{u}{4\pi }}{\displaystyle _S^{}}𝑑S^{}\left[\widehat{𝐧}^{}^{}G(𝐫,𝐫^{})\right]\varphi (𝐫^{}),𝐫S,`$ (3)
where “prime” denotes a restricted integration which excludes $`𝐫^{}=𝐫`$. The (surface) integral equation (3) can be solved for $`\varphi (𝐫S)`$.

## III Application to a periodic interface

Here we apply the integral equation formalism to a periodic interface. Suppose the interface profile depends only on $`x`$, described by $`y=f(x)`$, where $`f(x)`$ is a periodic function of $`x`$ with period $`L`$: $`f(x+L)=f(x)`$. Without loss of generality, we will let $`L=1`$ in subsequent studies. Thus medium 1 occupies the space $`y<f(x)`$ while medium 2 occupies the space $`y>f(x)`$ separated by the interface at $`y=f(x)`$. The external field is $`𝐄_0`$ and $`\varphi _0(𝐫)=𝐄_0𝐫`$ is the potential. For a periodic system, $`\varphi (𝐫)`$ is a periodic function of the lattice vector $`𝐓`$. In what follows, we adopt similar treatment as the Korringa, Kohn and Rostoker (KKR) method  and rewrite the integral equation as:

$`\left(1{\displaystyle \frac{u}{2}}\right)\varphi (𝐫)=𝐄_0𝐫+{\displaystyle \frac{u}{4\pi }}{\displaystyle _S^{}}𝑑S^{}\stackrel{~}{G}(𝐫,𝐫^{})\varphi (𝐫^{}),`$ (4)
where the integration is performed within a unit cell. The structure Green’s function (Greenian) is given by :

$`\stackrel{~}{G}(𝐫,𝐫^{})={\displaystyle \underset{𝐓}{}}\widehat{𝐧}^{}^{}G(𝐫,𝐫^{}+𝐓)={\displaystyle \underset{m}{}}{\displaystyle \frac{2\widehat{𝐧}^{}(𝐫𝐫^{}m\widehat{𝐱})}{|𝐫𝐫^{}m\widehat{𝐱}|^2}}.`$ (5)
We were able to evaluate the Greenian analytically :

$`\stackrel{~}{G}(x,y;x^{},f(x^{}))={\displaystyle \frac{2\pi [f^{}(x^{})\mathrm{sin}2\pi (xx^{})\mathrm{sinh}2\pi (yf(x^{}))]}{\mathrm{cos}2\pi (xx^{})\mathrm{cosh}2\pi (yf(x^{}))}}.`$ (6)
Eq.(6) is a truly remarkable result – the analytic expression is valid for an arbitrary interface profile. If the point of observation $`(x,y)`$ is located at the interface, the Greenian has a finite limit as $`x^{}x`$:

$`\stackrel{~}{G}(x,f(x);x^{},f(x^{})|x^{}x)={\displaystyle \frac{f^{\prime \prime }(x)}{1+(f^{}(x))^2}}.`$ (7)
We first solve Eq.(4) for the potential $`\varphi (x,f(x))`$ right at the interface:

$`\left(1{\displaystyle \frac{u}{2}}\right)\varphi (x,f(x))=E_0f(x)+{\displaystyle \frac{u}{4\pi }}{\displaystyle 𝑑x^{}\stackrel{~}{G}(x,f(x);x^{},f(x^{}))\varphi (x^{},f(x^{}))}.`$ (8)
Then we use Eq.(2) to find the potential at any arbitrary point $`(x,y)S`$, using the potential at the interface.

$`(1u)\varphi (x,y)=E_0y+{\displaystyle \frac{u}{4\pi }}{\displaystyle 𝑑x^{}\stackrel{~}{G}(x,y;x^{},f(x^{}))\varphi (x^{},f(x^{}))},`$ (9)
$`\varphi (x,y)=E_0y+{\displaystyle \frac{u}{4\pi }}{\displaystyle 𝑑x^{}\stackrel{~}{G}(x,y;x^{},f(x^{}))\varphi (x^{},f(x^{}))},`$ (10)
for $`(x,y)V_1`$ and $`(x,y)V_1`$ respectively.

## IV Application to corrugated membranes

Here we extend the formalism to a bilayer membrane. Consider two interface profiles described by $`y=f_t(x)`$, where $`t=0,1`$ denote the lower and upper interface profiles respectively. Again $`f_t(x)`$ is a periodic function of $`x`$. Thus medium 1 occupies the space $`f_0(x)<y<f_1(x)`$ while medium 2 occupies the space $`y<f_0(x)`$ and $`y>f_1(x)`$. For the upper (lower) profile, $`\widehat{𝐧}\widehat{𝐲}>0`$ ($`\widehat{𝐧}\widehat{𝐲}<0`$), thus the Greenian becomes

$`\stackrel{~}{G}_{tt^{}}(x,f_t(x);x^{},f_t^{}(x^{}))={\displaystyle \frac{2\pi (1)^t^{}[f_t^{}^{}(x^{})\mathrm{sin}2\pi (xx^{})\mathrm{sinh}2\pi (f_t(x)f_t^{}(x^{}))]}{\mathrm{cos}2\pi (xx^{})\mathrm{cosh}2\pi (f_t(x)f_t^{}(x^{}))}}.`$ (11)
The effective dielectric constant $`ϵ_e`$ of the bilayer membrane satisfies the relation:

$$ϵ_eE_0V=_Vϵ(𝐫)E_y𝑑V=ϵ_2E_0Vuϵ_2_{V_1}E_y𝑑V,$$

where $`V_1`$ is the volume of the embedded medium. As $`𝐄=\varphi `$, the volume integration can be converted into a surface integration by the Green’s theorem. Moreover, for the upper (lower) profile, $`\widehat{𝐧}\widehat{𝐲}>0`$ ($`\widehat{𝐧}\widehat{𝐲}<0`$), thus the effective dielectric constant becomes

$`ϵ_e=ϵ_2{\displaystyle \frac{uϵ_2V_1}{E_0V}}{\displaystyle _0^1}𝑑x\left[\varphi (x,f_1(x))\varphi (x,f_0(x))\right].`$ (12)
To solve the integral equation, we express the potential at an arbitrary point into a mode expansion:

$`\varphi (x,y)={\displaystyle \underset{jk}{}}C_{jk}\psi _j(x)\xi _k(yf_t(x)),`$ (13)
where $`\psi _j(x)`$ and $`\xi _k(y)`$ are mode functions. The potential on the interfaces suffices:

$`\varphi (x,f_t(x))={\displaystyle \underset{jk}{}}C_{jk}\psi _j(x)\xi _k(0)={\displaystyle \underset{j}{}}A_j\psi _j(x),`$ (14)
where $`A_j=_kC_{jk}\xi _k(0)`$. Here we make a few remarks on the choice of the mode functions. The choice of the mode function is somewhat arbitrary in theory. In practice, these functions should be simple and easy to use. Common choice ranges from extended mode functions like the Fourier series expansions to localized mode functions like the step and triangular functions . Substituting the mode expansion Eq.(14) into Eq.(4), the coefficients $`A_i`$ satisfy the matrix equation:

$`\left[\left(1{\displaystyle \frac{u}{2}}\right)𝖡{\displaystyle \frac{u}{4\pi }}𝖬\right]𝐀=E_0𝐕,`$ (15)
where

$`B_{ij}`$ $`=`$ $`{\displaystyle 𝑑x\psi _i(x)\psi _j(x)},`$ (16)
$`M_{ij}`$ $`=`$ $`{\displaystyle 𝑑x𝑑x^{}\psi _i(x)\stackrel{~}{G}_{tt^{}}(x,f_t(x);x^{},f_t^{}(x^{}))\psi _j(x^{})},`$ (17)
$`V_i`$ $`=`$ $`{\displaystyle 𝑑x\psi _i(x)f_t(x)}.`$ (18)
It should be remarked that the mode functions need not be orthonormal and the matrix B is non-diagonal in general.

## V Numerical results

As a model bilayer membrane, we adopt the interface profiles:

$$f_t(x)=a\left(\mathrm{cos}2\pi x0.1\mathrm{sin}4\pi x\right)+(t0.5),t=0,1,$$

where $`a`$ is the amplitude of corrugation, and the sine function is added to upset the reflection symmetry about $`x=0`$ . The width of the bilayer membrane is thus unity. We adopt the step functions for the mode expansions:

$$\theta _i^h(x)=1,x_i\frac{h}{2}<x<x_i+\frac{h}{2}\mathrm{and}\theta _i^h(x)=0,\mathrm{otherwise},$$

where $`h`$ is the width of the step function. In what follows, we adopt 100 step functions both for the lower and upper profiles, equally spaced in the unit interval $`x(1/2,1/2]`$. The integrals Eqs.(16)–(18) can be readily performed.

To study the dielectric behavior, we apply an ac field at a frequency $`\omega `$. The embedded medium inside the membrane has a complex dielectric constant

$$\stackrel{~}{ϵ_1}=ϵ_1+\frac{\sigma _1}{i\omega },$$

where $`ϵ_1`$ and $`\sigma _1`$ are the dielectric constant and conductivity of the embedded medium respectively, with $`\omega `$ being the frequency of the applied field. We adopt the following parameters in the calculation: $`ϵ_1=40,\sigma _1=160`$, while $`ϵ_2=80`$. Also let $`V_1/V=1`$. The Maxwell-Wagner relaxation time of a planar interface is given by:

$$\tau =ϵ_1/\sigma _1.$$

In Fig.1, we plot (a) the real and (b) imaginary parts of the complex effective dielectric constant $`ϵ_e`$ normalized to $`ϵ_2`$ as function of frequency for various amplitude of corrugation $`a`$ ranging from 0.1 to 1.0. As is evident from Fig.1, there is a giant dielectric dispersion as the amplitude of corrugation becomes large ($`a1.0`$).

In summary, we have employed the Green’s function formalism to study the dielectric behavior of a corrugated membrane. The integral equation is solved and the dielectric dispersion spectrum is obtained for a periodic corrugated membrane. We should remark that the present formalism can readily be generalized to multi-layers systems as well as corrugations in two dimensions .

## Acknowledgments

This work was supported by the Research Grants Council of the Hong Kong SAR Government under grant CUHK 4245/01P. K. W. Yu acknowledges useful conversation with Professor Hong Sun.

Fig.1/Yu


text:
# Radiative capture of polarized neutrons by polarized protons

## Abstract

A model-independent irreducible tensor approach to $`\stackrel{}{p}(\stackrel{}{n},\gamma )d`$ is presented and an explicit form for the spin-structure of the matrix $`𝑴`$ for the reaction is obtained in terms of the Pauli spin-matrices $`𝝈(n)`$ and $`𝝈(p)`$. Expressing the multipole amplitudes in terms of the triplet $``$ triplet and singlet $``$ triplet transitions, we point out how the initial singlet and triplet contributions to the differential cross section can be determined empirically.

The study of radiative capture of neutrons by protons is of interest not only to nuclear physics as a testing ground for theories of $`NN`$ interaction, but also to astrophysics, where the fusion reaction is part of the proton-proton chain responsible for the generation of solar energy and production of elements in the early universe astro . The cross section for the process was directly measured cs  for the first time at neutron energy of 550 keV, although measurements cs1  with thermal neutrons have been carried out earlier. The 10% discrepancy noted between such measurements and theory was accounted for with surprising accuracy by the inclusion mec  of meson exchange currents (MEC). Measurements cs2  at higher energies were used to test effects of $`\rho `$ and $`\omega `$ exchanges and relativistic corrections to the impulse approximation ia . Although a thorough review of the inverse reaction $`d(\gamma ,n)p`$ over a wide range of energies is found in inv , the energy region just above threshold does not seem to have received much attention. Photodisintegration experiments at low energies focused attention inv-expts  on the relative $`M1`$ and $`E1`$ contributions. Experiments employing polarized photons have been reported between 5 to 10 MeV phot  and more recently phot-rec  at 3.58 MeV. The angular distribution as well as polarization of the neutron in $`d(\gamma ,p)n`$ were measured angular . Though the measured angular distribution of the neutron was found ang-neu  to be in good agreement with theory at 2.75 MeV, the angular distribution and neutron polarization in angular  differ from theory which includes the MEC contributions. Employing polarized neutrons at 6 MeV and 13.4 MeV, the analyzing powers ang-neu  in $`p(\stackrel{}{n},\gamma )d`$ was measured which differed from theory. Recent cross section calculations recent1 ; recent2 , which agree with each other within 5% deviation, are found to differ from the 1967 estimates of Fowler et al astro . The cross section was obtained in recent1  by fitting the existing data with a polynomial expansion, while the calculation in recent2  includes MEC’s, isobar currents and pair currents. The theory is in good agreement with cross section data for neutrons with energy above 14 MeV, but deviates by about 15% from older $`d(\gamma ,p)n`$ measurements between 2.5 MeV and 2.75 MeV old , which correspond respectively to neutron energies 550 keV and 1080 keV in $`p(n,\gamma )d`$. Apart from the fact that uncertainties in $`p(n,\gamma )d`$ cross sections at energies below about 600 keV lead to dominant uncertainties in the determination of the relative abundances of elements in the early universe, it is interesting to study the fusion reaction for its own sake, especially since a beginning has also been made cs  to study $`p(n,\gamma )d`$ precisely in the region of energies of interest to astrophysics. This recent study cs  has shown experimentally that the $`M1`$ transition from S-wave capture, which dominates at lower energies and the $`E1`$ transition from P-wave capture at higher energies are comparable at 550 keV.

It is interesting to observe that the neutron capture in $`p(n,\gamma )d`$ from the $`{}_{}{}^{3}S_{1}^{}`$ state is from the isospin $`I=0`$ state, which is the same for the deuteron, while that from the $`{}_{}{}^{1}S_{0}^{},I=1`$ state is characterized by a $`\mathrm{\Delta }I=1`$ transition. Likewise, capture from $`{}_{}{}^{1}P_{1}^{},I=0`$ state leads to a $`\mathrm{\Delta }I=0`$ transition, whereas those from $`{}_{}{}^{3}P_{j=0,1,2}^{},I=1`$ states are characterized by $`\mathrm{\Delta }I=1`$ transitions (see TABLE 1). The cross section at very low energies is dominated by the $`\mathrm{\Delta }I=1`$ amplitude for the $`M1`$ transition from the $`{}_{}{}^{1}S_{0}^{}`$ state, since the $`NN`$ scattering length in the $`{}_{}{}^{1}S_{0}^{}`$ channel is large and the $`M1`$ amplitude is proportional to the isovector magnetic moment of the nucleon, which is more than five times larger than the isoscalar magnetic moment. It was, however, noted quite early breit-rustgi  that the spin-dependent effects in the fusion reaction are sensitive to the small $`\mathrm{\Delta }I=0`$ amplitudes $`M1`$ and $`E2`$. Hence a “polarized-target-beam test” was proposed for the $`{}_{}{}^{3}S_{1}^{}^3S_1`$ radiative transitions in thermal $`np`$ capture. The isoscalar $`M1`$ and $`E2`$ amplitudes have also been studied theoretically recently J-W , using effective field theory, since the circular polarization bazhanov  of photons emitted in the capture of polarized neutrons by unpolarized protons and the angular distribution muller  in the capture of polarized neutrons by polarized protons are sensitive to the presence of these amplitudes. Parameter-free predictions, employing the Weinberg wein  scheme of power counting, have also been made park  for the spin-dependent observables in $`\stackrel{}{p}(\stackrel{}{n},\gamma )d`$ at threshold energies. However, the recently muller  measured value of $`\eta =(1.0\pm \mathrm{\hspace{0.17em}2.5})\times 10^4`$ for $`\gamma `$-anisotropy at $`50.5\%`$ polarization of neutrons and protons in $`\stackrel{}{p}(\stackrel{}{n},\gamma )d`$ cannot test the effective field theory predictions at $`10^7`$. Although this measurement provides the first experimental value for a spin-dependent observable in $`\stackrel{}{p}(\stackrel{}{n},\gamma )d`$, spin-observables in elastic $`\stackrel{}{N}\stackrel{}{N}`$ scattering have been measured for more than two decades bug . More recently, attention was focused on such observables in $`\stackrel{}{N}\stackrel{}{N}`$ scattering, with a veiw to determine the exact strength of the tensor interaction and the case of $`\pi ^{}`$ production in $`\stackrel{}{p}\stackrel{}{p}`$ collisions pp  for which the threshold itself is high. Even at higher energies, it is clear that the singlet and triplet radiative captures from any arbitrary initial partial wave, $`\mathrm{}`$ lead to $`p(n,\gamma )d`$ characterized respectively by $`\mathrm{\Delta }I=0`$ and $`\mathrm{\Delta }I=1`$ or $`\mathrm{\Delta }I=1`$ and $`\mathrm{\Delta }I=0`$ depending on whether the initial parity $`(1)^{\mathrm{}}`$ is odd or even. Since the $`NN`$ interaction, as elicited from elastic scattering data, conserves channel spin, $`s`$ (apart from total angular momentum $`j`$) as well as isospin $`I`$, it would be of interest to study experimentally the relative strengths of the initial singlet and triplet contributions to the $`p(n,\gamma )d`$ reaction at any given energy, since the fusion reaction leads to transitions in the two-nucleon system, which change spin $`s`$ as well as isospin $`I`$. It is also worth noting that $`np`$ fusion near threshold bears some resemblance to the case of $`pd`$ fusion, where the validity of the so-called “no-quartet theorem” cjr  was questioned and led to several incisive theoretical studies friar  in later years. More recently, it was pointed out grpndsp  that a non-zero tensor analyzing power in $`\stackrel{}{p}(\stackrel{}{d},\gamma )^3He`$ by itself provides a clear signature to the contributions from the quartet amplitudes. Moreover, it was shown grpnd  that it is possible to determine empirically the individual doublet and quartet differential cross sections at any given energy through appropriate measurements of the relevant spin-observables in $`\stackrel{}{N}\stackrel{}{d}`$ fusion. In contrast to the difficult $`\stackrel{}{N}\stackrel{}{d}`$ fusion experiments suggested in grpnd , it should be technically more feasible to carry out the $`\stackrel{}{p}(\stackrel{}{n},\gamma )^2H`$ experiments muller .

The purpose of this Rapid Communication is to present a model-independent theoretical approach to the $`np`$ fusion process based on the irreducible tensor formalism and identify observables in $`\stackrel{}{p}(\stackrel{}{n},d)\gamma `$ experiments that facilitate determination of the singlet and triplet cross sections empirically for the fusion reaction at any energy.

Let $`𝒑(p,0,0)`$ and $`𝒌(k,\theta ,0)`$ denote respectively the neutron and photon momenta in a right-handed c.m. frame mad  and let $`m_d;𝒌\mu |𝑻|𝒑;m_n,m_p`$ denote the elements of the on-energy-shell $`𝑻`$-matrix for the fusion reaction $`p(n,\gamma )d`$ where $`m_n,m_p,m_d`$ denote respectively the spin-projections of the neutron, proton and deuteron and $`\mu =\pm 1`$ denote the left- and right-circular states of photon polarization, as defined by Rose rose . The unpolarized differential cross section for the reaction is then given by

$$\begin{array}{cc}\hfill \frac{\text{d}\sigma _0}{\text{d}\mathrm{\Omega }}=& \left(\frac{k}{2\pi E}\right)^2\frac{E_nE_pE_d}{p}\frac{1}{4}\underset{m_n,m_p=\frac{1}{2}}{\overset{\frac{1}{2}}{}}\underset{m_d=1}{\overset{1}{}}\hfill \\ & \times \underset{\mu =1,1}{}|m_d;𝒌\mu |𝑻|𝒑;m_n,m_p|^2,\hfill \end{array}$$
(1)
where $`E_p,E_n`$ and $`E_d`$ denote the c.m. energies of the proton, neutron and $`{}_{}{}^{2}H`$ respectively, when the reaction takes place at c.m. energy $`E`$. Expressing the initial $`np`$ system in terms of isosinglet and isotriplet states and making use of the standard multipole expansion rose  for the photon in the final state and the usual partial wave expansion for relative motion in the initial state with channel-spin $`s=0,1`$, the $`𝑻`$-matrix may be expressed as

$$𝑻=\underset{s=0,1}{}\underset{\lambda =|1s|}{\overset{(1+s)}{}}\underset{\mu =1,1}{}\left(S^\lambda (1,s)T^\lambda (\mu ,s)\right),$$
(2)
where $`S_\nu ^\lambda (1,s)`$ denote gr-msv  irreducible tensor operators of rank $`\lambda `$ in hadron spin-space. The irreducible tensor amplitudes $`T_\nu ^\lambda (\mu ,s)`$ in Eq. (2) are given explicitly by

$`T_\nu ^\lambda `$ $`(\mu ,s)=\sqrt{{\displaystyle \frac{2}{3}}}(2\pi )(1)^{1s}{\displaystyle \underset{\mathrm{}=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{L=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{j=|\mathrm{}s|}{\overset{(\mathrm{}+s)}{}}}(\text{i})^\mathrm{}L`$
$`\times (1)^{\lambda L}[\mathrm{}][L][j]^2d_{\nu \mu }^L(\theta )C(\mathrm{}L\lambda ;0\nu \nu )`$ (3)
$`\times W(\mathrm{}sL1;j\lambda )\left[T_{L;\mathrm{}s}^{j(\text{mag})}(E)\text{i}\mu T_{L;\mathrm{}s}^{j(\text{elec})}(E)\right],`$
where the energy dependence is carried entirely by the partial wave magnetic and electric $`2^L`$-multipole amplitudes $`T_{L;\mathrm{}s}^{j(\text{mag})}(E)`$ and $`T_{L;\mathrm{}s}^{j(\text{elec})}(E)`$ respectively, while the angular dependence is contained entirely in $`d_{\nu \mu }^L(\theta )`$. We use the shorthand $`[j]=\sqrt{(2j+1)}`$ and the rest of the notations follow rose . The multipole amplitudes are given in terms of the reduced matrix elements by

$$\begin{array}{cc}\hfill T_{L;\mathrm{}s}^{j(\text{mag})}(E)& =\frac{1}{2}\left\{(1)^{L+1}+(1)^{\mathrm{}}\right\}𝑻\hfill \\ \hfill T_{L;\mathrm{}s}^{j(\text{elec})}(E)& =\frac{1}{2}\left\{(1)^L+(1)^{\mathrm{}}\right\}𝑻,\hfill \end{array}$$
(4)
where $`𝑻`$ denotes

$$\begin{array}{cc}\hfill & ||𝑻||=_{I=0,1}(1)^I[I]^{\frac{1}{2}}[1(1)^{\mathrm{}+s+I}]\hfill \\ & \times C(\frac{1}{2}\frac{1}{2}I;\frac{1}{2}\frac{1}{2}0)(1L)j;I_d=0\left|\left|𝑻\right|\right|(\mathrm{}s)j;I.\hfill \end{array}$$
(5)
Expressing the irreducible tensor operators $`S_\nu ^\lambda (1,s)`$ of rank $`\lambda `$ in terms of the Pauli spin-matrices $`𝝈(n)`$ and $`𝝈(p)`$ and unit $`2\times 2`$ matrices $`\sigma _0^0(n),\sigma _0^0(p)`$ for neutron, proton respectively following gr-msv , we may rewrite Eq. (1) in the elegant form

$`{\displaystyle \frac{\text{d}\sigma _0}{\text{d}\mathrm{\Omega }}}={\displaystyle \frac{1}{4}}\text{Tr}\left(𝑴𝑴^{}\right),`$ (6)
where $`\text{Tr}(_{m_d}_\mu )`$ denotes the Trace or Spur and the matrix $`𝑴`$ for $`p(n,\gamma )d`$ has the form

$$𝑴=\underset{\lambda _1,\lambda _2=0}{\overset{1}{}}\underset{\lambda =|\lambda _1\lambda _2|}{\overset{(\lambda _1+\lambda _2)}{}}\left(\left(\sigma ^{\lambda _1}(n)\sigma ^{\lambda _2}(p)\right)^\lambda M^\lambda (\lambda _1,\lambda _2;\mu )\right),$$
(7)
in terms of irreducible tensor amplitudes $`M_\nu ^\lambda (\lambda _1,\lambda _2;\mu )`$ of rank $`\lambda `$ given by

$$\begin{array}{cc}\hfill M_\nu ^\lambda (\lambda _1,\lambda _2;\mu )=& \left(\frac{k}{2\pi E}\right)\left(\frac{E_nE_pE_d}{p}\right)^{\frac{1}{2}}\frac{3}{2}[\lambda _1][\lambda _2]\hfill \\ & \times \underset{s=0}{\overset{1}{}}[s]\left\{\begin{array}{ccc}\frac{1}{2}& \frac{1}{2}& 1\\ \frac{1}{2}& \frac{1}{2}& s\\ \lambda _1& \lambda _2& \lambda \end{array}\right\}T_\nu ^\lambda (\mu ,s),\hfill \end{array}$$
(8)
where $`\{\}`$ denote Wigner-$`9j`$ symbols rus . We may explicit Eq. (7) as

$$\begin{array}{cc}\hfill 𝑴=& A+B\left(𝝈(n)𝝈(p)\right)+\left(𝝈(n)+𝝈(p)\right)𝑪\hfill \\ & +\left(𝝈(n)𝝈(p)\right)𝑫+\left(𝝈(n)\times 𝝈(p)\right)𝑬\hfill \\ & +\left(\left(𝝈(n)𝝈(p)\right)^2F^2\right),\hfill \end{array}$$
(9)
where the coefficients are related to (8) through

$$\begin{array}{cc}\hfill A& =M_0^0(0,0;\mu );B=\frac{1}{\sqrt{3}}M_0^0(1,1;\mu );\hfill \\ \hfill C_\nu ^1& =\frac{1}{2}\left[M_\nu ^1(1,0;\mu )+M_\nu ^1(0,1;\mu )\right];\hfill \\ \hfill D_\nu ^1& =\frac{1}{2}\left[M_\nu ^1(1,0;\mu )M_\nu ^1(0,1;\mu )\right];\hfill \\ \hfill E_\nu ^1& =\frac{\text{i}}{\sqrt{2}}M_\nu ^1(1,1;\mu );F_\nu ^2=M_\nu ^2(1,1;\mu )\hfill \end{array}$$
(10)
and $`C_\nu ^1,D_\nu ^1,E_\nu ^1`$ denote respectively the spherical components of $`𝑪,𝑫,𝑬`$. Comparing Eq. (9) with $`𝑴`$ for elastic $`NN`$ scattering wolf  shows clearly that the fourth and fifth terms containing $`\left(𝝈(n)𝝈(p)\right)`$ and $`\left(𝝈(n)\times 𝝈(p)\right)`$ are the ones which induce transitions from the initial spin-singlet state to the final spin-triplet state of the deuteron. To estimate the singlet and triplet cross sections empirically we observe that the differential cross section for $`\stackrel{}{p}(\stackrel{}{n},\gamma )d`$ is given by

$$\frac{\text{d}\sigma }{\text{d}\mathrm{\Omega }}=\text{Tr}\left(𝑴\rho 𝑴^{}\right),$$
(11)
where the density matrix

$$\rho =\frac{1}{4}\left[1+\left(𝝈(n)𝑷(n)\right)\right]\left[1+\left(𝝈(p)𝑷(p)\right)\right]$$
(12)
describes the initial spin state if $`𝑷(n)`$ and $`𝑷(p)`$ denote respectively the neutron and proton polarizations. Rewriting Eq. (11) as

$$\frac{\text{d}\sigma }{\text{d}\mathrm{\Omega }}=\frac{1}{4}\underset{\alpha ,\beta =0,x,y,z}{}P_\alpha (n)P_\beta (p)B_{\alpha \beta },$$
(13)
where $`P_0(n)=1,P_0(p)=1`$ and

$$B_{\alpha \beta }=\text{Tr}\left(𝑴\sigma _\alpha (n)\sigma _\beta (p)𝑴^{}\right),$$
(14)
we readily see that the unpolarized differential cross section (6) is given by $`(1/4)B_{00}`$. Noting prc  further that

$`𝝅(1,1)={\displaystyle \frac{1}{4}}\left[1+\sigma _z(n)+\sigma _z(p)+\sigma _z(n)\sigma _z(p)\right]`$ (15)
$`𝝅(1,0)={\displaystyle \frac{1}{4}}\left[1+\sigma _x(n)\sigma _x(p)+\sigma _y(n)\sigma _y(p)\sigma _z(n)\sigma _z(p)\right]`$ (16)
$`𝝅(1,1)={\displaystyle \frac{1}{4}}\left[1\sigma _z(n)\sigma _z(p)+\sigma _z(n)\sigma _z(p)\right]`$ (17)
$`𝝅(0,0)={\displaystyle \frac{1}{4}}\left[1\sigma _x(n)\sigma _x(p)\sigma _y(n)\sigma _y(p)\sigma _z(n)\sigma _z(p)\right]`$ (18)
are the projection operators $`|smsm|,s=0,1;m=+s,\mathrm{},s`$, we readily identify

$`{\displaystyle \frac{\text{d}\sigma _{1,1}}{\text{d}\mathrm{\Omega }}}`$ $`={\displaystyle \frac{1}{16}}\left[B_{00}+B_{z0}+B_{0z}+B_{zz}\right]`$ (19)
$`{\displaystyle \frac{\text{d}\sigma _{1,0}}{\text{d}\mathrm{\Omega }}}`$ $`={\displaystyle \frac{1}{16}}\left[B_{00}+B_{xx}+B_{yy}B_{zz}\right]`$ (20)
$`{\displaystyle \frac{\text{d}\sigma _{1,1}}{\text{d}\mathrm{\Omega }}}`$ $`={\displaystyle \frac{1}{16}}\left[B_{00}B_{z0}B_{0z}+B_{zz}\right]`$ (21)
$`{\displaystyle \frac{\text{d}\sigma _{0,0}}{\text{d}\mathrm{\Omega }}}`$ $`={\displaystyle \frac{1}{16}}\left[B_{00}B_{xx}B_{yy}B_{zz}\right]`$ (22)
as the triplet and singlet contributions $`(\text{d}\sigma _{s,m})/(\text{d}\mathrm{\Omega })`$ which add up to $`(\text{d}\sigma _0)/(\text{d}\mathrm{\Omega })`$ given by Eq. (1). Eq. (11) can also be expressed in the form

$`{\displaystyle \frac{\text{d}\sigma }{\text{d}\mathrm{\Omega }}}=`$ $`{\displaystyle \underset{k_n,k_p,k}{}}\left(\left(P^{k_n}(n)P^{k_p}(p)\right)^kB^k(k_n,k_p)\right)`$ (23)
$`=`$ $`{\displaystyle \frac{\text{d}\sigma _0}{\text{d}\mathrm{\Omega }}}[1+𝑷(n)𝑨(n)+𝑷(p)𝑨(p)`$
$`+{\displaystyle \underset{k=0}{\overset{2}{}}}((𝑷(n)𝑷(p))^kA^k)],`$ (24)
where the irreducible bilinear amplitudes in Eq. (23)

$`B_q^k(k_n,k_p)={\displaystyle \frac{3}{2}}\left({\displaystyle \frac{k}{2\pi E}}\right)^2{\displaystyle \frac{E_nE_pE_d}{p}}(1)^{k_n+k_pk}[k_n][k_p]`$
$`\times {\displaystyle \underset{s,s^{},\lambda ,\lambda ^{},\mu }{}}(1)^\lambda [s][s^{}][\lambda ][\lambda ^{}]W(s^{}k1\lambda ;s\lambda ^{})`$
$`\times \left\{\begin{array}{ccc}\frac{1}{2}& \frac{1}{2}& s\\ \frac{1}{2}& \frac{1}{2}& s^{}\\ k_n& k_p& k\end{array}\right\}\left(T^\lambda (\mu ,s)T^^\lambda ^{}(\mu ,s^{})\right)_q^k`$ (25)
are related to the spherical components of the analyzing powers in Eq. (Radiative capture of polarized neutrons by polarized protons) through

$`{\displaystyle \frac{\text{d}\sigma _0}{\text{d}\mathrm{\Omega }}}A_q^1(n)=B_q^1(1,0)`$ (26)
$`{\displaystyle \frac{\text{d}\sigma _0}{\text{d}\mathrm{\Omega }}}A_q^1(p)=B_q^1(0,1)`$ (27)
$`{\displaystyle \frac{\text{d}\sigma _0}{\text{d}\mathrm{\Omega }}}A_q^k=B_q^k(1,1),`$ (28)
with $`B_0^0(0,0)=(\text{d}\sigma _0)/(\text{d}\mathrm{\Omega })`$, the unpolarized differential cross section. From Eq. (Radiative capture of polarized neutrons by polarized protons), we observe that the $`T_\nu ^\lambda (\mu ,s)`$ satisfy the symmetry property

$`T_\nu ^\lambda (\mu ,s)=(1)^{\lambda \nu }T_\nu ^\lambda (\mu ,s).`$ (29)
This in turn implies through Eq. (Radiative capture of polarized neutrons by polarized protons) that the bilinear amplitudes satisfy the property

$`B_q^k(k_n,k_p)=(1)^{kq}B_q^k(k_n,k_p)`$ (30)
so that

$`B_0^1(k_n,k_p)=0`$ (31)
and from Eqs. (26) and (27),

$`A_z(n)=A_z(p)=0.`$ (32)
Using these results, we can now write the channel-spin differential cross sections $`(\text{d}\sigma _{s,m})/(\text{d}\mathrm{\Omega })`$ in terms of the analyzing powers as

$`{\displaystyle \frac{\text{d}\sigma _{1,1}}{\text{d}\mathrm{\Omega }}}`$ $`={\displaystyle \frac{\text{d}\sigma _{1,1}}{\text{d}\mathrm{\Omega }}}={\displaystyle \frac{1}{4}}{\displaystyle \frac{\text{d}\sigma _0}{\text{d}\mathrm{\Omega }}}\left[1\frac{1}{\sqrt{3}}A_0^0+\sqrt{\frac{2}{3}}A_0^2\right]`$ (33)
$`{\displaystyle \frac{\text{d}\sigma _{1,0}}{\text{d}\mathrm{\Omega }}}`$ $`={\displaystyle \frac{1}{4}}{\displaystyle \frac{\text{d}\sigma _0}{\text{d}\mathrm{\Omega }}}\left[1\frac{1}{\sqrt{3}}A_0^02\sqrt{\frac{2}{3}}A_0^2\right]`$ (34)
$`{\displaystyle \frac{\text{d}\sigma _{0,0}}{\text{d}\mathrm{\Omega }}}`$ $`={\displaystyle \frac{1}{4}}{\displaystyle \frac{\text{d}\sigma _0}{\text{d}\mathrm{\Omega }}}\left[1+\sqrt{3}A_0^0\right].`$ (35)
The analyzing powers, moreover, are measurable readily in an experimental setup such as in muller , where $`𝑷(n)=𝑷(p)=P\widehat{𝑷}`$ with $`P=50.5\%`$ and $`\widehat{𝑷}`$ could be longitudinal as well as transverse. With $`\widehat{𝑷}`$ chosen parallel and antiparallel to the beam, it is clear that

$`{\displaystyle \frac{\text{d}\sigma _{+z}}{\text{d}\mathrm{\Omega }}}+{\displaystyle \frac{\text{d}\sigma _z}{\text{d}\mathrm{\Omega }}}=2{\displaystyle \frac{\text{d}\sigma _0}{\text{d}\mathrm{\Omega }}}\left[1+{\displaystyle \frac{P^2}{\sqrt{3}}}\left(\sqrt{2}A_0^2A_0^0\right)\right]`$ (36)
which yields $`\left(\sqrt{2}A_0^2A_0^0\right)`$.

Likewise, if $`\widehat{𝑷}`$ is chosen parallel and antiparallel to a direction (say $`\widehat{x}`$ of $`\widehat{y}`$) perpendicular to the beam, we have

$`{\displaystyle \frac{\text{d}\sigma _{+x}}{\text{d}\mathrm{\Omega }}}+{\displaystyle \frac{\text{d}\sigma _x}{\text{d}\mathrm{\Omega }}}`$ $`=2{\displaystyle \frac{\text{d}\sigma _0}{\text{d}\mathrm{\Omega }}}\left[1{\displaystyle \frac{P^2}{\sqrt{3}}}\left(A_0^0+\frac{1}{\sqrt{2}}A_0^2\right)+A_2^2\right]`$ (37)
$`{\displaystyle \frac{\text{d}\sigma _{+y}}{\text{d}\mathrm{\Omega }}}+{\displaystyle \frac{\text{d}\sigma _y}{\text{d}\mathrm{\Omega }}}`$ $`=2{\displaystyle \frac{\text{d}\sigma _0}{\text{d}\mathrm{\Omega }}}\left[1{\displaystyle \frac{P^2}{\sqrt{3}}}\left(A_0^0+\frac{1}{\sqrt{2}}A_0^2\right)A_2^2\right],`$ (38)
where we have made use of the fact that $`A_2^2=A_2^2`$ from Eqs. (28) and (30). By adding Eqs. (37) and (38), one readily obtains $`\left(A_0^0+\frac{1}{\sqrt{2}}A_0^2\right)`$. Since $`\left(\sqrt{2}A_0^2A_0^0\right)`$ is already known from Eq. (36), one can determine $`A_0^0`$ as well as $`A_0^2`$. Hence all the $`(\text{d}\sigma _{s,m})/(\text{d}\mathrm{\Omega })`$ for $`m=s,\mathrm{},s`$ and $`s=0,1`$ are determinable from experiment empirically at any given energy. We may perhaps add that these measurements need not have to be carried out at accuracies of $`10^7`$ to determine $`(\text{d}\sigma _{s,m})/(\text{d}\mathrm{\Omega })`$. Moreover, the discussion presented here aims to extend the program of Breit and Rustgi breit-rustgi  to actual individual determinations of the three triplet and one singlet cross sections for the important astro  fusion process at the differential level itself rather that estimate the relative importance of the triplet amplitude vis-a-vis the singlet amplitude. It would, therefore, be desirable to extend the recent experiment muller  on $`\stackrel{}{p}(\stackrel{}{n},\gamma )d`$ to measure the observables (36) to (38) in order to determine the differential cross sections $`(\text{d}\sigma _{s,m})/(\text{d}\mathrm{\Omega })`$ individually and hence study the fusion reaction more incisively at any given energy.

It is perhaps not out of place to mention here that the highly successful $`NN`$ potential models like Nijmegen1 and Nijmegen2 potentials Nijmegen , the Argonne potential Argonne  and the Bonn potential Bonn  which reproduce the elastic $`NN`$ scattering data with a remarkable $`\chi ^2/\text{datum}1`$, have been subjected to an interesting study of the deuteron properties by Polls et al. polls , who say, It is well known that the off-shell behavior of NN potentials cannot be pinned down by on-shell data…Moreover, within a given model, pinning down the on-shell point limits the range of variation for the off-shell behavior. The four modern high-precision models considered in this study pin down the on-shell $`T`$-matrix as much as by all means possible since they fit the NN scattering data with the perfect $`\chi ^2`$/ datum $``$ 1. Furthermore, it may be reasonable to believe that the four models cover about the range of diversity that there is to realistic physical models for the NN interaction. Based upon these premises, one may then conclude that the off-shell uncertainties revealed in this study are what we have to deal with, at the current status of theoretical nuclear physics. At this time, we do not know of any other objectively verifiable aspects that could further reduce the off-shell uncertainties. Tighter constraints for the off-shell behavior of NN may emerge in the future… In this context, it is perhaps pertinent to point out that $`npd\gamma `$ involves necessarily off-shell $`NN`$ interactions. Therefore it would also be of interest to examine the fusion reaction with a view to throw light on the off-shell uncertainties. Further work is in progress.

###### Acknowledgements.

One of us (GR) thanks Professor B. V. Srikantan for much encouragement and Professor Ramanatha Cowsik for inviting him to the Indian Institute of Astrophysics, while the other (PND) thanks the Council of Scientific and Industrial Research (CSIR), India for support through the award of a Senior Research Fellowship.


text:
# Exchange couplings in the magnetic molecular cluster Mn12Ac

## I Introduction

Molecular nanomagnets consist of a few paramagnetic ions coupled by exchange interactions and their properties are lying at the borderline between quantum and classical behavior. The best studied cluster Lis  so far is certainly Mn12Ac. At low fields and low temperature its ground state can be described as a single spin S=10. The reversal of this spin occurs via a macroscopic relaxation time at T=2 K. This discovery Sessoli93  has prompted a detailed examination of its magnetic behavior because of the potential use for ultimate information storage. If we concentrate only on magnetic ions, then Mn12Ac is made of an external ring of eight Mn3+ ions with spin S=2 and ring encloses a tetrahedron of four Mn4+ ions with spin S=3/2. These spins are coupled by exchange so as to lead to a S=10 ground state. This means that there is a ferrimagnetic arrangement of the spins. It is also known that the S=10 manifold is split by anisotropy into sublevels with $`10S^z+10`$. The properties of this subset of levels have been studied in detail by various experimental techniques because of the possibility of macroscopic quantum tunneling of the spin. The anisotropy responsible for the zero-field splitting of the S=10 manifold is smaller than the exchange interactions that determine first of all the ferrimagnetic ground state structure. For example the degenerate $`S^z=\pm 10`$ states are separated from the states with $`S^z=\pm 9`$ (belonging to the same S=10 manifold) by a gap determined by inelastic neutron scattering MirebeauL99  $`14`$ K. So a reasonable strategy to study magnetic properties of this cluster is to ignore first any anisotropy terms in a model spin Hamiltonian by use of the simplest Heisenberg exchange and then, once these parameters are known, refine the treatment by including higher-order anisotropy terms.

Recently the full magnetization curve M(H) of the cluster Mn12Ac has been obtained explosive  by an experimental technique using explosive compression of the magnetic flux to access to very high fields in the megagauss range. Many if not all transitions between levels with different total spin have been measured. This kind of measurement is a direct probe of the inner magnetic structure of the cluster i.e. of the higher-energy scale couplings.

In this paper we determine the values of the exchange couplings between the Mn ions by using a Heisenberg spin Hamiltonian to reproduce the high-field magnetization data. We also add as a constraint the value of the gap between S=10 and S=9. With these values we are able to predict the energies and multiplicities of some of the low-lying states above the S=10 ground state. This explains some features of the existing inelastic neutron scattering data that were not understood previously. In section II, we expose the details of the model Hamiltonian and the numerical technique we have used. In section III, the details of the fit are given. Finally section IV contains our conclusions.

## II The spin Hamiltonian

The core of the magnetic cluster of formula Mn12O12(CH3COO)16(H2O)4 is depicted in Fig. 1. There is an outer ring of eight ions Mn3+ with spin S=2, $`𝐒_5,\mathrm{}𝐒_{12}`$ that encircles four ions Mn4+ with S=3/2, $`𝐒_1,\mathrm{}𝐒_4`$. In Fig. 2 we draw a simplified representation of the cluster taking into account only the connectivity of exchange interactions Caneschi91 ; Sessoli93 . There are four most important exchange paths (bonds in the following). They are J1 which relate each spin 3/2 to a spin S=2, there is also a J2 coupling between the outer ring and the inner tetrahedron involving four triangles i.e. eight bonds, then J3 between the four spins S=3/2 with a tetrahedral structure so there are six such bonds, and J4 between the S=2 ions which means eight bonds. The simplest guess for the exchange is thus a Heisenberg spin Hamiltonian :

$`=`$ $`J_1`$ $`\{𝐒_1𝐒_7+𝐒_2𝐒_9+𝐒_3𝐒_{11}+𝐒_4𝐒_5\}`$
$`+`$ $`J_2`$ $`\{𝐒_6(𝐒_1+𝐒_4)+𝐒_8(𝐒_1+𝐒_2)+𝐒_{10}(𝐒_2+𝐒_3)+𝐒_{12}(𝐒_3+𝐒_4)\}`$
$`+`$ $`J_3`$ $`\{𝐒_1𝐒_2+𝐒_2𝐒_3+𝐒_3𝐒_4+𝐒_4𝐒_1+𝐒_1𝐒_3+𝐒_2𝐒_4\}`$
$`+`$ $`J_4`$ $`\{𝐒_5𝐒_6+𝐒_6𝐒_7+𝐒_7𝐒_8+𝐒_8𝐒_9+𝐒_9𝐒_{10}+𝐒_{10}𝐒_{11}+𝐒_{11}𝐒_{12}+𝐒_{12}𝐒_5\},`$ (1)
where we use the labeling of the spins given in Fig. 2. In this work we ignore the effect of anisotropy and use the Hamiltonian Eq.(1). As a consequence the total spin is a good quantum number. The total dimension of the Hilbert space of the magnetic degrees of freedom is $`10^8`$ which is too huge for brute force diagonalization. Previous works Sessoli93 ; Barrat97 ; Zvezdin96  have used various approximate schemes usually based on some assumptions about the relative order of magnitude of the J’s. In the Florentine coupling scheme Sessoli93  one starts from a large antiferromagnetic value for J1. Hence the cluster is described as a first approximation by four dimers Mn3+-Mn4+ with spin S=1/2 and four remaining Mn3+ ions with spin S=2. It is then feasible to treat the system of four spins S=1/2 and four spins S=2. However it is not clear that this hypothesis is able to explain the experimental results at hand. There are many results from neutron scattering MirebeauL99 , magnetization Barbara98  and heat capacity Gomes98  that point to the presence of excited states not too far from the S=10 ground state. It is thus desirable to obtain values of the exchange parameters without recourse to any a priori assumption of their relative strength. A recent study by exact diagonalization Raghu2001  using the point symmetry group of the magnetic cluster concluded that J1=215 K, J2=85 K, J3=85 K, J4=-64.5 K. Since these data refer to transitions involving the lowest energy levels for a given spin value, the Lánczos algorithm appears to be particularly well suited.

The recent megagauss magnetization data explosive  provide further information on the strength of the exchange interactions. We have used the Lánczos algorithm Golub  applied to the Hamiltonian Eq.(1). If we use conservation of spin projection on a given axis, say it z, then the Hilbert space dimensions are much reduced. The largest subspace corresponds to $`S^z=0`$ with dimension 8581300, while the subspace in which we have to find the ground state has dimension 817176 for $`S^z=10`$. The corresponding dimensions are given in table 1. Although the use of total spin would reduce further the dimensionality, it is much more complicated to program in an efficient way.

In a given subspace we use the iterative Lánczos algorithm that brings the Hamiltonian into tridiagonal form. This is done typically in at most one hundred iterations which is enough to obtain the ground state energy with a precision of $`10^{10}`$. The advantage of the Lánczos algorithm is that it requires only to perform the product of the Hamiltonian on a vector. Since spin Hamiltonians are very sparse matrices, one needs only to store in fact two vectors to use the Lánczos algorithm. Even with the huge dimensions that appear in table I, this is feasible on present day computers.

If we now consider the effect of a magnetic field Golinelli , then due to the rotational symmetry of the Hamiltonian, the magnetic field couples to a conserved quantity and we just have to shift energies :

$$E_{S^z}(B)=E_{S^z}(B=0)g\mu _\mathrm{B}BS^z,$$
(2)
where $`E_{S^z}(B=0)`$ is an eigenenergy of (1) in the sector with spin projection $`S^z`$. As a function of the applied field there will be crossings of levels with different values of $`S^z`$. The magnetization curve of the cluster Mn12Ac is a series of discrete jumps and plateaus due to the finite size of this spin system. It is only for the thermodynamic limit that one gets smooth magnetization curves for spin systems. The critical field corresponding to the transition $`S^z`$ to $`S^z+1`$ is given by :

$$B_{crit}(S^zS^z+1)=\frac{1}{g\mu _\mathrm{B}}(E_{0,S^z+1}(B=0)E_{0,S^z}(B=0)),$$
(3)
where $`E_{0,S^z}(B=0)`$ is the ground state energy in the sector with given $`S^z`$. We have used g=2 in this paper.

## III Reproducing Magnetization data

To determine the exchange parameters, we have used the data coming from various experimental sources. First it is well established that the ground state has total spin S=10. These states are exactly degenerate in our case since we make the simplifying assumption of exact rotational symmetry. Then the gap between this S=10 manifold and the first excited state with S=9 is known from magnetic susceptibility measurements Caneschi92  to be 35 K. The remaining piece of high-energy information comes from the megagauss experiment of ref. explosive . In this experiment one measures the differential susceptibility $`dM/dH`$ vs. $`H`$. For a system with a series of discrete jumps as Mn12Ac, a spike in the differential susceptibility corresponds to a change of the spin of the ground state. The lowest-lying spike lies at $`B_1=`$ 382 T and is interpreted Zvezdin96  as the transition from the ground state with S=10 to an S=11 state and then there are three spikes at $`B_2=`$ 416 T, $`B_3=`$ 448 T, and $`B_4=`$ 475 T corresponding to the crossing of S=12, 13, and 14 states respectively. At higher fields there is a huge spike centered at 530 T which presumably corresponds to several unresolved crossings, maybe S=15 and S=16. Above this field value it is difficult to locate the other remaining transitions so no other values were determined in Ref. explosive . We have used the four values that are determined with good accuracy $`B_1,\mathrm{},B_4`$. We compute the corresponding theoretical values $`B_i^{\mathrm{theo}}`$ and measure the quality of the fit by the following quantity :

$$ϵ=\underset{i=1}{\overset{4}{}}\left(\frac{B_iB_i^{\mathrm{theo}}}{B_i}\right)^2.$$
(4)
We have performed calculations of the levels with Sz=9 up to Sz=22 on a grid of values of the ratios of the exchange couplings J2/J1, J3/J1 and J4/J1 between +2 and -2 by steps of 0.1. In this range of parameters, we first reject values for which Sz=10 is not the ground state. Then we look for regions with small $`ϵ`$ parameter defined in Eq.(4). In these regions, we fix the absolute scale of energy via J1 by requiring that the critical fields should be equal to the experimental value $`B_1=`$ 382 T. This leads to approximate values for the other critical fields $`B_2,\mathrm{},B_4`$. We next refine the search by including only the regions in which the gap S=10-S=9 is close to its experimental value of 35 K. This is now done by adjusting the four dimensionfull couplings J1,…,J4. We observe that the biggest effect on the overall spectrum structure is due to J1 and J2. The best values are close to J1 $`119`$ K, J2 $`118`$ K. This corresponds to $`ϵ10^4`$. In Fig. 3 we plot the variation of the fitting parameter $`ϵ`$ in the J1-J2 plane close to the best values of these parameters. If we search for a region of small $`ϵ`$ by tuning these two parameters only then we find that this region is quite insensitive to the choice of J3 and J4. This is illustrated in Fig. 4 where $`ϵ`$ is given in the J3-J4 plane.

Our next observation is that once J1 and J2 are determined, the remaining couplings J3 and J4 may be varied to obtain more precisely the S=10-S=9 gap of 35 K. The agreement with the gap is given in Fig. 5 where we have plotted the gap value in the plane J3-J4. We find that J4 is significatively antiferromagnetic J$`{}_{4}{}^{}23`$ K while it is difficult to give a precise estimate for J3 : it barely differs significatively from zero. A tentative value is J$`{}_{3}{}^{}8`$ K. Taking into account the experimental uncertainties on the critical fields, we estimate that J3 and J4 are determined with an error bar of $``$ 6 K. The determination of J3 and J4 does not affect much the values of J1 and J2, as can be seen in Fig. 6.

As proposed in previous works, we find that J1 is larger than J3 and J4. However, we clearly need a second coupling J2 which should be close to J1 (taking into account the uncertainties on all the experimental results). The part of the spectrum relevant to magnetization data corresponding to this preferred set of parameters is given in Table II.

It is interesting to compare our results with those already available in the literature. In the first works Sessoli93 , the assumption of large J1 combined with a perturbative treatment lead to J1=225 K, J2=90 K, J3=90 K, J4=0 K. In fact this set of parameters does not lead to a S=10 ground state, as first pointed out by Raghu et al. Raghu2001 . There is another set of parameters J1=215 K, J2=85 K, J3=-85 K, J4=-45 K suggested by Chudnovsky Chudnovsky96 . However this set which has a correct S=10-S=9 ordering leads to a gap of 223 K, much too large. A recent study by exact diagonalization Raghu2001  concluded that J1=215 K, J2=85 K, J3=85 K, J4=-64.5 K. This set has a correct ordering of 9-10 levels with a gap which is adjusted to the experimental value of 35 K but we find that it does not lead to a satisfactory magnetization curve. This set gives critical fields equal to $`B_1=`$ 192 T, $`B_2=`$ 239 T, $`B_3=`$ 356 T and $`B_4=`$ 406 T (the $`ϵ`$ parameter is 0.4 instead of our value of $`10^4`$).

With our preferred set of exchange parameters we can then compute some of the excited levels that are above the first excited manifold with S=9. We find that there is a S=8 manifold sitting at 37 K above the S=10 ground state so very close to the S=9 states that have been used to constrain our fit. Then there is a S=9 multiplet which is found at 67 K above the ground state. This picture is close to what is found from susceptibility measurements Caneschi92  where two close S=9 levels are necessary to reproduce the low-temperature behavior. However the scheme is not exactly the same : the multiplicities do not coincide.

Our new picture has some interesting consequences for the interpretation of neutron data. There is some evidence from inelastic neutron scattering Hennion97  for a mode at 1.2 THz i.e. 70 K, this may correspond to the second S=9 multiplet that we find at 67 K. This is allowed by neutron scattering selection rule starting from the S=10 ground state. It is then normal that neutrons do not see the S=8 states at 37 K = 0.73 THz because of the selection rules however it not yet clear why neutrons do not see the first excited state with S=9 at 35 K = 0.72 THz.

## IV Conclusions

We have studied the energy levels of a Heisenberg spin model Hamiltonian appropriate to describe the magnetic cluster Mn12Ac. We obtain a determination of four exchange couplings : J1=119 K, J2=118 K, J3=-8 K, J4=23 K. Such a set of parameters reproduces the magnetization curve observed in megagauss experiments, leads to a S=10 ground state and a gap of 35 K to a first excited level with S=9 as measured experimentally. The numerical method we have used does not rely upon any approximations, so the main source of uncertainty in the values we quote comes from the measurement of the critical fields in the magnetization process. We predict that there is an excited level with S=8 at 37 K only slightly above the S=9 excited state which lies at 35 K and the next excited state is a S=9 multiplet at 67 K above the S=10 ground state. These findings explain part of the existing neutron scattering data.

###### Acknowledgements.

We wish to acknowledge support from the European community (Molnanomag HPRN-CT 1999 0012), the French Ministry of Research and the Region Ile-de-France (Chaire Blaise Pascal 2001). We thank also Andreas Honecker and Bruce Normand for a careful reading of the manuscript.


text:
# Charge and stability of strange quark matter at finite temperature

## I Introduction

It has been found from the theoretical study that the introduction of strange quarks into a plasma with two flavors could lower the Fermi energy of the system and thus the mass of quark matter . If its mass is lower than the mass of hyperonic matter with the same strangeness fraction, or even lower than the non-strange nucleonic matter, the strange quark matter (SQM) would be the true ground state of QCD. Therefore, to search the SQM in the laboratory is an important task for recent heavy ion collision experiments at ultrarelativistic energy .

The essential problem for detectability of SQM in heavy ion experiments is to study its stability during the formation of quark gluon plasma (QGP). Many properties of stability of SQM, for example, the ”stability window” , the charge-dependence  and the strangeness fraction-dependence  of stability have been widely discussed by many authors. They employed different models to investigate this problem. But unfortunately, many conclusions are different, in particular, part of them, say, the charge dependence is contrary.

The electric charge is of vital importance to the experimental searches of SQM. By means of MIT bag model, Jaffe and his co-workers  suggested that the SQM is slightly positively electric charged. They considered strong and weak decay by nucleon and hyperon emission together and concluded that the stable SQM will have a low but positive charge to mass ratio. Using the same model as that of Jaffe et al. but considering the initial condition of possible strangelet production in relativistic heavy ion collision carefully, Greiner and his co-workers  argued that strangelets are most likely highly negatively charged. They also emphasized that the stability of infinite SQM depends on the bag constant $`B`$. An infinite SQM, treated as a gas of noninteracting quarks, is absolutely stable for $`B^{1/4}145\text{MeV}`$, and metastable for $`B^{1/4}150200\text{MeV}`$, but unstable for larger bag constants. Of course, this conflict of conclusions is limited in the framework of MIT bag model and at zero temperature.

Employing another model, namely, the quark mass density-dependent (QMDD) model, Peng and his co-workers reconsidered this problem . They introduced a critical density and proved that SQM cannot maintain its flavor equilibrium, and thus no longer be stable under this critical density. They found that proper negative charges can lower the critical density and thus make the SQM become more stable. Their conclusion is not surprised if one notice the confinement mechanism of QMDD model. In QMDD model, the mass of $`u,d`$ and $`s`$ quarks are given by

$`m_q`$ $`=`$ $`{\displaystyle \frac{B}{3n_B}},(q=u,d,\overline{u},\overline{d}),`$ (1)
$`m_{s,\overline{s}}`$ $`=`$ $`m_{s0}+{\displaystyle \frac{B}{3n_B}},`$ (2)
where $`n_B`$ is the baryon number density, $`m_{s0}`$ is the current mass of the strange quark matter and $`B`$ is the vacuum energy density inside bag (bag constant). The mechanism of confinement can be mimicked through the hypothesis Eqs.(1) and (2). The masses of quarks become infinitely large as the volume increases to infinity or the density decreases to zero. So that the vacuum is unable to support it . Noting that the boundary condition of confinement of MIT bag model corresponds to the zero quark mass inside the bag but infinity at the boundary or outside the bag, we see that the confinement mechanisms are similar to these two models. Due to this correspondence, as pointed out by Benrenuto and Lugones , in almost all cases, the properties of SQM given by QMDD model are nearly the same as those obtained in the MIT bag model. The conclusion of Peng et al. supports the result given by Greiner et al.

This paper evolves from an attempt to extend these discussions to finite temperature. The reasons of this extension are: at first, the SQM might be possible and formed during the phase transition from hadronic matter to a deconfined QGP. From lattice QCD calculations, the phase transition of deconfinement will happen at critical temperature. Secondly, as was pointed out by our previous papers , QMDD model cannot be used to explain the process of the quark deconfinement phase transition because the quark confinement is permanent in this model. Instead of this model, a quark mass density- and temperature-dependent (QMDTD) model in which the quark confinement is not permanent has been suggested by us . As its application, we hope to employ QMDTD model to discuss the charge and strangeness fraction dependences of SQM at finite temperature.

The organization of this paper is as follows: In next section, we review the QMDTD model and give the main formulae which are necessary for studying the thermodynamical behaviors of SQM. In section 3, we give the result of numerical calculations for QMDTD model and investigate the stability of hot SQM in detail. Our conclusions are summarized in the last section.

## II Quark mass density- and temperature-dependent models

To overcome the difficulties of QMDD model which appear at finite temperature, in refs., based on the Friedberg-Lee model, we suggested a QMDTD model. The basic difference between our model and QMDD model is that, instead of a constant $`B`$ in Eqs.(1) and (2), we argued that $`B`$ must be a function of temperature and choose $`B(T)`$ as

$`B`$ $`=`$ $`B_0\left[1\left({\displaystyle \frac{T}{T_c}}\right)^2\right],0TT_c`$ (3)
$`B`$ $`=`$ $`0,\text{ }T>T_c,`$ (4)
where $`B_0`$ is bag constant at zero temperature and $`T_c`$ is the critical temperature of quarks deconfinement. Substituting Eq.(3) into Eqs.(1) and (2), we find the masses of $`u,d,s`$ quarks and their corresponding anti-quarks as

$`m_q`$ $`=`$ $`{\displaystyle \frac{B_0}{3n_B}}\left[1\left({\displaystyle \frac{T}{T_c}}\right)^2\right],(q=u,d,\overline{u},\overline{d}),`$ (5)
$`m_{s,\overline{s}}`$ $`=`$ $`m_{s0}+{\displaystyle \frac{B_0}{3n_B}}\left[1\left({\displaystyle \frac{T}{T_c}}\right)^2\right],`$ (6)
when $`0TT_c`$. When $`TT_c`$, the masses of quarks are independent of temperature: $`m_{u,d,\overline{u},\overline{d}}=0`$ and $`m_{s,\overline{s}}=m_{s0}`$. Due to the difference between Eqs.(1),(2) and Eqs.(5),(6), as was proved in ref., we found that when $`n_B`$ approaches to zero, the temperature $`T`$ $`T_c`$ in our model, but diverges in QMDD model. In QMDTD model, the masses of $`u,d`$ quarks become zero (or constant for strange quarks) when the temperature approaches $`T_c`$. This means that the quarks can be deconfined. Here, the vacuum energy density $`B`$ is a decreasing function as temperature increases. As was pointed out by Greiner et al., the chosen values of $`B`$ affect the stability of SQM remarkably . Therefore, it is of interest to consider the effect of $`B(T)`$.

Treating the infinite SQM as a gas of noninteracting quarks and electrons , we find that the thermodynamical potential density of SQM is

$$\mathrm{\Omega }=\underset{i}{}\mathrm{\Omega }_i=\underset{i}{}\frac{g_iT}{(2\pi )^3}_0^{\mathrm{}}𝑑k\frac{dN_i}{dk}\mathrm{ln}\left(1+e^{\beta (\epsilon _i(k)\mu _i)}\right),$$
(7)
where $`i`$ stands for $`u,d,s`$ (or $`\overline{u},\overline{d},\overline{s}`$ ) and the electron $`e`$($`e^+`$), $`g_i=6`$ for quarks and antiquarks, $`g_i=2`$ for $`e`$ and $`e^+`$. $`\epsilon _i(k)=\sqrt{m_i^2+k^2}`$ is the single particle energy and $`m_{i\text{ }}`$is given by Eqs.(5),(6) for quarks and antiquarks. $`\mu _i`$ is the chemical potential (for antiparticle $`\mu _{\overline{i}}=\mu _i`$). The number density of each particle can be obtained from eq.(7) by means of

$$n_i=\frac{1}{V}\frac{\mathrm{\Omega }}{\mu _i}|_{T,n_B},$$
(8)
and we find

$$\mathrm{\Delta }n_i=n_in_{\overline{i}}=\frac{g_i}{(2\pi )^3}_0^{\mathrm{}}d^3k\left(\frac{1}{\mathrm{exp}[\beta (\epsilon _i\mu _i)]+1}\frac{1}{\mathrm{exp}[\beta (\epsilon _i+\mu _i)]+1}\right).$$
(9)
For SQM, the baryon number density satisfies

$$n_B=\frac{1}{3}(\mathrm{\Delta }n_u+\mathrm{\Delta }n_d+\mathrm{\Delta }n_s).$$
(10)
The electric charge density $`Q`$ of SQM reads

$$Q=\frac{2}{3}\mathrm{\Delta }n_u\frac{1}{3}\mathrm{\Delta }n_d\frac{1}{3}\mathrm{\Delta }n_s\mathrm{\Delta }n_e,$$
(11)
and the light quarks are converted to strange quarks through the weak processes

$$u+du+s,su+e^{}+\overline{\nu }_e,du+e^{}+\overline{\nu }_e,u+e^{}d+\nu _e.$$
(12)
Neglecting the contribution of neutrinos , we obtain the conditions of chemical equilibrium from the reactions (12) as

$`\mu _s`$ $`=`$ $`\mu _d,`$ (13)
$`\mu _s`$ $`=`$ $`\mu _u+\mu _e.`$ (14)
Noting that the chemical potentials of particles and anti-particles are determined by Fermi distribution Eq.(9), we see that Eqs.(5), (6), (9), (10), (11), (13) and (14) form an equation group and can be solved self-consistently. The solution of this equation group determines the chemical stability of SQM. The results of numerical calculation will show in the next section. Here we want to point out that at a fixed temperature $`T`$ and for a fixed value of charge density $`Q`$ this equation group which determines the configuration of SQM has no solution under a critical baryon number density $`n_{Bc}`$. This means that the SQM is unstable when $`n_B<n_{Bc}`$. With the values of $`Q`$ changed and the temperatures limited between $`0TT_c`$, the condition $`n_Bn_{Bc}`$ determines a stable area which might be useful for the detectability of SQM in experiments of high energy heavy ion collisions.

## III Conclusion and discussion

The numerical calculations of the equation group mentioned above have been done with the parameters set:

$$B_0=170\text{MeV}\text{fm}^3,m_{s0}=150\text{MeV},T_c=170\text{MeV}.$$
(15)
Firstly, we study a SQM system with neutrality charge $`Q=0`$, for a given $`n_B`$, the fractions of different quarks are given by

$$F_{q(q=u,d,s)}=\frac{\mathrm{\Delta }n_q}{3n_B}.$$
(16)
At a fixed temperature $`T=100\text{MeV}`$, the $`u,d,s`$ quarks fraction are drawn in Fig.1 by solid line, dashed line and dotted line, respectively. At high baryon number density, all of the $`u,d,s`$ quarks tend to become a triplicate. When the density becomes lower, $`u`$ fraction increases and $`s`$ fraction decreases monotonously. We see from this figure that the fractions of $`u,d,s`$ quarks satisfy an inequality:

$$F_d>F_u>F_s.$$
(17)
In particular, we would like to emphasize that, as shown in Fig.1, three different quark fraction curves stop at a critical baryon number density ($`n_{Bc}=8.28\times 10^3\text{fm}^3`$). One can not find real roots for $`F_d,F_u`$ and $`F_s`$ provided $`n_B<n_{Bc}`$. At $`n_{Bc}`$, $`F_s=0.1070`$. In ref., the authors defined the critical density by the condition $`F_s=0`$. Of course, the SQM can not maintain the chemical equilibrium if no strange quarks exist. But as shown in Fig.1, this is only a sufficient but is not the necessary condition. Because, one cannot find solution in the ranges $`0<F_s<0.107`$.

To compare the QMDD model in which $`B`$ is a constant and our QMDTD model in which $`B`$ satisfies Eqs.(3) and (4), we plot the critical density $`n_{Bc}`$ vs temperature $`T`$ with neutral charge density $`Q=0`$ in figure 2, where the dashed line and the solid line refer to QMDD model and QMDTD model, respectively. We see that the curve given by QMDD model can never be zero even the temperature increases to $`T>T_c`$. It can easily be understood if we notice that $`m_q`$ approaches to infinity when $`n_B0`$ from Eqs.(1) and (2). To excite a particle with infinite mass must spend infinite energy and then infinite temperature. Because of $`n_{Bc}0`$, when $`n_B<n_{Bc}`$, the system would still be unstable even at very high temperature for QMDD model. But for QMDTD model, we see from the solid line that $`n_{Bc}=0`$ at $`T=T_c=170\text{MeV}`$. It means that the SQM is completely stable and may be detected from high energy heavy ion collision.

Now we are in a position to address the charge dependence of stability for SQM. The critical density vs temperature curves with various charges are shown in figure 3, where the solid line, dotted line and dashed line refer to $`Q=0,0.1e\text{fm}^3,0.2e\text{fm}^3`$, respectively. We see that the critical density $`n_{Bc}`$ decreases rapidly and the stable area expands widely when $`Q`$ decreases to more negative. So the negative charge is in favor of the stability of SQM. Our result supports the conclusion given by Greiner et al.  and extends their result to finite temperature.

To study the effect of the strangeness fraction, we show the same curves as that of Fig.1 but for $`Q=0.2e\text{fm}^3`$ in figure 4. Comparing Fig.1 with Fig.4, we find that all of three curves tend to become a triplicate for $`Q=0.2e\text{fm}^3`$ at high density, too. And there are two major differences between these two figures, firstly, the critical density decreases from $`8.28\times 10^3\text{fm}^3`$ to $`5.02\times 10^4\text{fm}^3`$ when charge density $`Q`$ changes from neutrality to negativity $`0.2e\text{fm}^3`$; secondly, Eq.(17) becomes

$$F_d>F_s>F_u.$$
(18)
The fraction of $`s`$ quarks is higher than the fraction of $`u`$ quarks. This is of course reasonable because $`s`$ quark has negative charge and $`u`$ quark positive. The critical density $`n_{Bc}`$ decreases as the $`F_s`$ increases. The higher strangeness fraction is in favor of the stability of SQM. To illustrate this result more transparently, using the expression of strangeness fraction

$$f_s=\frac{\mathrm{\Delta }n_s}{n_B},$$
(19)
which is usually employed by many references, we draw the critical strangeness fraction $`f_{sc}`$ vs temperature $`T`$ curves in figure 5, where the solid line and dashed line refer to $`Q=0`$ and $`0.2e\text{fm}^3`$, respectively. The $`f_{sc}`$ is strangeness fraction corresponding to the critical density $`n_{Bc}`$ at the same conditions of charge and temperature. We see from Fig. 5 the critical strangeness fraction increases with increasing temperature. The stable area corresponds to the larger negative charge and higher strangeness fraction.

Finally, we hope to investigate the stability of SQM from another view point. Instead of considering the chemical equilibrium, we study the behavior of chemical potentials for electron and $`u`$ quark, respectively. The reason for this choice is that the electron ($`u`$ quark) is the particle with minimum mass but negative (positive) charge in our model. The smaller the particle mass is, the lower the chemical potential will be. When chemical potential tends to be negative, the system becomes unstable because the Fermi surface is negative in this case. Therefore, we can use the conditions $`\mu _e0`$, $`\mu _u0`$ to determine the stability of SQM.

The chemical potential of electron $`\mu _e`$ vs baryon number density $`n_B`$ at a fixed temperature $`T=100\text{MeV}`$ with various positive charge density $`Q=1\times 10^4,1\times 10^3,3\times 10^3,5\times 10^3`$ and $`7\times 10^3e\text{fm}^3`$ are drawn in figure 6. We see that $`\mu _e`$ becomes negative and the system becomes unstable when $`Q<0.001e\text{fm}^3`$. Similarly, the chemical potential of $`u`$ quark $`\mu _u`$ vs baryon number density $`n_B`$ at a fixed temperature $`T=100\text{MeV}`$ with various negative charge density $`Q=0.2,0.322,0.4,`$ $`0.6`$ and $`0.8`$ $`e\text{fm}^3`$ are drawn in figure 7. We see from this figure that when $`Q<`$ $`0.322`$, the curve will pass through a negative area and thus SQM is unstable. Therefore, at temperature $`T=100\text{MeV}`$, the charge density must be limited to $`0.001e\text{fm}^3>Q>0.322e\text{fm}^3`$ for stable SQM.

## IV Summary

In summary, by means of the QMDTD model, we find that the negative charge and the higher strangeness fraction are in favor of the stability of SQM at finite temperature. This result is in agreement with that given by Greiner et al. at zero temperature, but we extend their conclusion to finite temperature. A critical baryon number density $`n_{Bc}`$ under which the SQM is unstable has been obtained. We find that $`n_{Bc}`$ decreases rapidly when charge density $`Q`$ becomes more negative and approaches to zero at the deconfinement temperature $`T_c`$. By using the chemical potentials of electron and $`u`$ quark, we suggest a method to determine the limit of charge density for the stability of SQM. We hope the stable areas given by our model can have impact on the detectability of SQM in future heavy ion experiments.

## V Acknowledgment

This work was supported in part by the NNSF of China under contract Nos. 19975010, 10047005.

## VI Figure Captions

Fig.1 The quark fraction $`F_q`$ as a function of baryon number density $`n_B`$ with neutral charge, $`F_u,F_d`$ and $`F_s`$ are represented by solid, dashed and dotted lines, respectively. Here $`F_d>F_u>F_s.`$

Fig.2 The critical density $`n_B`$ as a function of temperature $`T`$ with neutral charge, the dashed line is for QMDD model and the solid line is for QMDTD model.

Fig.3 The critical density $`n_B`$ as a function of temperature $`T`$ with various charge density $`Q=0,0.1`$ and $`0.2e\text{fm}^3`$, they are represented by solid, dotted and dashed lines, respectively.

Fig.4 The quark fraction $`F_q`$ as a function of baryon number density $`n_B`$ with a negative charge density $`Q=0.2e\text{fm}^3`$, $`F_u,F_d`$ and $`F_s`$ are represented by solid, dashed and dotted lines, respectively. Now $`F_d>F_s>F_u.`$

Fig.5 The critical strangeness fraction $`f_{s\text{ }}`$as a function of temperature $`T`$. The solid line is for $`Q=0`$ and the dashed line is for $`Q=0.2e\text{fm}^3`$.

Fig.6 The chemical potential of electron $`\mu _e`$ as a function of baryon number density $`n_B`$ with the positive charge density $`Q=1\times 10^4,1\times 10^3,3\times 10^3,5\times 10^3`$ and $`7\times 10^3e\text{fm}^3`$. At a fixed temperature $`T=100\text{MeV}`$, they are represented by solid lines respectively.

Fig.7 The chemical potential of $`u`$ quark $`\mu _u`$ as a function of baryon number density $`n_B`$ with the negative charge density $`Q=0.2,0.322,0.4,0.6,`$and $`0.8e\text{fm}^3`$. At a fixed temperature $`T=100\text{MeV}`$, they are represented by solid lines respectively.