text:
# The 𝑈_𝐴(1) Problem on the Lattice Talk presented by G.C. Rossi

## 1 Introduction

In a nut-shell the $`U_A(1)`$ problem lies in the fact that, in the absence of $`U_A(1)`$ anomaly contributions, the unphysical bound $`m_{\mathrm{octet},I=0}<\sqrt{3}m_\pi `$ holds . In 1976 ’t Hooft  pointed out that the resolution of this problem had to be related to the existence of topologically non-trivial gauge field configurations in Euclidean QCD. Dwelling on this observation, it is possible to obtain an explicit formula for the $`\eta ^{}`$ mass either in the ’t Hooft limit ($`N_c\mathrm{}`$, with $`g^2N_c`$ and $`N_f`$ held fixed ), as was done in , or by assuming that anomalous flavour-singlet axial WTIs retain their validity order by order in an expansion in $`uN_f/N_c`$ around $`u=0`$, as argued in . In both cases one can derive the “leading-order” Witten–Veneziano (WV) relation

$$m_\eta ^{}^2=\frac{2N_f}{F_\pi ^2}A,A=d^4xQ(x)Q(0)|_{\mathrm{YM}},$$
(1)
where $`F_\pi 94\mathrm{MeV}`$ is the pion decay constant, $`A`$ the “topological susceptibility” and $`Q=\frac{g^2}{32\pi ^2}ϵ_{\mu \nu \rho \sigma }\mathrm{tr}[F_{\mu \nu }F_{\rho \sigma }]`$ the topological charge density. The notation $`\mathrm{}|_{\mathrm{YM}}`$ means that the $`QQ`$-correlator is to be computed in the pure Yang–Mills (YM) theory, i.e. in the absence of quarks.

The formal relations in eq. (1), when translated in any regularized version of QCD, such as lattice QCD, become more involved and quite a number of subtleties have to be dealt with in order to be able to determine the correct field theoretical expression of $`A`$, which should be used in the above equations. Two problems need to be solved to make the formulae in (1) precise and of practical use. One has to i) find a properly normalized lattice definition of the topological charge density, $`Q_L`$; ii) subtract from the $`Q_L(x)Q_L(0)`$ product appropriate contact terms, as required to make it an integrable (operator-valued) distribution.

A rigorous derivation of a formula for the $`\eta ^{}`$ mass, which could be unambiguously employed in numerical simulations, can be given , if use is made of the (lattice) regularized anomalous flavour-singlet axial WTIs of full QCD in order to construct a properly defined $`Q_LQ_L`$ correlation function. The remarkable result of this analysis is that no subtraction is needed to this end in the chiral limit, if the definition of lattice topological charge density suggested by fermions obeying the Ginsparg–Wilson (GW) relation  (e.g. overlap fermions ) is employed.

## 2 GW fermions and the $`\eta ^{}`$ mass formula

Regularizing the fermionic part of the QCD action using GW fermions  offers the great advantage over the standard Wilson discretization  that global chiral transformations can be defined, which are an exact symmetry of the massless theory, as in the formal continuum limit. This follows from the relation $`\gamma _5D+D\gamma _5=D\gamma _5D`$ , obeyed by the GW Dirac operator, $`D`$ 111We set the lattice spacing, $`a`$, equal to 1. In all lattice formulae the continuum limit is, however, understood.. In this regularization the $`U_A(1)`$ anomaly is recovered à la Fujikawa , because of the non-trivial Jacobian that accompanies the change of fermionic integration variables induced by a $`U_A(1)`$ lattice transformation. As a consequence the (anomalous) flavour-singlet WTIs in the presence of $`N_f`$ massless fermions take the form

$`_\mu A_\mu ^0(x)\widehat{O}={\displaystyle \frac{2N_f}{2}}\mathrm{Tr}[\gamma _5D(x,x)]\widehat{O}+\delta _A^x\widehat{O},`$
where $`\widehat{O}`$ is any finite (multi)local operator, $`\delta _A^x\widehat{O}`$ is its local variation and $`_\mu A_\mu ^0(x)`$ is the divergence of the singlet axial current.

Neither $`\mathrm{Tr}[\gamma _5D(x,x)]`$ nor $`A_\mu ^0(x)`$ are finite operators, but finite linear combinations, $`\widehat{Q}`$ and $`\widehat{A}_\mu ^0`$, can be easily constructed  by writing

$`\widehat{Q}(x)={\displaystyle \frac{1}{2}}\mathrm{Tr}[\gamma _5D(x,x)]{\displaystyle \frac{Z}{2N_f}}_\mu A_\mu ^0(x)`$ (2)
$`\widehat{A}_\mu ^0(x)=(1Z)A_\mu ^0(x).`$ (3)
where $`Z`$ is logarithmically divergent to lowest order in perturbation theory and vanishes as $`u0`$.

With the above definitions the renormalized singlet axial WTI becomes

$$_\mu \widehat{A}_\mu ^0(x)\widehat{O}=2N_f\widehat{Q}(x)\widehat{O}+\delta _A^x\widehat{O}.$$
(4)
This equation shows that $`\widehat{Q}`$ and $`\widehat{A}_\mu ^0`$ are correctly normalized and that there is no power-divergent mixing of $`\mathrm{Tr}[\gamma _5D(x,x)]`$ with the pseudo-scalar quark density, which would bring in a dangerous lineraly divergent mixing coefficient.

A formula for the $`\eta ^{}`$ mass can be obtained from the above WTIs, by observing that the $`\eta ^{}`$ mass must vanish as $`u0`$ for massless QCD to have no $`\theta `$-dependence. For that one starts by defining the lattice Green function 222A contact term, $`\mathrm{CT}(p)`$, should be added to the r.h.s. of eq. (5) to make it finite at $`p0`$. $`\mathrm{CT}(p)`$ is a polynomial of degree 4 in $`p`$, which vanishes at $`p=0`$, because the r.h.s. of eq. (5) is certainly finite (actually zero) at $`p=0`$. $`\mathrm{CT}(p)`$ plays no rôle in the argument below, since we will be finally only interested in the value of $`\chi _{tL}`$ at $`p=0`$. For brevity we will not indicate it explicitly in the following.

$$\chi _{tL}(p)=\frac{1}{2N_f}d^4x\mathrm{e}^{ipx}_\mu \widehat{A}_\mu ^0(x)\widehat{Q}(0),$$
(5)
with $`\widehat{Q}`$ given by eq. (2). In the full theory, where the $`\eta ^{}`$ is massive, one has $`\chi _{tL}(0)=0`$. On the other hand, since, as we obseved above, at the chiral point, the $`\eta ^{}`$ mass vanishes as $`u0`$, only the $`\eta ^{}`$ pole will contribute to $`\chi _{tL}(p)`$ in this limit, leading to the relation

$$\underset{p0}{lim}\underset{u0}{lim}\chi _{tL}(p)=\frac{F_\pi ^2}{2N_f}m_\eta ^{}^2|_{u=0}.$$
(6)
The l.h.s. of eq. (6) can be evaluated by using the WTI (4) with $`\widehat{O}=\widehat{Q}`$. Recalling that $`Z`$ vanishes when $`u0`$ and $`\delta _A^x\widehat{Q}=0`$, we obtain

$`{\displaystyle \frac{F_\pi ^2}{2N_f}}m_\eta ^{}^2|_{u=0}=\underset{p0}{lim}\underset{u0}{lim}{\displaystyle d^4x\mathrm{e}^{ipx}}`$
$`{\displaystyle \frac{1}{2}}\mathrm{Tr}[\gamma _5D(x,x)]{\displaystyle \frac{1}{2}}\mathrm{Tr}[\gamma _5D(0,0)].`$ (7)
The limits $`u0`$ and $`p0`$ in this expression can be readily performed in the order indicated, if one can assume that taking the first limit simply amounts to setting the fermion determinant equal to unity. One gets in this way

$`{\displaystyle \frac{F_\pi ^2}{2N_f}}m_\eta ^{}^2|_{u=0}=`$
$`={\displaystyle d^4x\frac{1}{2}\mathrm{Tr}[\gamma _5D(x,x)]\frac{1}{2}\mathrm{Tr}[\gamma _5D(0,0)]}|_{\mathrm{YM}}.`$ (8)
The restriction to pure YM theory, indicated in eq. (8), is an obvious consequence of the fact that, for a Green function of only gluonic operators, neglecting the fermion determinant is tantamount to limiting the functional integral to the pure gauge sector of the theory.

A relevant question at this point is to ask at which value of $`N_c`$ eq. (8) is supposed to be valid. The answer depends on the behaviour of QCD with $`N_f`$. Various scenarios are envisegeable.

1) In the most favorable situation, in which the limit $`u0`$ of $`\chi _{tL}(p)`$ exists and is equal to the value it takes at $`N_f=0`$ (fermion determinant equal to 1), formula (8) is valid for any value of $`N_c`$ and for each $`N_c`$ it yields the mass of the $`\eta ^{}`$ meson (at leading order in $`u=0`$) in the world with the corresponding number of colours .

2) If quenching can be attained only in the limit in which the number of colours goes to infinity, then eq. (8) will yield a formula for the $`m_\eta ^{}^2`$ valid up to O($`1/N_c`$) corrections .

3) Finally it may happen that taking the limit $`u0`$ does not correspond to quenching. In this case one cannot pass from the fairly complicated eq. (7) to the more useful formula (8).

Remembering that $`\frac{1}{2}\mathrm{Tr}[\gamma _5D(x,x)]`$ should be identified with the topological charge density , we conclude that eq. (8) can be rewritten in the very suggestive form

$$\frac{F_\pi ^2}{2N_f}m_\eta ^{}^2|_{u=0}=\underset{V\mathrm{}}{lim}\frac{(n_Rn_L)^2}{V},$$
(9)
where $`(n_Rn_L)^2`$ is the expectation value of the square of the index of the GW fermion operator, $`D`$, and $`V`$ is the physical volume of the lattice.

In the form (9) the $`\eta ^{}`$ mass formula can be directly compared with overlap fermion simulation data. Existing numerical results in 2 (abelian Schwinger model ) and in 4 (QCD ) dimensions agree quite well with theoretical expectations and/or experimental numbers .

## 3 Conclusions

We have shown that a formula for the $`\eta ^{}`$ mass can be rigorously derived in lattice QCD, exploiting the anomalous flavour-singlet axial WTIs of the theory. If fermions obeying the GW relation are used, there exists a natural definition of topological density which is correctly normalized and for which the naive form of the WV formula holds without the need of introducing any subtraction.

In the literature two sorts of approaches have been proposed to deal with the problem of computing $`A`$ on the lattice, which have led to numerical values as good as those obtained with the present method . The first one is based on a direct field theoretical definition of $`A`$  that takes into account the need for the renormalization of $`Q_L`$ and the subtraction of the operators $`F^2`$ and $`11`$ in the short distance expansion of $`Q_LQ_L`$. The second one makes use of the notion of “cooling”  to carry out the necessary operations of renormalization and subtraction. Both methods are reliable to the extent that they are able to capture the topological properties of the gauge field configurations that determine the number of zero modes of the Dirac operator. Simulations based on the geometrical definition of $`Q_L`$ of ref.  have not yet led to comparably good results . Acknowledgements \- I wish to thank the Organizers of Lattice2001 for the wonderful and exciting atmosphere of the Conference.


format: text

text:
# 1 Calorimeters geometry considered

## 1 Calorimeters geometry considered

The pion and electron responses are obtained with GEANT simulations with a detailed description of CMS calorimeter geometry (CMSIM 115), version TDR-2/1/.

The ECAL is the PbWO4 one (readout number 1).

An additional 1 cm scintillator layer is placed in front of HB to compensate for the energy loss in the cables, electronics and cooling system on the back side of the ECAL, which is simulated by a uniform slab of 0.2 $`\lambda `$ of dead material (readout number 2).

The sampling thicknesses of copper alloy (90%Cu + 10%Zn) are 5 cm in the barrel segments, except the inner and outer plates which are 7 cm stainless steel. The 9 mm gaps between the absorber plates are filed with 4 mm scintillator planes (readout number 3).

To improve the energy measurement for late developing hadron showers, tail catcher layers are inserted in front of and behind the first return yoke iron layer (RY1) (readout number 4).

The regions behind the ECAL are empty except for the scintillator layer and equivalent dead material uniform absorber layer. The effects of the CMS 4T field are fully included.

## 2 Data sample

We have performed a GEANT simulation of the response to pions for a set of incident energies from 5 to 3000 GeV and to electrons from 5 to 500 GeV at pseudorapidity set to 0. GHEISHA was used as a hadron interaction simulator. The sum of the energies deposited by the showers in the active elements of each readout has been stored in a disk file so that the showers could be analysed later.

## 3 Standard calibration

The standard calibration uses sampling coefficients for calorimeter parts with different sampling ratio. The reconstructed energy $`E^{rec}`$ of simple shower is given by the weighted sum of the energies deposited in the readouts:

$$E^{rec}=\underset{i=1,4}{}c_iE_i,$$

where:

$`E_i`$ \- amplitude of the signal from the calorimeter longitudinal segmentation (readouts);

$`c_i`$ \- calibration coefficients, determined by the minimisation of the width of the energy distributions.

The Gaussian part of the reconstructed energy distributions at various incident energies are then fitted to obtain the calorimeter energy resolution. The energy resolution is parametrised by the expression:

$$\sigma /E=a/\sqrt{E}b.$$

The energy resolution obtained by the standard calibration with simultaneously fit of pion and electron events is shown in fig. 1a) and 1b) and the residuals of the reconstructed energy in fig. 2a) and 2b) (open points).

## 4 Non-linear technique

The linear behaviour of the calorimeter response could be restored by application of the non-linear method which improve the linearity of the calorimeter response and the energy resolution in the broad energy range. Non-linear technique is the selection of some additional energy depended parameters which provide correct energy reconstruction of the showers.

Thus the reconstructed energy $`E_{nl}^{rec}`$ is parametrised as:

$$E_{nl}^{rec}=\underset{i=1,4}{}f_i(\stackrel{}{A},E_i)E_i,$$

where $`f_i`$ are non-linear functions of N unknown parameters $`A_j`$, j = 1, N and measured readout signals $`E_i`$.

The system is solved with autoregularized Newton type method/2/ by minimisation of the $`\chi ^2`$ functional (REGN code/3/)

$$\chi ^2=\underset{j=1,M}{}W_j(E_j^{in}E_{nl,j}^{rec})^2,$$

where M is number of generated set of energies multiplied by the number of generated events at each energy and $`W_j`$ are statistical weights.

In the REGN computer code the $`\chi ^2`$ is one of the different criteria available for solving of the system and for testing the mathematical model. The other criteria permit one to chose uniquely between several model functions the best one/4,5,6/.

This technique was applied for single pions calibration/7/. The obtained energy resolution after applying the non-linear calibration method for simultaneously fit of pion and electron beams is shown on fig. 1a) and 1b) (full points). Resolution for pions becomes better in comparison to standard linear calibration. A small degradation of electron energy resolution is observed. The comparison of standard calibration method to the non-linear technique gives a significant improvement in the linearity of the reconstructed energy for both pions and electrons (fig. 2a) and 2b)) (full points).

The unification of our technique with the most clear physical approaches8,9,10 can give better energy resolution because we can test strait different energy dependences of unknown parameters which depends of interaction cross-sections and calorimeters construction details.

## 5 Summary

The CMS calorimeters simultaneously pion and electron calibration was done using linear and non-linear approaches. The standard calibration method gives huge deviation from linearity of the reconstructed energy. The calibration using non-linear technique improves behaviour of the energy resolution for pions and restores the calorimeter linearity. Such non-linear technique can be tested successfully for pions, electrons and jets cases simultaneously.


text:
# Automated optical identification of a large complete northern hemisphere sample of flat spectrum radio sources with 𝑆_{6cm}>200 mJy

## 1 Introduction

Over the last decade, a comprehensive catalogue of compact flat spectrum radio sources, the Jodrell Bank $``$ VLA Astrometric Survey (JVAS), has been constructed (Patnaik et al, 1992; Browne et al., 1997; Wilkinson et al., 1998). The main aim of this survey was to provide the astronomical community with a network of bright radio sources with accurate positions ($`1015`$ mas) primarily intended for use as phase calibrators for the Jodrell Bank MERLIN, the VLA and VLBI networks. In addition to its astrometric goals, JVAS proved an effective means of finding radio sources exhibiting strong gravitational lensing, with six systems confirmed to date (King et al. 1999).

Its virtually complete northern-sky coverage makes JVAS uniquely suited for studying the high luminosity end of the flat spectrum radio source population. Our group has been particularly involved in the search for distant quasars (Hook et al. 1996, Hook et al. 1998, Hook & McMahon 1998). As part of this project, sources in the JVAS survey were optically identified using the catalogue output of the APM (the Automated Plate Measurement Facility at Cambridge) scans of the POSS-I plates. This catalogue is ideal for identifying large samples of objects, since they are generated in an automated way allowing an objective assessment of their reliability and homogeneity.

This paper describes the automated optical identification procedure of JVAS, as used for the search for distant quasars. It can be used as a guidance for future projects, involving the optical identification of comprehensive samples of flat spectrum radio sources, like those in the Cosmic Lens All Sky Survey (CLASS; Myers et al. 2001). Since JVAS is not complete, a new sample is constructed using selection criteria similar to those of JVAS, but with the emphasis on completeness, and with the use of more accurate selection surveys. By comparing this sample with JVAS, the completeness of JVAS is assessed. This is described in section 2. Section 3 and 4 describe the APM-POSS-I catalogue and the optical identification procedure of JVAS and its complete counterpart. Section 5 gives the results and discusses possible relations between optical and radio properties using the complete sample.

## 2 A complete sample of flat spectrum radio sources

### 2.1 The Jodrell Bank $``$ VLA Astrometric Survey (JVAS)

The JVAS catalogue has been presented in three separate papers, for the regions $`+35^{}\delta +75^{}`$ (Patnaik et al., 1992), $`+0^{}\delta +20^{}`$ (Browne et al., 1997), and $`+20^{}\delta +35^{}`$ and $`+75^{}\delta +90^{}`$ (Wilkinson et al., 1998). Their initial source list was constructed using the Greenbank surveys conducted at 1.4 and 5 GHz by Condon & Broderick (1985, 1986) and Condon, Broderick & Seielstad (1989), selecting all sources with spectral indices larger than $`\alpha =0.5`$ ($`\alpha `$ defined as $`S_\nu \nu ^\alpha `$) and $`S_{5GHz}200`$ mJy, excluding the region $`|b|<2^{}.5`$. The 5 GHz flux densities used to construct the sample were determined directly from the maps and have been found to be systematically higher by $`10\%`$ than the flux densities determined by Gregory and Condon (87GB, 1991) from the same maps. Since the main goal of JVAS was to construct a grid of phase calibrator sources, and the distribution of the sample above contained some undesirable large holes ($`5^{}`$) in some areas, a few additional sources were selected with $`S_{5GHz}150`$ mJy to fill these holes. Where possible, the potential calibrator sources were cross-checked against lower frequency surveys to make sure they had genuine flat spectra. The resulting source list was observed with the VLA at 8.4 GHz between 1990 and 1992 (Patnaik et al, 1992; Browne et al., 1997; Wilkinson et al., 1998). The resulting catalogue includes sources with observed 8.4 GHz peak brightness $`50`$ mJy/beam.

For sources that are in the regions $`+0^{}\delta +20^{}`$ and $`+35^{}\delta +75^{}`$ the rms position error was estimated to be approximately 10 mas in each coordinate (Patnaik et al., 1992; Browne et al., 1997). For sources in the regions $`+20^{}\delta +35^{}`$ and $`+75^{}\delta +90^{}`$ this was estimated to approximately be 40 mas (Wilkinson et al., 1998). The resulting catalogue contains 2121 radio sources.

### 2.2 The selection of a complete flat spectrum sample

The JVAS catalogue was primarily intended to provide phase calibrators with a uniform sky distribution, and was not intended as a statistically complete flat spectrum sample. Therefore, a similar, complete sample has to be defined to be able to conduct statistically meaningful optical/radio studies. This sample, which we call JVAS++, will then also be used to assess the completeness of JVAS.

The JVAS++ sample was constructed using similar selection criteria as used for the Cosmic Lens All Sky Survey (CLASS; Myers et al. 2001), and is based on the recently available and more accurate radio surveys, the GB6 at 5 GHz (Gregory et al. 1996), and the NVSS at 1.4 GHz (Condon et al. 1998). All the relevant parameters of the radio surveys involved are given in table 1. In the selection procedure, the 1.4 GHz flux density of a GB6 radio source is defined as all the NVSS flux within 70′′ of the GB6 position (as used for CLASS). The sample has the following selection criteria:

* GB6 5 GHz flux density, $`S_{5GHz}>200`$ mJy
* declination, $`0^{}<\delta <75^{}`$, galactic latitude, $`|b|>30^{}`$.
* radio spectral index, $`\alpha _{1.45GHz}>0.5`$.

In this way, 915 sources were selected. The distribution on the sky is shown in figure 1. The high galactic latitude cut-off was chosen to reduce the possible influence of galactic foreground extinction. Furthermore, the APM catalogue is only available at these latitudes. About 16% (153) sources are not in the JVAS sample and have been observed with the VLA in B configuration at 8.4 GHz as part of the CLASS survey, in a similar manner as the targets in the JVAS sample (Myers et al. 2001). These provided us with radio positions at an accuracy similar to those of the JVAS sources. Fifty-eight sources are fitted by multiple components at 8.4 GHz. For these, the positions of the brightest components are used. Fifteen objects are not detected with the VLA, all exhibiting extended structure on arcminute scales.

### 2.3 The completeness and reliability of JVAS

JVAS contains 762 sources which overlap with JVAS++. In the same area of sky, another 371 sources are part of JVAS, but are not in JVAS++. In addition, 153 objects are in JVAS++, but are not part of JVAS.

The large differences between the JVAS and JVAS++ are due to several effects:

* the use of different selection surveys for JVAS.
* the use of a difference flux density scale for JVAS, which was found to be offset by 10% at 5 GHz.
* the use of a variety of spectral data in the target selection, the non-detection of some sources, and the exclusion of some sources for which the observed VLA position was found to be offsetted by more than 40 arcseconds from the pointing position.

The use of a different flux density scale at 5 GHz, implies that the actual selection criteria for JVAS are, $`S_{5GHz}>180`$ mJy, and $`\alpha _{1.45}>0.58`$. This effect, in combination with the inclusion of sources with 150 mJy $`<S_{5GHz}<`$ 180 mJy, to fill holes in the sky distribution, causes the low reliability of JVAS. Note that since these sources are fainter than 200 mJy, they do not influence the reliability above this flux density, which is typically $`>95\%`$ (see fig. 2, left). However, they are included, when the reliability is plotted as function of spectral index (fig. 2, right), causing it to be typically $`<75\%`$.

The use of different selection surveys has three effects. Firstly, the Greenbank 1400 MHz survey, used for the selection of JVAS, has a 12 FWHM resolution, compared to 45′′ for the NVSS and a search radius of 70′′ used for the revised ‘complete’ sample. Therefore, sources with extended structure at arcminute scales may have been excluded from JVAS with $`\alpha <0.5`$, but included in the comparison sample. This is reflected in Figure 2 (left), showing that the completeness of JVAS is decreasing to $`50\%`$ at the steep spectral index cut-off. Indeed 49 objects exhibit extended structure at arcminute scales. These are shown in the appendix, where their NVSS images are overlaid with data from the Digitized Sky Survey (Lasker et al. 1990).

Secondly, the use of different selection surveys may cause sources to be included or excluded due to their possible intrinsic variability. Thirdly, measurement uncertainties may cause differences in the selection. To show the influences of variability and measurement uncertainties on the selection process, the Greenbank 1400 MHz observations were compared with the NVSS, and the 87GB data were compared with the GB6. Figure 3 shows the normalised differences between the GB1400 and NVSS flux densities for the sources in the complete sample, with no radio spectral selection on the left, ‘steep’ spectrum ($`\alpha <0.2`$) sources in the middle, and ‘flat’ spectrum ($`\alpha >0.2`$) sources in the right panel. About 12% of the sources do not have a GB1400 flux density since they are too faint to be in the catalog by White & Becker (1992). These are excluded from the analysis. The NVSS flux densities were decreased by 3% to match the GB1400 flux density scale. The distributions are overlaid by Gaussians with $`\sigma =`$ 0.13, 0.13, and 0.25 respectively. The distribution for the flat spectrum sources is clearly much broader ($`\sigma =0.25`$) than that for steep spectrum objects ($`\sigma `$ 0.13). This indicates that the differences between GB1400 and the NVSS for flat spectrum sources is dominated by variability, which is about 2 times larger than the measurement uncertainties. From a similar distribution of normalised flux density differences between NVSS and the FIRST survey (White et al. 1997), we found that the measurements errors in both FIRST and NVSS are about 5%. This means that, assuming that the distribution of the steep spectrum sources is solely due to measurement errors, the GB1400 survey has an uncertainty in flux density of 12%. Note that this uncertainty has a proportional term and a constant term, the latter due to noise and confusion. In FIRST and NVSS, this constant term is unimportant for our sample. However for the GB1400 survey with a noise level of $`25`$ mJy, this factor probably dominates the uncertainty in flux density for the faintest sources.

A similar analysis has been carried out for the two selection surveys at 5 GHz, which is shown in figure 4. The GB6 and 87GB are not independent; GB6 is based half of on the 87GB dataset, and half on a similar dataset, taken a year earlier. Since we are interested in the variability aspect, we assumed that the flux densities of the sources in the second dataset is $`2\times GB687GB`$. The normalised differences between the flux densities taken at these epochs are shown in figure 4, with on the left with no spectral selection, in the middle for steep spectrum sources, and at the right for flat spectrum sources. The distributions are overlaid with a Gaussian with a $`\sigma `$ of 0.14, 0.14, and 0.17 respectively. Assuming that the distribution of the steep spectrum sources is solely due to measurement uncertainties, each dataset has a typical error of 10%, and the uncertainty in flux density of the GB6 is 7%. The constant noise term of GB6 is expected to contribute for only 2% for the faintest sources in the sample. Note that the flat spectrum sources seem to be only slightly variable, compared to the data at 1.4 GHz, although flat spectrum sources are known to be increasingly variable towards higher frequencies. This is due to the fact that the time-baseline at 5 GHz ($``$ 1 year) is much shorter than at 1.4 GHz ($``$ 10 years).

It is difficult to estimate the completeness and reliability of JVAS and JVAS++, due to all the effects described above. Most of the differences between the two samples are not caused by incompleteness or unreliability of JVAS, but due to a difference in selection criteria. To assess the completeness of JVAS and JVAS++, we performed a simple simulation. We treated the NVSS and GB6 data as infinitely accurate and added Gaussian distributed offsets (with $`\sigma `$’s as estimated above) to the flux densities, and determined how many sources enter and leave the sample. In this way, the completeness and reliability of JVAS++ were estimated both to be 95% respectively. For JVAS, we changed the selection criteria to $`\alpha _{1.45GHz}>0.58`$ and $`S_{5GHz}>180`$ mJy after the offsets were added, to mimic the use of a different flux density scale. In addition it was assumed that all the sources with $`S<160`$ mJy were added artificially in JVAS to fill in holes in the sky distribution ($``$100 objects). In this way the completeness and reliability of JVAS was estimated to be 90% and 70% respectively.111Note that all the sources in JVAS are real and that their positions are reliable

## 3 The automated optical identification procedure

### 3.1 The APM $``$ POSS-I catalogue

Optical identification of sources in the JVAS and JVAS++ samples was carried out using the output catalogue of the APM scans of the Palomar Sky Survey (POSS-I) photographic plates in the $`e`$ (red) and $`o`$ (blue) passbands (Bunclark & Irwin, 1983). Plates were scanned with $`0.5^{\prime \prime }`$ pixels. The APM covers the Northern sky with $`|b|20^{}30^{}`$. For each detected object a magnitude was measured and the image classified for the $`e`$ and $`o`$ plates separately. The images were classified as galaxies, stellar, merged objects or noise based on combinations of various image parameters such as moments, ellipticities etc. The spread in the distribution of each parameter with respect to the stellar locus was used to weight that parameter when combined to form a final classification parameter. The distribution of the combined parameter has a mean and spread which is dependent on magnitude, and a statistical correction was calculated to convert these for the stellar objects to Gaussian distributions with zero mean and rms=1. This correction was then applied to all images and the standard deviation from a stellar profile $`N\sigma _c`$ was recorded. For images about 2 magnitude above the plate limit, this classification is $`90\%`$ accurate.

The $`e`$ and $`o`$ catalogues were merged, assuming that sources within 2′′ are one and the same object. Objects with centroids on the $`e`$ and $`o`$ plates which differ in position between 2 and 4′′ were considered to be “nearly” matches. In such case only the position on the E plate were preserved. The final positions of the objects were corrected for subtle distortions of the plates which are most prominent near the plate edges.

The zero-point for photometry in each field was calculated by assuming that the $`e`$ plate has a limiting magnitude of 20.0 and by assuming a universal, magnitude independent position in the colour-magnitude plane for the stellar locus, as described by McMahon & Irwin 1992. The zero-point of the magnitude system has an rms accuracy of 0.3 magnitude.

The APM $``$ POSS-I catalogue can be accessed via the internet at www.ast.cam.ac.uk/$``$apmcat.

### 3.2 The optical-radio correlation

Three source lists were constructed:

* Sources which are part of both JVAS++ and the JVAS sample (762 objects)
* Sources which are only part of JVAS++ (153 objects)
* Sources which are only part of JVAS (371 objects)

Hence, the combination of sources in lists I and II make up JVAS++, and the sources in lists I and III make up the JVAS sample (in the selected area of sky). The three lists were correlated with the APM$``$POSS-I catalogue, selecting all optical objects within $`30^{\prime \prime }`$ of each 8.4 GHz radio position. For the objects with extended radio structure at arcminute scales, the radio-optical overlays, as shown in the appendix, were used to search for a possible bright optical identification. These were added manually. The distributions of $`\mathrm{\Delta }\alpha _{optrad}`$ and $`\mathrm{\Delta }\delta _{optrad}`$ are shown in figure 5, for sources in the JVAS++ sample, excluding objects which exhibit multiple components at 8.4 GHz and/or extended structure at arcminute scales. From these plots the accuracy of the APM positions was determined, by fitting the $`\mathrm{\Delta }\alpha `$ and $`\mathrm{\Delta }\delta `$ distributions with a Gaussian on top of a flat distribution, as expected for a combined population of genuine identifications and random background objects. It shows that the APM positions have an uncertainty of rms=0.4′′ in both $`\alpha `$ and $`\delta `$, and show no systematic offsets to the radio positions.

Figure 6 shows the distribution of sources around the JVAS++ positions. The solid line is a fit to the background source density of 0.0006 arcsec-2. We found that this background level was not enhanced by possibly related cluster objects around the radio sources, by counting the number of objects within a 30′′ radius, offset 5 north from each source. It is therefore expected that $`2\%`$ of the sources have a random background object within 3 arcsec from their radio position. Such a background object will be incorrectly identified as the optical identification unless the genuine identification is detected and has an APM position closer to the radio position.

Unfortunately, for several reasons the errors in the optical positions are not Gaussian distributed, as assumed above:

* Bright galaxy identifications: The central positions of bright extended objects, e.g. nearby galaxies, can be determined significantly less accurate in the APM than those of fainter stellar objects.
* Merged Objects: Due to the limited resolution of the APM scans, two or more nearby sources are sometimes merged into one single object. In such a case, the central position of the merged object can be several arcseconds away from the position of one of the individual objects (i.e. the identification).
* Bright Stars: Bright Stars, can obliterate other objects in their vicinity, or can produce spurious images.
* Anomalies in the APM Catalogue:. In very rare cases the APM can produce nonsensical results unrelated to any of the former effects, e.g. plate defects such as dirt and scratches, aeroplane tracks and asteroids.

These effects produce a broad non-Gaussian tail in the error distribution at a $`510\%`$ level. To identify these problems for the individual objects, all the APM images were checked by eye and compared with images from the Digital Sky Survey as retrieved from Skyview. In the course of this process the images were also searched for possible surrounding groups or clusters of galaxies. In addition, if no identification in the APM was found, it was checked whether a possibly faint identification was visible in the DSS. All the JVAS and JVAS++ sources which were found to have one of these problems are shown in the appendix.

### 3.3 Completeness and reliability of the optical identifications

The likelihood method (de Ruiter, Arp and Willis, 1977) was used to quantify the identification procedure. First, for each candidate identification the following dimensionless measure of the uncertainty in position difference, $`r`$, was calculated:

$$r=\sqrt{\frac{\mathrm{\Delta }\alpha ^2}{\sigma _{\alpha opt}^2}+\frac{\mathrm{\Delta }\delta ^2}{\sigma _{\delta opt}^2}}$$
(1)
where $`\mathrm{\Delta }\alpha `$ and $`\mathrm{\Delta }\delta `$ are the offsets between the optical and radio positions, and $`\sigma _{\alpha opt}`$ and $`\sigma _{\delta opt}`$ the uncertainties in the optical right ascension and declination positions. The error in the radio position is small $`<0.05^{\prime \prime }`$ compared to the optical error ($`0.55^{\prime \prime }`$) and therefore neglected. Given the normalised position difference $`r`$ for a certain radio-optical pair, the likelihood ratio is defined by (de Ruiter, Arp and Willis, 1977)

$$LR(r)=\frac{dp(r|id)}{dp(r|c)}=\frac{re^{\frac{r^2}{2}}}{2\lambda re^{\lambda r^2}}=\frac{1}{2\lambda }e^{\frac{r^2(2\lambda 1)}{2}}$$
(2)
where $`\lambda =\pi \sigma _{\alpha opt}\sigma _{\delta opt}\rho _{bg}=0.5\rho _{bg}`$, where $`\rho _{bg}`$ is the background source density which is stored for each POSS-plate in APM calibration files. This is the ratio of the probability that a given object found between $`r`$ and $`r+dr`$ is the correct identification $`p(r|id)`$, divided by the probability that it is a contaminating object $`p(r|c)`$. The cumulative distribution of likelihood ratios is shown in figure 7.

The likelihood ratio cutoff used for this sample is 1.0 which means the probability that the given object is the correct identification is at least equal to the probability that it is a background object. Using this cutoff we find an identification rate, $`\theta `$, of 83% (i.e. 756 out of 915). This identification fraction can be related to the two a posteriori probabilities that the object found at an angular distance $`r`$ from the radio source position is a genuine identification, $`p(id|r)`$, or a confusing object, $`p(c|r)`$, using Bayes’ theorem

$$p(id|r)=\frac{\theta LR(r)}{\theta LR(r)+1\theta }p(c|r)=\frac{1}{\theta LR(r)+1\theta }$$
(3)
It was assumed that the sources which were found to be merged objects or bright galaxies have an infinite likelihood ratio, resulting in $`p(id|r)=1`$ and $`p(c|r)=0`$. The likelihood ratios of the objects, for which the manual check has shown they have a genuine optical identification but a large radio-optical position offset, are set to $`2\times 10^3`$, as if the offset was zero, and therefore as 100% reliable identifications. The completeness, $`C`$, of the identifications, which is the number of accepted identifications over the number of correct id’s, and the reliability, $`R`$, which is the fraction of identified sources for which the identification is correct, are given by:

$$C=1\left(\underset{L<1}{}p(id|r)\right)/N_{id}$$
(4)
$$R=1\left(\underset{L>1}{}p(c|r)\right)/N_{id}$$
(5)
where, $`N_{id}`$ is the total number of identifications. These correspond to a completeness $`C=99\%`$ and a reliability $`R=99\%`$ for LR$`1`$ used.

## 4 Results and discussion

Three tables were produced for the three source lists, which are shown in the appendix, with a detailed description of the columns.

### 4.1 Radio spectral properties and optical identifications

The JVAS++ sample is uniquely suited to study the high luminosity end of the flat spectrum radio source population. Figure 8 shows the $`e`$-magnitudes of the optical identifications as function of their $`oe`$ colours, classified as stellar (filled circles) and extended (open circles) objects on the red plates. Objects with only limits on the colours are not shown. The colours of the extended objects are more biased towards the red than those of the stellar identifications, as expected for quasars and galaxies, which are found to have generally blue and red optical colours respectively, This can also be seen in figure 9 where the colour-distribution of the two classes of objects are shown. Only at $`e>19`$, near the plate limit, does this picture get blurred, since no reliable classification can be made and the large majority of objects are classified as stellar.

Figure 10 shows the optical identification fraction as function of GB6 flux density. The upper line indicates the total identification fraction (including extended, stellar and merged objects). The middle line indicates the fraction of JVAS++ sources identified with APM objects classified as stellar in the red band. The lower line indicates the fraction of JVAS sources identified with APM objects classified as extended in the red band. There is a hint that both the total and stellar identification fraction decrease with decreasing flux density. These two effects are likely caused by a change-over to galaxy identifications at faint flux density levels, which leads to more loss against the plate limit. A similar effect at a similar flux density level has been seen by Shaver et al. (1997) and Falco, Kochanek & Muñoz (1998). Note that due to the uncertain classification near the plate limit, the true decrease in identification fraction for quasars with flux density may even be more pronounced.

Figure 11 shows the optical identification fraction as function of 5.0 to 1.4 GHz spectral index. The lines indicate the specific identification fractions as in figure 10. Clearly, the total and stellar identification fractions decrease towards steeper spectral indices, while the fraction of sources identified with extended objects appears to be independent of spectral index. This trend can be explained by assuming that at $`\alpha <0.3`$, the unbeamed population of ‘steep’ spectrum galaxies contributes significantly to the total radio source population. Most of these radio galaxies will be too faint to appear on the POSS, and will be unidentified radio sources. This results in a drop in the total identification fraction and the identification fraction of quasars.

For the optical identifications with available redshifts in the literature, the $`oe`$ colours are shown as function of redshift (figure 12). This figure should be interpreted with great care, since it includes a strong observational bias, e.g. towards bright galaxies and red quasars. Therefore, not surprisingly, most of the extended optical objects (galaxies) are observed to be at $`z<0.5`$. They are clearly redder than the stellar objects at similar redshifts. Note however that a few objects which are classified as extended, are found at much larger redshifts. Their optical spectra show that they are actually quasars, wrongly classified as extended objects in the APM. There is a clear trend that the optical colours of the quasars become redder towards high redshifts. This is a well known effect and the forms the basis of the colour selection of candidate high redshift quasars (Hook et al. 1996), and is due to intervening Ly$`\alpha `$ absorption systems. Note however, that the majority of the objects in this redshift regime was actually selected for spectroscopic follow up on the basis of their red colour ($`oe1.0`$), which may have strengthen this effect. The solid line indicates the expected $`oe`$ colour of an average quasar reddened by intergalactic absorption. As a quasar spectrum, a power law with spectral index of $``$0.5 was used with emission lines taken from the composite quasar spectrum as constructed by Francis et al. (1991). For the intergalactic absorption the model of Madau (1995) was used. The model follows the trend of redder colours towards high redshift well. At low redshift the data are systematically redder than the model. This is probably due to a contribution of underlying galaxy light.

It is interesting to investigate the relation between the optical apparent magnitude of a quasar and its radio flux density, since both should be related to the power output of the central engine. Bolton & Wall (1970) showed that objects at fainter flux densities tend to have fainter optical magnitudes, but that there is a large spread. This trend is also present in figure 13 (left), where the e-magnitude is shown as function of 5 GHz flux density for sources which are classified as stellar on the red plates, and which have a blue $`oe<1.0`$ colour. This can also be seen in figure 14, where for the same subsample, the magnitude distributions are shown for bright (GB6$`>`$1 Jy), intermediate (350 mJy $`>`$ GB6 $`>`$1 Jy), and faint (GB6 $`<`$ 350 mJy) radio sources.

It is interesting to investigate what may cause this large spread in the ratio of optical to radio luminosity. We noticed that many of the optically faintest quasars were known to be High Polarised Quasars (HPQ); objects which exhibit optical polarisation at levels $`>3\%`$. In the right panel of figure 13 the $`e`$ magnitude versus 5 GHz flux density is shown for objects in the sample which have their optical polarisation properties studied by Fugmann & Meisenheimer (1988) or Wills et al. (1992). Indeed, it shows that the HPQ (’H’) prefer to have high radio to optical flux ratios while the low polarised quasars (LPQ; ’L’) show low ratios. Objects at $`z<0.5`$ were excluded from this plot, to avoid possible host galaxy contamination, and to exclude the population of nearby BL Lac objects and optically violent variable (OVV) quasars which also show high optical polarisation, but often combined with low or rapidly varying radio to optical flux density ratios.

The distribution of HPQs and LPQs is consistent with the idea that the optical quasar light is a combination of a polarised synchrotron component plus an unpolarised component (e.g. Smith et al. 1994). If the unpolarised component is significantly brighter than the synchrotron component (LPQ), then the radio flux (also synchrotron) will be relatively faint compared to the optical, resulting in a ‘flat’ optical to radio spectral index. If the unpolarised component is fainter (HPQ), and a constant spectral index is assumed for the synchrotron component, then the radio flux will be relatively bright, resulting in an overall ‘steep’ radio-to-optical spectral index. Indeed, the dashed line in figure 13 (right) indicates an optical-to-radio spectral index of $``$0.75, with most of the LPQs and HPQs situated above and below this line respectively.

The most straightforward explanation of why the ratio of the polarised to unpolarised components varies so much from quasar to quasar is Doppler boosting of the synchrotron component (eg. Wills et al. 1992). We therefore hypothesise that the large spread in optical to radio luminosity ratios is caused by source to source variations of Doppler boosting of the radio flux, leaving the unpolarised component of the optical emission unaffected. The fact that in our sample we see the HPQs fainter than the LPQs is counter-intuitive since these are the boosted objects. However, in this scheme, the HPQs are not the boosted counterparts of the LPQs in the sample, but are the boosted counterparts of a population which is fainter in the optical and radio. If this explanation is correct, then a correlation is expected between the Doppler factor and the optical to radio spectral index. Lähteenmäki & Valtaoja (1999) estimated the Doppler factors of eighty of the brightest flat spectrum objects in the sky, using total flux density variation monitoring data at 22 and 37 GHz. Thirty-nine objects in their sample overlap with JVAS. Nineteen of those are located at $`z>0.5`$, and have APM identifications with $`oe<1.0`$ (to avoid possible influence of the host galaxy and extinction). For the optical and radio flux densities, the average of $`o`$ and $`e`$ magnitudes, and the average of 1.4, 5.0 and 8.4 GHz flux densities are used, to minimise the influence of variability. The estimated Doppler boosting, $`D_{var}^3`$, is plotted against the optical-to-radio spectral index in figure 15. It shows that the Doppler factor is indeed correlated (97% significance) with the radio-to-optical spectral index, and that the slope of the relation is as expected as for the hypothesis that the large spread in optical to radio luminosity ratios is caused by source to source variations of Doppler boosting of the radio flux leaving most of the optical emission unaffected. Evidently, this scheme is too simplistic. First of all, the spectral index of the synchrotron component is most likely steeper in the optical than in the radio, causing the boosting to be stronger in the optical than in the radio, changing the observed optical-to-radio spectral index of the synchrotron component. Furthermore, in the optical we may see a related, but younger part of the synchrotron component, which possibly exhibits different boosting properties, complicating the simple picture scetched above.

However, the main idea is supported by two previous studies. Firstly, Yee & Oke (1978) and Shuder (1981) showed that the emission line luminosity for a range of AGN types is proportional to the luminosity of the underlying optical continuum over four orders of magnitude. Secondly, Rawlings & Saunders (1991) found that the emission line luminosity of an unbiased sample of FRII radio galaxies is approximately proportional to the total jet kinetic power, which is closely coupled to the power of the central engine. This implies that the optical luminosity should also be a direct indicator of the jet power, hardly affected by Doppler boosting. Indeed, Wills & Brotherton (1995) use the ratio of the radio core to optical continuum luminosity, $`R_V`$, which is the equivalent to the radio-optical spectral index as used in this paper, as an improved measure of quasar orientation over the ratio of radio-core to lobe flux density, $`R`$. They show that the use of $`R_V`$, rather than $`R`$, results in a significantly improved inverse correlations with the beaming angle as deduced from apparent superluminal velocities and inverse-Compton-scattered X-ray emission, and with the FWHM of a quasar’s broad H$`\beta `$ emission line. The use of optical continuum luminosity rather than extended radio luminosity to represent the unbeamed jet power probably works better because the latter is more affected by source-to-source differences in the intergalactic medium (Wills & Brotherton, 1995), and source age.

## 5 Summary

We have described the automated optical identification procedure of the sources from the Jodrell Bank $``$ VLA Astrometric Survey (JVAS), and a similar, complete radio sample, JVAS++, using the APM scans of the POSS-I plates. It yields an identification rate of 83%, with a completeness and reliability of both 99%. About 20% is identified with extended objects, eg. galaxies. The identification rate appears to drop towards lower flux densities, and towards steeper radio spectra, especially for the stellar classifications. Furthermore, the optical fluxes of quasars with faint radio flux densities appear to be biased towards fainter magnitudes, although there is a large spread in the optical-to-radio spectral index. It is shown that this large spread in radio-to-optical spectral index may be caused by source to source variations in the Doppler boosting of the synchrotron emission.

## Acknowledgements

We thank Robert Sharp for calculating the optical colours of a template quasar as function of redshift. We thank the referee, Jasper Wall, for valuable comments. We acknowledge the use of NASA’s SkyView facility (http://skyview.gsfc.nasa.gov) located at NASA Goddard Space Flight Center. This research was in part funded by the European Commission under contract ERBFMRX-CT96-0034 (CERES).

## Appendix A

Tables 2,3, 4 show the optical-radio catalogue for the source lists I, II, and III, as defined in section 2. For each table, column 1 gives the name (JVS=JVAS, CLS=CLASS), columns 2 the radio position (J2000), columns 3, 4, and 5 give the optical-radio offset in right ascension, declination and in total. Column 6 gives the logarithm of the likelihood ratio. Columns 7 to 12 give the APM magnitude, classification (-1=stellar, 1=extended, 2=merged) and psf for the red and blue plate respectively. Column 13 gives the APM colour, columns 14, to 17 give the radio flux densities at 5 GHz (GB6), 1.4 GHz (NVSS and GB1400), and at 8.4 GHz (JVAS). Column 18 gives the redshift as found in the NASA/IPAC Extragalactic Database (NED). Column 19 gives a possible comment. The comments are explained in table 4. The optical parameters are not shown if $`\mathrm{\Delta }\alpha ,\mathrm{\Delta }\delta >3.0^{\prime \prime }`$, unless the check by eye has shown that a larger optical-radio position offset is the result of a bright galaxy or a merged object identification, or due to extended radio emission.

Figure 17 show all the sources for which the APM did not give a good representation. The three panels show $`4^{}\times 4^{}`$ representations of the Digitized Sky Survey, The APM-red, and the APM blue data. The number in the left corner indicates the optical-radio position offset.

Figure 16 shows contour plots of NVSS data of objects in the complete sample showing extended structure on arcminute scale. The greyscales represent optical Digitized Sky Survey data. Image sizes are 12’. The dashed circle indicate a 70” radius around the GB6 position.


text:
# Metastable states in the Blume-Emery-Griffiths spin glass model

## Abstract

We study the Blume-Emery-Griffiths spin glass model in presence of an attractive coupling between real replicas, and evaluate the effective potential as a function of the density overlap. We find that there is a region, above the first order transition of the model, where metastable states with a large density overlap exist. The line where these metastable states appear should correspond to a purely dynamical transition, with a breaking of ergodicity. Differently from what happens in $`p`$-spin glasses, in this model the dynamical transition would not be the precursor of a 1-step RSB transition, but (probably) of a full RSB transition.

The Blume-Emery-Griffiths (BEG) spin glass model  can be seen as a generalization of the Sherrington Kirkpatrick (SK) model, in which each site carries, beside to the spin variable $`S_i=\pm 1`$, a lattice gas variable $`n_i=0,1`$. Its Hamiltonian is given by

$$=\underset{i<j}{}J_{ij}S_iS_jn_in_j\frac{K}{N}\underset{i<j}{}n_in_j\mu \underset{i}{}n_i$$
(1)
where $`J_{ij}`$ are quenched Gaussian variables with zero mean and variance $`J^2/N`$. The case $`K/J=0`$ goes also under the name of Ghatak Sherrington model . It has been found that in general, fixed the value $`K/J`$, there exists a tricritical point in the phase space. For higher temperature or chemical potential the model presents a second order transition from a paramagnetic to a spin glass phase. For lower temperature or chemical potential the transition becomes first order, with a jump in the mean density $`d=n_i`$ and overlap $`q_{\alpha \beta }=S_i^\alpha n_i^\alpha S_i^\beta n_i^\beta `$ between two replicas $`\alpha `$ and $`\beta `$ (see Fig. 1 for the case $`K/J=0`$). Although the overlap parameter can become discontinuous, the transition seems to be quite different from the transition of other discontinuous spin glasses, like the $`p`$-spin models. A first obvious difference is that the transition of the BEG spin glass is first order also in the Ehrenfest sense, that is with a jump in the mean energy and entropy. Moreover, while in $`p`$-spin models the first step of replica symmetry breaking (RSB) turns out to be exact, in the BEG spin glass the question of the order of RSB appears much more subtle. The question was studied by many authors in the case of the Ghatak Sherrington model . It was found that the replica symmetric solution is unstable against RSB in the whole spin glass region, with the appearance of complex eigenvalues in the stability matrix. Recently , it has been shown that this kind of instability persists also in the 1-step RSB solution. Therefore it seems reasonable that only a full replica symmetry breaking may probably give the correct solution, as it happens in the SK model.

In the last years, a version of the model called frustrated Ising lattice gas (FILG), has been extensively studied as a model of a structural glass . In this version of the model the disordered interactions are of the type $`J_{ij}=\pm J`$, and $`K=J`$. Indeed, the FILG can be interpreted as a model for a system of asymmetrical molecules, with the variables $`n_i`$ representing the presence of a molecule on the site $`i`$, and the variables $`S_i`$ representing the spatial orientation of the molecule. Two neighbor particles feel a repulsion of intensity $`2J`$ if their relative orientation is “wrong”, that is $`J_{ij}S_iS_j<0`$. Numerical simulations on finite dimensional lattices show that at low temperature the model develops a two step relaxation , as observed in supercooled liquids and predicted by the mode coupling theory (MCT) of the glass transition .

Given these links between the model and the glass transition in structural liquids, it would be interesting to study (at least in mean field) if the discontinuous transition of the model is preceded at a higher temperature by a purely dynamical transition, from an ergodic to a non-ergodic phase, as it happens for $`p`$-spin glasses . This would require the solution of the dynamical equations, as done for the $`p`$-spin, to see if they show a singularity, and of what kind, at some temperature greater than the static one. In this paper we address the problem from another point of view, exploiting the effective potential theory, introduced some time ago  as a tool to identify the presence of metastable states, and a transition from an ergodic to a non-ergodic phase. We will refer to the so-called “annealed version” of the method. One considers two “real” replicas of the system (as opposed to the replicas generated by the replica trick), labeled 1 and 2, and defines a “degree of similarity”, that is an overlap $`q`$, between them. Then one studies the system composed by the coupled replicas, with Hamiltonian $`=_1+_2ϵq`$. If $`F(ϵ)`$ is the free energy of the coupled system, then the Legendre transform $`V(q)=\underset{ϵ}{\mathrm{max}}[F(ϵ)+ϵq]`$ will represent the “effective potential” as a function of the overlap $`q`$, that is the work needed to keep the two replicas, in the configuration space, at the distance $`q`$. What one finds in $`p`$-spin models  and in structural glasses in the HNC approximation , and expects generally whenever metastable states are present in the system, that is for temperatures between the static and dynamical transition, is that the potential $`V(q)`$ presents two minima. A lower minimum, for low values of the overlap $`q`$, corresponds to the “unbounded state”, in which the two replicas are in different “valleys” of the configuration space. A second higher minimum, for high values of $`q`$, corresponding to the “bounded state”, in which the two replicas are constrained to stay in the same valley. As the free energy of each single system is independent from which valley the system is in, the gap $`\mathrm{\Delta }V`$ between the two minima can be interpreted as being equal to $`T\mathrm{\Sigma }`$, where $`\mathrm{\Sigma }`$ is the logarithm of the number of valleys, and is called complexity or configurational entropy.

We apply the effective potential theory to the model defined by (1) with $`K/J=0`$. The phase space is depicted in Fig. 1, where the continuous curve represents the line of second order transition, the black dot the tricritical point, and the dashed curve the line of first order transitions. We now introduce two real replicas, coupled in the density overlap $`q_d=\frac{1}{N}_in_i^1n_i^2`$. The Hamiltonian of the system is

$$=\underset{a=1,2}{}\left(\underset{i<j}{}J_{ij}S_i^aS_j^an_i^an_j^a+\mu \underset{i}{}n_i^a\right)ϵ\underset{i}{}n_i^1n_i^2$$
(2)
where $`a=1,2`$ is the label of the real replica. By the usual replica trick, we can write the free energy of the system as

$$F=\frac{\beta J^2}{2}(d^2+q_s^2)+\frac{\beta J^2}{n}\underset{\alpha <\beta }{}(q_{\alpha \beta }^2+p_{\alpha \beta }^2)\frac{T}{n}\mathrm{log}\text{T}r\mathrm{exp}(_{\text{eff}})$$
(3)
where $`d=n^{a\alpha }`$ is the mean density, $`q_s=S^{1\alpha }S^{2\alpha }n^{1\alpha }n^{2\alpha }`$ is the spin overlap between two directly coupled replicas, $`q_{\alpha \beta }=S^{a\alpha }S^{a\beta }n^{a\alpha }n^{a\beta }`$ and $`p_{\alpha \beta }=S^{1\alpha }S^{2\beta }n^{1\alpha }n^{2\beta }`$ are respectively the overlap matrices between replicas with the same and with different “real label”, and $`_{\text{eff}}`$ is the effective single site Hamiltonian

$`_{\text{eff}}=`$ $`\left({\displaystyle \frac{\beta ^2J^2}{2}}d+\beta \mu \right){\displaystyle \underset{a\alpha }{}}n^{a\alpha }\beta ^2J^2{\displaystyle \underset{\alpha <\beta }{}}q_{\alpha \beta }{\displaystyle \underset{a}{}}S^{a\alpha }S^{a\beta }n^{a\alpha }n^{a\beta }`$ (5)
$`2\beta ^2J^2{\displaystyle \underset{\alpha <\beta }{}}p_{\alpha \beta }S^{1\alpha }S^{2\beta }n^{1\alpha }n^{2\beta }\beta ^2J^2q_s{\displaystyle \underset{\alpha }{}}S^{1\alpha }S^{2\alpha }n^{1\alpha }n^{2\alpha }\beta ϵ{\displaystyle \underset{\alpha }{}}n^{1\alpha }n^{2\alpha }`$
In the replica symmetric approximation $`q_{\alpha \beta }=q`$ and $`p_{\alpha \beta }=p`$. We make the further approximation $`q=p`$, and obtain for the free energy

$`F=`$ $`2T\mathrm{log}2+{\displaystyle \frac{\beta J^2}{2}}(d^2+q_s^22q^2)`$ (7)
$`T{\displaystyle 𝒟z\mathrm{log}\left\{1+2e^\mathrm{\Xi }\mathrm{cosh}(\beta J\sqrt{q}z)+e^{2\mathrm{\Xi }+\beta ϵ}\left[e^\mathrm{\Omega }\mathrm{cosh}^2(\beta J\sqrt{q}z)\mathrm{sinh}\mathrm{\Omega }\right]\right\}}`$
where $`\mathrm{\Xi }=\frac{\beta ^2J^2}{2}(dq)+\beta \mu `$, $`\mathrm{\Omega }=\beta ^2J^2(q_sq)`$, and $`𝒟z=\frac{dz}{\sqrt{2\pi }}e^{z^2/2}`$.

The physical states of the system will be given by the saddle points of the free energy, which are given by the equations

$`d`$ $`=`$ $`{\displaystyle 𝒟z\frac{e^\mathrm{\Xi }\mathrm{cosh}(\beta J\sqrt{q}z)+e^{2\mathrm{\Xi }+\beta ϵ}\left[e^\mathrm{\Omega }\mathrm{cosh}^2(\beta J\sqrt{q}z)\mathrm{sinh}\mathrm{\Omega }\right]}{1+2e^\mathrm{\Xi }\mathrm{cosh}(\beta J\sqrt{q}z)+e^{2\mathrm{\Xi }+\beta ϵ}\left[e^\mathrm{\Omega }\mathrm{cosh}^2(\beta J\sqrt{q}z)\mathrm{sinh}\mathrm{\Omega }\right]}}`$ (9)
$`q_s`$ $`=`$ $`{\displaystyle 𝒟z\frac{e^{2\mathrm{\Xi }+\beta ϵ}\left[e^\mathrm{\Omega }\mathrm{cosh}^2(\beta J\sqrt{q}z)\mathrm{cosh}\mathrm{\Omega }\right]}{1+2e^\mathrm{\Xi }\mathrm{cosh}(\beta J\sqrt{q}z)+e^{2\mathrm{\Xi }+\beta ϵ}\left[e^\mathrm{\Omega }\mathrm{cosh}^2(\beta J\sqrt{q}z)\mathrm{sinh}\mathrm{\Omega }\right]}}`$ (10)
$`q`$ $`=`$ $`{\displaystyle 𝒟z\left(\frac{e^\mathrm{\Xi }\mathrm{sinh}(\beta J\sqrt{q}z)\left[1+e^{\mathrm{\Xi }+\beta ϵ+\mathrm{\Omega }}\mathrm{cosh}(\beta J\sqrt{q}z)\right]}{1+2e^\mathrm{\Xi }\mathrm{cosh}(\beta J\sqrt{q}z)+e^{2\mathrm{\Xi }+\beta ϵ}\left[e^\mathrm{\Omega }\mathrm{cosh}^2(\beta J\sqrt{q}z)\mathrm{sinh}\mathrm{\Omega }\right]}\right)^2}`$ (11)
Notice that the saddle points are neither maxima nor minima of the free energy, because $`F`$ has to be minimized with respect to $`d`$ and $`q_s`$, and maximized with respect to $`q`$. If $`ϵ=0`$ then $`q_s=q`$, $`\mathrm{\Omega }=0`$, and one obtains the saddle point equations of the uncoupled system. Differentiating the free energy $`F`$ with respect to $`ϵ`$, we can evaluate the density overlap $`q_d=n^{1\alpha }n^{2\alpha }`$,

$$q_d=\frac{F}{ϵ}=𝒟z\frac{e^{2\mathrm{\Xi }+\beta ϵ}\left[e^\mathrm{\Omega }\mathrm{cosh}^2(\beta J\sqrt{q}z)\mathrm{sinh}\mathrm{\Omega }\right]}{1+2e^\mathrm{\Xi }\mathrm{cosh}(\beta J\sqrt{q}z)+e^{2\mathrm{\Xi }+\beta ϵ}\left[e^\mathrm{\Omega }\mathrm{cosh}^2(\beta J\sqrt{q}z)\mathrm{sinh}\mathrm{\Omega }\right]}$$
(12)
In the above expression we can put $`ϵ=\mathrm{\Omega }=0`$, and find the expression of the density overlap for a system without coupling, which cannot be obtained directly from the free energy of the uncoupled system.

Now we want to compute the effective potential as a function of the density overlap $`q_d`$. The procedure to do so was previously sketched, and consists in evaluating $`V(q_d)=\underset{ϵ}{\mathrm{max}}[F(ϵ)+ϵq_d]`$. We evaluate $`V(q_d)`$ in three points of the phase space (marked by asterisks in Fig. 1). One point at the edge of the region in which spin glass solutions of the saddle point equations exist, one on the line of first order transitions, and one in between. The result is shown in Fig. 2. As soon as one enters in the region where spin glass solutions exist, a secondary minimum in the effective potential appears, signaling the existence of metastable states. We expect therefore that the dotted line of Fig. 1 will represent a line of purely dynamical transitions, below which ergodicity is broken and the relaxation function $`n_i(0)n_i(t)`$ will not decay to its equilibrium value, which in the paramagnetic phase is the density squared, but will stick at the value given by the secondary minimum of the potential.

The presence of metastable states, and of a dynamical transition at a temperature higher than the static transition, in a model in which the low temperature spin glass phase is characterized by a full RSB, seems rather peculiar. It is reasonable to think that the usual picture, valid for $`p`$-spin models, according to which at the dynamical transition the phase space is splitted into an extensive number of valleys having all the same “distance”, that is the same mutual overlap, is modified in this case. Here in fact the spin glass phase is characterized by a complex ultrametric structure, typical of the SK model. In particular one may ask if in this case the gap $`\mathrm{\Delta }V`$ in the effective potential can still be interpreted as a configurational entropy, that is the logarithm of the number of valleys in the configuration space. Finally, it would be interesting to understand how and to what extent this model can tell us something about real systems, that is structural glasses.

I thank A. Coniglio, M. Nicodemi, M. Sellitto and J.J. Arenzon for many interesting discussions. This work was partially supported by the European TMR Network-Fractals (Contract No. FMRXCT980183), MURST-PRIN-2000 and INFMPRA(HOP).


text:
# Ab-initio design of perovskite alloys with predetermined properties: The case of Pb(Sc0.5Nb0.5)O3

## Abstract

A first-principles derived approach is combined with the inverse Monte Carlo technique to determine the atomic orderings leading to prefixed properties in Pb(Sc0.5Nb0.5)O3 perovskite alloy. We find that some arrangements between Sc and Nb atoms result in drastic changes with respect to the disordered material, including ground states of new symmetries, large enhancement of electromechanical responses, and considerable shift of the Curie temperature. We discuss the microscopic mechanisms responsible for these unusual effects.

Perovskite alloys with the general form (AA′′…)(BB′′…)O3 are of great technological interest because of their large piezoelectric and dielectric responses . Examples of their applications are piezoelectric transducers and actuators, non-volatile ferroelectric memories, and microelectronic devices. Interestingly, little is known and understood about the microscopic mechanisms responsible for their very convenient properties.

One particular aspect that seems very promising both to gain fundamental knowledge and to produce even more efficient materials is the effect of atomic ordering on the properties of these alloys . For instance, one may wonder whether some specific atomic arrangements can generate drastic and useful changes with respect to the usually-grown disordered material. Examples of such changes are a large enhancement of the electromechanical responses of the system or moving its Curie temperature into a temperature range of interest. If such arrangements exist, identifying the microscopic origin of their unusual properties would definitely be of great fundamental and technological importance. Addressing these issues in an efficient manner requires (1) a computational technique able to accurately determine the properties of sufficiently large alloy configurations and providing at the same time a deep microscopic insight; and (2) the ability to find the right atomic ordering leading to a given prefixed property, i.e. to solve the so-called inverse problem .

The purpose of this Letter is to demonstrate that such an ideal situation is nowadays possible. More precisely, we have used the recently-proposed first-principles derived effective-Hamiltonian method of Ref.  to tackle the inverse problem of designing and understanding ferroelectric alloys with prefixed properties. We focused on Pb(Sc0.5Nb0.5)O3 (PSN) alloy because Sc and Nb are chemically very different. (Note that Sc and Nb belong to different columns of the Periodic Table, which makes PSN an heterovalent alloy). As a result, we expected that atomic ordering in PSN might lead to significant changes with respect to the disordered case. As a matter of fact, we found that some arrangements between the Sc and Nb atoms (i) greatly enhance the electromechanical responses, (ii) lead to currently unobserved ground states of orthorhombic and monoclinic symmetries (while the disordered material adopts a well-known rhombohedral ground state), and (iii) can considerably shift the Curie temperature. Furthermore, the analysis of our findings suggests that the mechanism responsible for these anomalous features in PSN can be understood by simple electrostatic considerations and/or simple phenomenological models. Some promising directions for future experimental and theoretical research are also provided.

Modeling PSN.– We first describe the effective-Hamiltonian approach of PSN. One starts by using the so-called “virtual crystal approximation”  to create a simple system Pb$``$B$``$O3 in which $``$B$``$ is a virtual atom involving a kind of average between Sc and Nb atoms. The relevant lattice distortions of Pb$``$B$``$O3 are described via a first-principles derived effective-Hamiltonian $`H_{\mathrm{eff}}^\mathrm{B}`$ as the one used in Ref.  to simulate the finite-temperature properties of barium titanate. This Hamiltonian is then corrected to account for the presence of the two different B atoms in the real alloy. The complete Hamiltonian has the form:

$$\begin{array}{c}H_{\mathrm{eff}}(\{𝐮_i\},\{𝐯_i\},\{\eta _l\},\{\sigma _i\})=\hfill \\ H_{\mathrm{eff}}^\mathrm{B}(\{𝐮_i\},\{𝐯_i\},\{\eta _l\})+\mathrm{\Delta }E(\{𝐮_i\},\{𝐯_i\},\{\sigma _j\}),\hfill \end{array}$$
(1)
where $`𝐮_i`$ is the local soft mode in cell $`i`$, $`\{𝐯_i\}`$ are vectors related to the inhomogeneous strains, and $`\{\eta _l\}`$ are the homogeneous strains in Voigt notation. The soft mode $`𝐮`$ is defined as the supercell average of the local modes: $`𝐮=\frac{1}{N}_i𝐮_i`$, where $`N`$ is the number of cells in the system . The variables $`\{\sigma _j\}`$ characterize the atomic configuration: $`\sigma _j=+1`$ (resp. $`1`$) indicates there is a Nb (resp. Sc) atom in cell $`j`$.

Essentially, the correction term for PSN has the form:

$$\mathrm{\Delta }E=\underset{i}{}\underset{j(i)}{}\sigma _jQ_{ji}\widehat{𝐞}_{ji}𝐮_i,$$
(2)
where $`i`$ runs over all the cells and $`j`$ over the first three nearest-neighbor shells of cell $`i`$. Here, $`\widehat{𝐞}_{ji}`$ is the unit vector joining the B-site $`j`$ to the B-site $`i`$; $`\{Q_{ji}\}`$ are the alloy-correction coefficients, which only depend on the distance between sites $`i`$ and $`j`$. Interestingly, Eq. (2) can be rewritten as the interactions between the dipole moments $`Z^{}𝐮_i`$ and radial fields $`_{ji}`$ — the latter being generated by the $`\sigma _j`$, being dependent on $`\{Q_{ji}\}`$, and acting on the neighboring cells $`i`$ —:

$$\mathrm{\Delta }E=\underset{i}{}Z^{}𝐮_i\left[\underset{j(i)}{}_{ji}\right]=\underset{i}{}Z^{}𝐮_i_i,$$
(3)
where $`_i`$ can be viewed as the total field acting on $`𝐮_i`$. We find that the first-principles calculated $`Q_{ji}`$ parameters all have negative signs. As a result, the field $`_{ji}`$ created on site $`i`$ by a Sc (resp. Nb) atom sitting on a site $`j`$ is attractive (resp. repulsive). This is consistent with an electrostatic picture of PSN resulting from the fact that Sc+3 (resp. Nb+5) atoms are negatively (resp. positively) charged with respect to the average B-atom valence of $`+4`$; i.e. to some extent $`_{ji}`$ can be viewed as the electric field created in site $`i`$ by $`\sigma _j`$.

Equations (1)-(3) indicate that there is a competition in PSN between the global tendency of the system towards the ferroelectric ground state (GS) given by $`H_{\mathrm{eff}}^\mathrm{B}`$, which is of rhombohedral symmetry with the local modes all oriented along the same $`111`$ direction , and the local tendency of the modes to align parallel to $`_i`$, which typically favors an inhomogeneous local-mode distribution and therefore a decrease of the magnitude of $`𝐮`$. Recent Monte Carlo simulations using this effective-Hamiltonian description of PSN have demonstrated its high accuracy to reproduce direct first-principles results in disordered, rocksalt-ordered (RS), and other non-trivially ordered structures . These simulations also reproduce very well the finite-temperature experimental results for disordered PSN, once the simulation temperatures are linearly rescaled so that the calculated transition temperature coincides with the measured one . On the other hand, accurately predicting the properties of samples with local compositions significantly different from the 50% average of Sc and Nb may require the use of higher-order correction terms in Eq. (2). Even so, we believe that all the trends discussed in this Letter are qualitatively correct since Eq. (2) includes the first-order (and thus the most important) correction terms.

Solving inverse problems.– The so-called inverse Monte Carlo is the most convenient computational scheme to search for the atomic ordering — i.e. the solution sample — leading to a prefixed property . It is a Monte Carlo annealing done on the compositional variables $`\sigma _j`$, and controlled by an energy functional $`𝒪[\{\sigma _j\}]`$ that is defined so that its minimum — i.e. the result of the annealing procedure — corresponds to the solution sample we pursue. For instance, one would choose $`𝒪=|d_{34}|`$ to determine the atomic configuration exhibiting the highest value of the $`d_{34}`$ piezoelectric coefficient. Note that we perform direct Monte Carlo simulations using the effective Hamiltonian of Eq. (1) to evaluate the energy functional $`𝒪[\{\sigma _j\}]`$ in every inverse Monte Carlo step. Such evaluations usually require a long computing time and, as a result, we typically restricted ourselves to the use of a $`4\times 4\times 4`$ supercell (320 atoms or 64 perovskite-like cells), which is smaller than the $`12\times 12\times 12`$ supercell of Refs.  but still much larger than those feasible by direct first-principles methods.

Tuning the GS structural properties of PSN.– We estimate the GS structural properties (denoted in the following by the “0” superscript) by computing the corresponding thermal-averaged values at 5 K. The inverse problems considered here are listed in Table I. We find that the RS sample maximizes the magnitude of the GS soft mode, which is consistent with the fact that the internal field $`_i`$ is null in any cell $`i`$ in this sample and, therefore, no inhomogeneous local-mode distribution resulting from Eq. (3) can exist. For all the other inverse problems indicated in Table I, the solution sample can be described by a single non-null $`\sigma _𝐤`$ coefficient defined as:

$$\sigma _𝐤=\frac{1}{N}\underset{j}{}\mathrm{exp}(i𝐤𝐑_j)\sigma _j,$$
(4)
where $`𝐤`$ is a vector in the first Brillouin zone and $`𝐑_j`$ is the lattice vector associated to cell $`j`$ . Our results show that the composition modulations described by $`k_\alpha =1/43/4`$ (which correspond to the “Nb-Nb-Sc-Sc” sequence along the $`\alpha `$ cartesian direction) are related to the decrease of the $`\alpha `$ component $`u_\alpha ^0`$ of the GS soft mode. This result is a direct consequence of the electrostatic-related Eq. (3), since the other one-dimensional modulations compatible with the considered $`4\times 4\times 4`$ supercell — i.e. $`k_\alpha =0`$ (“Nb-Nb-Nb-Nb” or “Sc-Sc-Sc-Sc”) and $`k_\alpha =1/2`$ (“Nb-Sc-Nb-Sc”) — result in zero internal field along the $`\alpha `$ direction and therefore in no reduction of $`u_\alpha ^0`$. Consequently, the sample defined by a single $`\sigma _𝐤`$ with $`𝐤=(1/4,1/4,1/4)`$ exhibits the rhombohedral GS with the smallest magnitude of the soft mode. Similarly, the modulations described by $`𝐤=(0,0,1/4)`$ and $`𝐤=(1/2,1/2,1/4)`$ present a “$``$” sequence of very strong internal fields $`_i`$ along the $`z`$ direction. As a result, the $`z`$ components of the local modes align parallel to these fields and the mean value $`u_z^0`$ is null. This indicates that PSN samples with an orthorhombic GS can be produced by conveniently arranging the Sc and Nb atoms. We also tried to produce a tetragonal GS but only obtained a partial success, namely a sample characterized by $`𝐤=(0,1/4,1/4)`$ with a monoclinic GS for which $`|u_x^0||u_y^0|=|u_z^0|`$. Interestingly, both orthorhombic and monoclinic phases have never been observed (yet) in PSN materials since the ground states of the two experimentally-grown samples (disordered and RS) are of rhombohedral symmetry .

Maximum electromechanical responses at low temperature.– We now consider the inverse problem of maximizing the electromechanical responses of PSN at 5 K . We focus on the piezoelectric coefficient $`d_{34}`$ and the dielectric susceptibility $`\chi _{33}`$ because we numerically find that they exhibit much larger values than any other coefficient. Interestingly, we also find that both of them reach their maximum values for the same atomic ordering. The solution sample exhibits $`\chi _{33}=8319`$ and $`d_{34}=4868`$ pC/N, to be compared with the values of 325 and 90 pC/N obtained in disordered PSN. A striking enhancement of the electromechanical response can thus be achieved by rearranging the Sc and Nb atoms. (We also obtained other samples with $`\chi _{33}>4000`$ and $`d_{34}>2000`$ pC/N, all sharing the structural features that we discuss in the following.) The GS of the solution sample is characterized by $`𝐮^0=(0.544,0.567,0.076)`$, where the main feature is the small value of the $`z`$ component of the soft mode. This reflects the fact that the $`\sigma _𝐤`$ coefficients that contribute most to the atomic ordering (around 70% of the total) have $`k_z=1/4`$. Our results thus suggest that samples with a “$`++`$” composition modulation in the $`z`$ direction, where “$`+`$” (resp. “$``$”) stands for a plane that is rich in Nb (resp. Sc), will simultaneously exhibit small values of $`u_z^0`$ and large values of $`\chi _{33}`$ and $`d_{34}`$ at low temperatures. However, this modulation cannot be too strong since, for instance, we numerically found that the atomic configuration characterized by the single $`𝐤=(0,0,1/4)`$ (for which $`u_z^0=0`$, see Table I) has small electromechanical responses.

We now discuss the relationship between the small value of $`|u_z^0|`$ and the large electromechanical responses. Fig. 1 shows the probability distribution of $`u_x`$ and $`u_z`$ corresponding to the Monte Carlo simulation at 5 K of the solution sample, as well as the free-energy map of the system around $`𝐮^0`$. This latter is calculated out of the probability distribution as done in Ref. . The two main features of Fig. 1 are (1) the broadness of the $`u_z`$ distribution, which is responsible for the large value of $`\chi _{33}`$ (see Ref. ) and (2) the flatness of the corresponding free-energy well, which indicates the ease of rotating the polarization around $`u_z^0`$. Interestingly, both features are typically associated with a small value of $`|u_z^0|`$ within the usual Devonshire-Landau description of ferroelectric perovskites . Moreover, we numerically find that $`u_z`$ and $`\eta _4`$ are strongly coupled in our effective Hamiltonian for PSN, which implies that their thermal fluctuations are very correlated. Consequently, according to Ref. , a large value of $`\chi _{33}`$ typically corresponds to a large value of $`d_{34}`$, as we indeed find.

Tuning the finite-temperature properties of PSN.– Studying PSN at finite (relatively-high) temperatures requires long Monte Carlo simulations in order to obtain good statistics. Thus, due to computational limitations we decided not to directly solve finite-temperature (FT) inverse problems. However, some links between the GS and FT properties are expected on the basis of well-known models that share some essential features with our effective-Hamiltonian description of PSN. More precisely, the so-called $`\varphi ^4`$ model  suggests the following relationship: the larger the value of $`|u^0|`$, the deeper the corresponding free-energy well at 0 K and, therefore, the higher the Curie temperature $`T_c`$. Hence, we thought it was pertinent to check whether the Curie temperature of PSN could be tuned by changing the value of the GS soft mode. Consequently, we studied the FT behavior of four samples that all undergo a single cubic to rhombohedral transition but exhibit different values for their GS soft modes. Our results, shown in Fig. 2, confirm that it is indeed possible to move the Curie temperature to a temperature range of interest by designing alloys with a prefixed value of $`|u^0|`$. It is remarkable that relatively-small changes in the GS soft mode result in such large variations in $`T_c`$. Interestingly, the transition of the sample with the smallest $`|u^0|`$ takes place near room temperature. It is important to notice that one can tune the FT electromechanical responses of the system by moving $`T_c`$, because these responses reach their maximum values around the paraelectric to ferroelectric transition.

In summary, we have implemented a method to solve inverse problems in perovskite alloys by using the first-principles effective-Hamiltonian approach of Ref. , and applied it to PSN. We predict that the system properties (GS symmetry, electromechanical responses, Curie temperature, etc.) can be greatly modified by varying the atomic ordering. Furthermore, all the anomalous features we predict can be essentially understood by simple electrostatic considerations and/or simple phenomenological models. As a result, (1) we have indicated which kinds of alloy structures may lead to currently-unobserved structural phases as well as to giant dielectric and piezoelectric responses, and (2) we expect similar results to occur in other heterovalent perovskite alloys. Finally, we hope that our predictions will stimulate attempts at experimental growth of such samples.

This work is supported by Office of Naval Research Grants No. N00014-00-1-0542 and No. N00014-01-1-0600, National Science Foundation Grant DMR-9983678 and Arkansas Science and Technology Authority Grant N99-B-21. We thank Aaron George, Alberto García, and J. M. Pérez-Mato for comments and discussions. J. I. further thanks the financial support of the Government of the Basque Country and the hospitality of all the members of the Physics Department of the University of Arkansas.