text:
# 1 The International Gravitational Event Collaboration

## 1 The International Gravitational Event Collaboration

In table I there is a list of the five cryogenic resonant antennas now in operation. These detectors are aluminum (niobium in the case of the Australian one) cylindrical bars, all equipped with resonant transducers, that give them two narrow bands of detection of a few hertz around the two near coupled modes of the bar-transducer system, at frequencies of about 900 Hz (about 700 Hz for Niobe). These are called the detection bands of the antenna; in some detectors they can constitute a single wider band.

The orientation of the bars was chosen in order to achieve the maximum parallelism of the bars.

In order to coordinate and correlate the data produced by these detectors, in July 1997, all the groups that operate them have undersigned a data exchange protocol named ”International Gravitational Event Center” (”IGEC”). It is based on some technical and policy rules. The main points are:

\- The data consist of ”candidate events” and service information.

\- The minimum information about each event consist of the time of the event maximum, the amplitude in units of standard burst strain, the duration and the mean detector noise at the time of the event.

\- Additional fields, giving further information on the candidate events (e.g. the shape and/or parameters relating the shape) may be provided.

\- Time accuracy will be at least 0.1 s.

\- Each group will set up a site (ftp, www,…) on which the data will be continuously available in agreed formats, with an updating rate not exceeding one day.

The full text of the agreement, together with other data exchange documents, is posted at http://grwav1.roma1.infn.it/gwda/de/de.htm .

In this paper I will show how the candidate events are produced by the Rome group (that is operating the two cryogenic resonant antennas Explorer and Nautilus) and how the IGEC data can be analyzed.

| Antenna | Institution | City |
| --- | --- | --- |
| Allegro | LSU | Baton Rouge |
| Auriga | LNL (INFN-AURIGA) | Legnaro |
| Explorer | CERN (INFN-ROG) | Geneva |
| Nautilus | LNF (INFN-ROG) | Frascati |
| Nïobe | UWA | Perth |

## 2 How the events are obtained

The simplified modelisation of the detector as a linear system is normally very good and the gravitational signal is simply added to the noise; nevertheless the detection of short gravitational pulses (with unknown shape) in the data of gravitational wave antennas is not an easy task. The main reasons are:

\- the low signal-to-noise ratio (SNR) and the rarity of the expected pulses

\- the ignorance of their shape

\- the non-stationarity of the noise of the detectors and the presence of many spurious events.

As regards the shape of the pulses, resonant antenna people normally consider as a standard pulse a delta function (or the part of that signal in the detector band) that, in the data at the output of the antenna, becomes a known waveform. The noise is considered gaussian and stationary (also if some consider it slowly non-stationary), completely described by the noise power spectrum.

With these assumptions, the problem becomes the classical problem of detecting a known waveform in gaussian noise, that is optimally solved by the ”matched filter”. This filter is a non-causal linear system that has the property that at its output the signal-to-noise ratio (i.e. the ratio between the square of the maximum of the signal and the variance of the noise) is maximized.

This means that, if we normalize the filter in order to have the same amplitude of the maximum of the input and output signal, the output noise (that is also gaussian because of the linearity of the filter), has the minimum variance, so the probability that the noise samples exceed the threshold is minimized.

There are many ways to express the equation of the matched filter. In the frequency domain it is

$$M(j\omega )=\frac{F^{}(j\omega )}{S(\omega )}$$
(1)
where $`F^{}(j\omega )`$ is the complex conjugate of the Fourier transform of the response of the antenna to the standard pulse and $`S(\omega )`$ is the noise power spectrum.

As a matter of fact, each group implement practically the filter in one or more different ways. This because:

\- there are different ways of taking data (high frequency sampling, aliased sampling, two lock-ins acquisition, single central lock-in,…)

\- the use of frequency domain procedures or time domain procedures

\- the use of adaptive or non-adaptive procedures

\- the threshold mechanism

\- particular procedures of event features estimation and/or spurious event detection.

In our case (Ref ), the most advanced procedure that we use is based on a high frequency sampling (5000 Hz), from which we extract a band of 40 Hz centered at 816 Hz, containing the detection bands.

Other bands are also analyzed in order to monitor the behavior of the antenna and the presence of local disturbances, but this will not be discussed here.

On these data we implement three adaptive matched filters (and also filters of other types), followed by a two state adaptive threshold mechanism.

## 3 The adaptive matched filters

The noise power spectrum of our antennas is not stationary: both the electric noise, that normally gives a wide-band and a narrow band contributions, and the seismic noise, that gives normally a narrow band contribution, change with time. Moreover often extra noise peaks appear in the spectrum, due to local disturbances. These bands change slowly in amplitude and in frequency. Almost all these non-stationarities have time constants of at least a few hours. The presence of the non-stationarities can be seen as the fact that the detector has a time-varying sensitivity.

Also the resonance frequencies of the antennas can slowly change, mainly because of the change in the transducer polarization voltage. In order to have a good filter in presence of these varying features of the detectors, we use an adaptive approach: we apply the filter in the frequency domain, using the expression of eq.1 and use, for the power spectrum, a first-order auto-regressive sum of the periodograms, with a time constant of one hour and with an updating rate of about one minute. The filter is so recomputed about every minute; this is our ”basic” adaptive filter.

In presence of a big disturbance, the basic filter estimated spectrum is ”dirtied” and can remain dirty for a long time (also many time constants) after the end of the disturbance, producing filters that are not optimal for the ”clean” data. To solve this problem, we have implemented a matched filter that uses only ”clean” periodograms for the estimation of the power spectrum.

We studied also another policy that gets good results in the case that the disturbed period can be seen as a low sensitivity period. It is based on the reduction of the auto-regressive time constant in the case of highly disturbed spectra.

In order to evaluate the performances of the different filters, we add to the acquired data fictitious ”theoretical” pulses and analyze the resulting SNR; this is done every half an hour, disturbing the real data for less than the $`0.1\%`$ of the time. We chose as ”official filter” the filter that gives the better results at that time. This procedure provides also a sort of Monte Carlo check procedure for the whole system and can be used also to analyze the behavior in case of non-delta-pulse events.

## 4 The threshold mechanism

### 4.1 The adaptive threshold

Because of the non-stationarity of the antenna noise, also at the output of one matched filter (or of any linear filter) the noise is not stationary, i.e. the variance of the noise changes (slowly) with the time; this means that the sensitivity of the detector changes with the time. If we put a fixed threshold on these data, we have more (spurious) candidate events when higher is the noise (and lower is the sensitivities) and this can highly worsen the statistics. So we must change the threshold with the time.

This can be done in various ways, for example choosing a fixed number N and taking the N highest events of each hour.

We use a different procedure, based on an ”adaptive threshold” defined in the following way.

Let $`x_i`$ be the filtered data samples. We estimate the background statistics by computing the auto-regressive mean of the absolute value and of the square of $`x_i`$

$$m_i=(1w)|x_i|+wm_{i1}$$
(2)
$$q_i=(1w)x_i^2+wq_{i1}$$
(3)
with

$$w=e^{\frac{\mathrm{\Delta }t}{\tau }}$$
(4)
where $`\mathrm{\Delta }t`$ is the sampling time and $`\tau `$ is the ”memory” of the auto-regressive mean (we normally chose $`\tau =600s`$).

Then we define the standard deviation

$$\sigma _i=\sqrt{q_im_i^2}$$
(5)
and the threshold $`\theta `$ is not set on $`|x_i|`$, but on the critical ratio of $`|x_i|`$ given by

$$z_i=\frac{|x_i|m_i}{\sigma _i}$$
(6)
This procedure was developed for the general case of $`x_i`$; if $`x_i`$ is simply zero mean, as is in the case of the matched filter, the algorithm for estimating the adaptive threshold can be simplified, using just the estimation $`q_i`$ and computing the critical ratio as $`z_i=|x_i|/\sqrt{q_i}`$.

In case of very large value of $`\sigma _i`$, we can reduce the memory $`\tau `$ in order to reduce the ”blinding” effect at the end of large disturbances.

### 4.2 The event definition: the two state event machine

In order to define candidate events, let us suppose that we have chosen an adaptive threshold $`\theta `$ and a dead time (that is the minimum time between two different events; it depends on the apparatus, the noise and the expected signal: we normally use 3 s).

Then we use an easy two-state (0 and 1) mechanism that we call the ”event machine”. The algorithm performed is the following.

\- The machine is normally set in the state $`0`$.

\- When the signal goes over the threshold, it changes to state $`1`$ and an event begins.

\- The state changes to $`0`$ after the signal has remained below the threshold for a time longer than the dead time.

\- The ”duration” of the event is the duration of the state $`1`$, subtracted the dead time.

A simplified model of this algorithm (Ref ) is a two-state Markov chain (in the discrete time). Its transition matrix is

$$\left(\begin{array}{cc}1p_{01}& p_{01}\\ p_{10}& 1p_{10}\end{array}\right)$$
(7)
where $`p_{01}`$ and $`p_{10}`$ are the transition probability for the two states. We can easily compute the probabilities of the two states as

$$p_0=\frac{p_{10}}{p_{01}+p_{10}}$$
(8)
$$p_1=1p_0$$
(9)
and the mean length (in unit of sampling time)

$$\overline{L}=\frac{1}{p_{10}}$$
(10)
### 4.3 The event density

We can define the event density $`\lambda `$ as the number of the events per unit time and then the candidate events production can be modeled by a Poisson process with parameter $`\lambda `$.

Obviously the value of $`\lambda `$ depends on the value of the adaptive threshold $`\theta `$. The functional dependence of $`\lambda `$ on $`\theta `$ depends on the distribution of the filtered data (or, more precisely, on its tail), that only theoretically is gaussian; in practice it depends strongly on the local disturbances.

This is what finally limits the sensitivity of a single antenna. However, because the disturbances give heavy tails in the distributions, reducing the threshold (and so enhancing the sensitivity) in this region gives a not big increment on $`\lambda `$. It is not so for the gaussian distribution (that has light tails).

## 5 The coincidences

If we consider the case of $`N`$ antennas, each defined with a $`\lambda _i`$, and if we choose a coincidence window of duration $`t_w`$, in the hypothesis of uncorrelated data, that is that of the background noise event, we have an expected density of ”casual” coincidences given by

$$\lambda ^{(N)}=t_w^{N1}\underset{i=1}{\overset{N}{}}\lambda _i$$
(11)
$`t_w`$ is chosen depending on the apparatuses, the time precision and the light time delay between the antennas; good values ranges between $`0.1`$ and $`1s`$. In table 2 there are some values for the coincidence densities (expressed in number of coincidences per day) for the case that all $`\lambda _i`$’s are equal to the value $`\lambda `$ and $`t_w=0.1s`$.

|  | $`\lambda =\mathrm{𝟏𝟎𝟎}\mathrm{𝐞𝐯𝐞𝐧𝐭𝐬}/\mathrm{𝐝𝐚𝐲}`$ | $`\lambda =\mathrm{𝟏𝟎𝟎𝟎}\mathrm{𝐞𝐯𝐞𝐧𝐭𝐬}/\mathrm{𝐝𝐚𝐲}`$ |
| --- | --- | --- |
| $`𝐍=\mathrm{𝟐}`$ | $`1.1510^2`$ | $`1.15`$ |
| $`𝐍=\mathrm{𝟑}`$ | $`1.3410^6`$ | $`1.3410^3`$ |
| $`𝐍=\mathrm{𝟒}`$ | $`1.5510^{10}`$ | $`1.5510^6`$ |
| $`𝐍=\mathrm{𝟓}`$ | $`1.7910^{14}`$ | $`1.7910^9`$ |

In order to evaluate the statistical significance of the found coincidences, one should consider the appropriate Poisson statistics.

Anyway, the Poisson model is normally not a good approximation, because of all the non-stationarities in the detection process. In practice the $`\lambda `$’s are functions of the time. In particular, there are

\- periodicities (e.g. the solar day)

\- aggregation (the disturbances often are clustered in time)

\- holes in the data.

Also the signals (the ”true” events) are not expected to be uniformly distributed in time, because of the radiation pattern of the antennas and the non-uniformity of the space distribution of the biggest sources, that causes a sidereal day modulation of the detection probability.

To have an efficient evaluation method in presence of this problems, in the case of two antennas, Weber introduced a non-parametric procedure, that we call ”pulse correlation”.

### 5.1 The pulse correlation

The pulse correlation method is based on the evaluation of the casual (”background”) coincidence rate for two antennas obtained by adding a bias time $`\tau `$ to the eventts of one antenna. $`\tau `$ is normally set equal to $`kt_w`$, with $`k`$ integer and $`nkn`$ (e.g., $`n`$ can be set equal to 1000).

Taken a period of time $`T`$ (e.g. one day, one month,…) the number of coincidences $`C(\tau )`$ are computed. If we take the mean value $`\mu _C`$ of $`C(\tau )`$ excluding $`C(0)`$, the estimated number of ”true” coincidences is given by

$$\stackrel{~}{C}=C(0)\mu _C$$
(12)
if $`\stackrel{~}{C}0`$. We can define $`\lambda _C=\mu _C/T`$ as the estimated chance coincident event density (analogous to $`\lambda ^{(2)}`$ for the equation 11).

Because of the non-stationarities, the expected shape of $`C(\tau )`$ for $`\tau 0`$ is not uniform and the evaluation of $`\mu _C`$ deserves some cares.

I would note the analogy between the pulse correlation function and the classical cross-correlation function.

To evaluate the chance probability of having $`C(0)`$, one can use one of the two following methods:

a) use the Poisson statistics with parameter $`\mu _C`$

b) compute the histogram of $`C(\tau )`$ and evaluate the ”frequentistic” probability of having $`C(0)`$.

Both methods must be used with care, the first because of the uncertainty in the evaluation of $`\mu _C`$ and the second because it assumes that the event production is stationary on times of the order of $`nt_w`$.

A particular use of the pulse correlation is in correlating the events of an antenna with themselves (”pulse auto-correlation”). This method of analysis can be used to identify the non-stationarities of the events: in particular aggregation and periodicities.

The pulse correlation can be applied on data selected by some ”single event” rules (e.g. the amplitude or the time of occurrence, that can be a particular solar or sidereal hour) or ”coincidence” rules (e.g. consider a coincidence only if the amplitude and/or the duration of the two events is about the same).

In the evaluation of the probabilities, particular care must be taken on any choice made ”a posteriori”.

### 5.2 The pulse correlation in the case of more than two antennas

How can we generalize the pulse correlation in the case of $`N>2`$ antennas ? Consider the following two methods:

a) the couple pulse correlation, obtained summing the pulse correlation of all the $`N(N1)/2`$ couples of antennas. We obtain a function $`C_N(\tau )`$. With these method there is a ”natural” weighting of multiple coincidences. In fact a coincidence between only two antennas is considered just once, a coincidence between 3 antennas is considered 3 times , a 4 antennas coincidence it is considered 6 times and a 5 antennas coincidence is considered 10 times. The background is the sum of all the backgrounds.

b) the multiple pulse correlation $`C(\tau _1,\tau _2,\mathrm{},\tau _{N1})`$, obtained by adding to the events of each of the first $`N1`$ antennas a different bias time. The delay variables $`\tau _i`$ can be chosen equal to $`kt_w`$, with $`nkn`$ (e.g., $`n`$ can be set equal to 10; in this case, with five antennas, we have 194481 delays).

In searching for coincidences for more then two antennas, it must be taken into account that, in case of not too much big gravitational events, the probability that an event is overlooked by a detector is not negligible. This because of the presence of the additive noise (see, for example, Ref ), the difference in the frequencies, the not perfect parallelism and the difference in sensitivities, that, also if are basically similar, change in time depending on the local noises. So methods like the multiple pulse correlation are good only in case of huge events and, for small events, would produce a false dismissal probability near to 1.

More complex procedures were presented in Ref . With these procedures, with five detectors it is possible to estimate also the position of the source in the sky and the polarization of the gravitational pulse, together with a better estimate of its energy. Also the rejection of spurious events is enhanced. Unfortunately these methods get poor results with parallel antennas.

## 6 Another perspective

Because of the rarity of the gravitational events, it is important to not overlook the presence of real events in the data. So any method used to reduce the false detection probability, should not enhance too much the false dismissal probability. This means that we should reject a candidate event only after a careful analysis.

On the other hand, also very low probability coincidence events must be carefully checked for consistency, e.g. to compare the amplitudes and the lengths in the different antennas or to check in the other antennas, that did not detect events at that time, if there is ”something” (i.e. an event under the threshold).

As it is shown in table 2, the multiple coincidence operation (2 or 3 antennas on 5) reduces strongly the number of candidate events; this can be done automatically. For the surviving candidates we should analyze carefully (i.e. not automatically) the outputs of all the antennas in operation, the signal shapes and the auxiliary channels.

It is important to have a data-base of all the events, with some characteristics (e.g. the ”color”, i.e. some spectral parameters, the length, the shape,…).

A very important information to take into account are other impulsive astrophysical events, e.g. the gamma ray bursts and the neutrinos bursts detected in underground experiments. These must constitute a complementary data base to be analyzed routinely together with the candidate events. Also supernova surveys can be useful, in order to correctly interpret the other detector results, also if they don’t give precise events time.

After all, the ”true decision” on the ”promotion” of a candidate event will be made not by an algorithm, but by a long discussion.


format: text

text:
# Partition Function Zeros and Finite Size Scaling of Helix-Coil Transitions in a Polypeptide

## Abstract

We report on multicanonical simulations of the helix-coil transition of a polypeptide. The nature of this transition was studied by calculating partition function zeros and the finite-size scaling of various quantities. New estimates for critical exponents are presented.

A common, ordered structure in proteins is the $`\alpha `$-helix and it is conjectured that formation of $`\alpha `$-helices is a key factor in the early stages of protein-folding . It is long known that $`\alpha `$-helices undergo a sharp transition towards a random coil state when the temperature is increased. The characteristics of this so-called helix-coil transition have been studied extensively , most recently in Refs. . They are usually described in the framework of Zimm-Bragg-type theories  in which the homopolymers are approximated by a one-dimensional Ising model with the residues as “spins” taking values “helix” or “coil”, and solely local interactions. Hence, in such theories thermodynamic phase transitions are not possible. However, in preliminary work  it was shown that our all-atom model of poly-alanine exhibits a phase transition between the ordered helical state and the disordered random-coil state. It was conjectured that this transition is due to long range interactions in our model and the fact that it is not one-dimensional: it is known that the 1D Ising model with long-range interactions also exhibits a phase transition at finite $`T`$ if the interactions decay like $`1/r^\sigma `$ with $`1\sigma <2`$ . Our aim now is to investigate this transition in the frame work of a critical theory by means of the finite size scaling (FSS) analysis of partition function zeros. Analysis of partition function zeros is a well-known tool in the study of phase transitions, but was to our knowledge never used before to study biopolymers.

For our project, the use of the multicanonical algorithm  was crucial. The various competing interactions within the polymer lead to an energy landscape characterized by a multitude of local minima. Hence, in the low-temperature region, canonical simulations will tend to get trapped in one of these minima and the simulation will not thermalize within the available CPU time. One standard way to overcome this problem is the application of the multicanonical algorithm  and other generalized-ensemble techniques  to the protein folding problem . For poly-alanine, both the failure of standard Monte Carlo techniques and the superior performance of the multicanonical algorithm are extensively documented in earlier work .

In the multicanonical algorithm  conformations with energy $`E`$ are assigned a weight $`w_{mu}(E)1/n(E)`$. Here, $`n(E)`$ is the density of states. A simulation with this weight will lead to a uniform distribution of energy:

$$P_{mu}(E)n(E)w_{mu}(E)=\mathrm{const}.$$
(1)
This is because the simulation generates a 1D random walk in the energy, allowing itself to escape from any local minimum. Since a large range of energies are sampled, one can use the reweighting techniques  to calculate thermodynamic quantities over a wide range of temperatures by

$$<𝒜>_T=\frac{{\displaystyle 𝑑x𝒜(x)w^1(E(x))e^{\beta E(x)}}}{{\displaystyle 𝑑xw^1(E(x))e^{\beta E(x)}}},$$
(2)
where $`x`$ stands for configurations.

It follows from Eq. 1 that the multicanonical algorithm allows us to calculate estimates for the spectral density:

$$n(E)=P_{mu}(E)w_{mu}^1(E).$$
(3)
We can therefore construct the partition function from these estimates by

$$Z(\beta )=\underset{E}{}n(E)u^E,$$
(4)
where $`u=e^\beta `$ with $`\beta `$ the inverse temperature, $`\beta =1/k_BT`$. The complex solutions of the partition function determine the critical behavior of the model. They are the so-called Fisher zeros , and correspond to the complex extension of the temperature variable.

Our investigation of the helix-coil transition for poly-alanine is based on a detailed, all-atom representation of that homopolymer, and goes beyond the approximations of the Zimm-Bragg model . The interaction between the atoms was described by a standard force field, ECEPP/2, (as implemented in the KONF90 program ) and is given by:

$`E_{tot}`$ $`=`$ $`E_C+E_{LJ}+E_{HB}+E_{tor},`$ (5)
$`E_C`$ $`=`$ $`{\displaystyle \underset{(i,j)}{}}{\displaystyle \frac{332q_iq_j}{ϵr_{ij}}},`$ (6)
$`E_{LJ}`$ $`=`$ $`{\displaystyle \underset{(i,j)}{}}\left({\displaystyle \frac{A_{ij}}{r_{ij}^{12}}}{\displaystyle \frac{B_{ij}}{r_{ij}^6}}\right),`$ (7)
$`E_{HB}`$ $`=`$ $`{\displaystyle \underset{(i,j)}{}}\left({\displaystyle \frac{C_{ij}}{r_{ij}^{12}}}{\displaystyle \frac{D_{ij}}{r_{ij}^{10}}}\right),`$ (8)
$`E_{tor}`$ $`=`$ $`{\displaystyle \underset{l}{}}U_l\left(1\pm \mathrm{cos}(n_l\chi _l)\right).`$ (9)
Here, $`r_{ij}`$ (in Å) is the distance between the atoms $`i`$ and $`j`$, and $`\chi _l`$ is the $`l`$-th torsion angle. Note that with the electrostatic energy term $`E_C`$ our model contains a long range interaction neglected in the Zimm-Bragg theory . Since one can avoid the complications of electrostatic and hydrogen-bond interactions of side chains with the solvent for alanine (a nonpolar amino acid), explicit solvent molecules were neglected. Chains of up to $`N=30`$ monomers were considered. We needed between 40,000 sweeps ($`N=10`$) and 500,000 sweeps ($`N=30`$) for the weight factor calculations by the iterative procedure described in Refs. . All thermodynamic quantities were estimated from one production run of $`N_{sw}`$ Monte Carlo sweeps starting from a random initial conformation, i.e. without introducing any bias. We chose $`N_{sw}`$=400,000, 500,000, 1,000,000, and 3,000,000 sweeps for $`N=10`$, 15, 20, and 30, respectively.

For our analysis of the partition function zeros we first divide the energy range into intervals of lengths $`0.5`$ kcal/mol. Equation 4 becomes now a polynomial in the variable $`u`$ and can be easily solved with MATHEMATICA to obtain all complex zeros $`u_j^0`$ ($`j=1,2,\mathrm{}`$). For the case of $`N=10`$ we also repeated the calculation of the zeros for energy bin sizes $`1.0`$ kcal/mol and $`0.25`$ kcal/mol. The changes in the zeros were smaller than the statistical errors. The effect of the energy bin size on the zeros is also discussed in Ref. . Figure 1 shows the distribution of the zeros for $`N=30`$ and provides already strong evidence for a singularity on the real axis: in the case of the (analytic) Zimm-Bragg theory the zeros would be located solely on the negative real $`u`$-axis . We summarize in Table I the leading zeros for each of the four chain lengths, where we have used the mapping $`u=e^{\beta /2}`$ due to our binning procedure.

The FSS relation by Itzykson et al.  for the leading zero $`u_1^0(N)`$,

$$u_1^0(N)=u_c+AN^{1/d\nu }[1+O(N^{y/d})],y<0$$
(10)
shows that the distance from the closest zero $`u_1^0`$, to the infinite-chain critical point $`u_c=e^{\beta _c/2}`$ on Re($`u`$) axis, scales with a relevant linear length $`L`$, which we translated as $`N^{1/d}`$ in the above equation. Here, $`\beta _c`$ is the inverse critical temperature of the infinite long polymer chain and $`y`$ is the correction to scaling exponent. We remark that, unlike in the Zimm-Bragg model, we have no theoretical indication to assume $`d`$ as a particular integer geometrical dimension and report therefore estimates for the quantity $`d\nu `$.

For sufficiently large $`N`$, the exponent $`d\nu `$ can be obtained from the linear regression

$$\mathrm{ln}|u_1^0(N)u_c|=\frac{1}{d\nu }\mathrm{ln}(N)+a.$$
(11)
This relation requires an accurate estimate for $`u_c`$. Therefore, we prefer to calculate our estimates for $`d\nu `$ from the corresponding relation with $`|u_1^0u_c|`$ replaced by its imaginary part $`\mathrm{Im}u_1^0`$. Including chains of all lengths, $`N=1030`$, this approach leads to $`d\nu =0.93(5)`$, with a goodness of fit $`Q=0.48`$. Figure 2 displays the corresponding fit. Omitting the smallest chain, i.e. restricting the fit to the range $`N=1530`$, does not change the above result. We obtain now $`d\nu =0.93(7)`$, with $`Q=0.22`$. This indicates that the $`d\nu `$ determination is stable over the studied chains and therefore, the correction exponent $`y`$ can be disregarded in face of the present statistical error.

Considering the real part of the leading zeros given in Table I, $`\mathrm{Re}(\beta _1^0(N))=\mathrm{ln}\{[\mathrm{Re}u_1^0(N)]^2+[\mathrm{Im}u_1^0(N)]^2\}`$, we can derive the critical temperature through the following FSS fit ,

$$\mathrm{Re}(\beta _1^0(N))=\beta _c+bN^{1/d\nu }.$$
(12)
We obtain $`\beta _c=0.906(12)`$ ($`Q=0.005`$) for the range $`N=1030`$, and $`\beta _c=0.929(14)`$ ($`Q=0.63`$) for $`N=1530`$. The last and more acceptable estimate, corresponds to $`T_c(\mathrm{})=541(8)`$ K.

A stronger version for the relation (10) considers that the next zeros $`u_j^0(N)`$, should also satisfy a scaling relation ,

$$|u_j^0(N)u_c|\left(\frac{j}{N}\right)^{1/d\nu },$$
(13)
where $`j`$ labels the zeros in order of increasing distance from $`u_c`$. This relation is expected to be satisfied for large $`j`$ and allows for an independent check of the estimate for our exponent $`d\nu `$. The scaling plot in Fig. 3 for the roots closest to the critical point $`u_c`$ demonstrates that the assumed scaling relation is indeed observed for our data as $`N`$ increases and consistent with our estimate of the exponent $`d\nu `$.

Our results for the critical temperature and critical exponent can be compared with independent estimates obtained from FSS of the specific heat:

$$C_N(T)=\beta ^2(<E^2(T)><E(T)>^2)/N.$$
(14)
Defining the critical temperature $`T_c(N)`$ as the position where the specific heat $`C_N(T)`$ has its maximum, we can again calculate the critical temperature by means of Eq. 12. With the values in Table II we obtain $`T_c(\mathrm{})=544(12)`$ K, which is consistent with the value obtained from the partition function zeros analysis, $`T_c(\mathrm{})=541(8)`$ K. Choosing $`T_1(N)`$ and $`T_2(N)`$ such that $`C(T_1)=1/2C(T_c)=C(T_2)`$, we have the following scaling relation for the width $`\mathrm{\Gamma }_C(N)`$ of the specific heat ,

$$\mathrm{\Gamma }_C(N)=T_2(N)T_1(N)N^{1/d\nu }.$$
(15)
Using the above equation and the values given in Table II, we obtain $`d\nu =0.98(11)`$ $`(Q=0.9)`$ for chains of length $`N=15`$ to $`N=30`$, i.e. omitting the shortest chain. This value is in agreement with our estimate $`d\nu =0.93(5)`$, obtained from the partition function zero analysis. Including $`N=10`$ leads to $`d\nu =1.19(10)`$, but with a less acceptable fit ($`Q=0.1`$). The analysis of partition function zeros seems also to be more stable than one relying on Eq. 15. No significant change in $`d\nu `$ was observed when the data from ref.  (which relied on much smaller number of Monte Carlo sweeps) were used in the partition function zeros analysis, while Eq. 15 leads for this reduced statistics to an estimate for $`d\nu =1.9`$.

Through the scaling relation for the peak of the specific heat, we can evaluate yet another critical exponent, the specific heat exponent $`\alpha `$, by:

$$C_N^{max}N^{\alpha /d\nu }.$$
(16)
In particular, with the values for $`C_N^{max}`$ as given in Table II, we obtain $`\alpha =0.86(10)`$. The scaling plot for the specific heat is shown in Fig. 4: curves for all lengths of the poly-alanine chains nicely collapse on each other indicating the scaling of the specific heat and the reliability of our exponents. It worths to note that our estimates for $`d\nu `$ and $`\alpha `$, as obtained from the finite size scaling of the specific heat, obey within the errorbars the hyperscaling relation $`d\nu =2\alpha `$.

It is well-known that renormalization-group fixed point picture leads to a critical exponent $`d\nu =1`$, $`\alpha =1`$ and $`\gamma =1`$ for a first-order phase transition . Our estimate $`d\nu =0.93(5)`$ for the correlation exponent deviates from unity and rather indicates that the ‘helix-coil-transition” is a strong second order transition. However, the errorbars are such that a first order phase transition cannot be excluded. Our values for the specific heat exponent $`\alpha =0.86(10)`$ and the susceptibility exponent $`\gamma =1.06(14)`$ (data not shown) are consistent with a first order phase transition, but also not conclusive. A common way to evaluate the order of a phase transition is by means of the Binder energy cumulant ,

$$b(T,N)=1\frac{<E^4(T,N)>}{3<E^2(T,N)>^2}.$$
(17)
For a second order phase transition one would expect that the minimum of this quantity $`b(T_{min},\mathrm{})`$ approaches $`2/3`$. Here $`T_{min}`$ defines the temperature where the cumulant reaches its minimum value and $`b(T_{min},\mathrm{})=lim_N\mathrm{}b(T,N)`$. With the present values of Table II we find the infinite volume extrapolation $`b(T_{min},\mathrm{})=0.23(13)`$ $`(Q=0.12)`$, for the range $`N=1530`$, which is consistent with a first order phase transition. However, we cannot exclude the possibility of a second order phase transition because the energy cumulant scales with the maximum of specific heat , $`b(T,N)N^{\alpha /d\nu 1}`$, and the true asymptotic limit is reached only for rather large chains due to the value of $`\alpha /d\nu `$. In fact, the straight line fit for the range $`N=1030`$ is less consistent with our data ($`Q0.001`$). Hence, we conclude that our results seem to favor a (weak) first order phase transition, but are not precise enough to exclude the possibility of a second order phase transition.

To summarize, we have used a common technique for investigation of phase transitions, analysis of the finite size scaling of partition function zeros, to evaluate the helix-coil transition in an all-atom model of poly-alanine. Although our results are due to the complexity of the simulated model not precise enough to determine the order of the phase transition, we have demonstrated that the transition can be described by a set of critical exponents. Hence, we have shown for this example that structural transitions in biological molecules can be described within the frame work of a critical theory.

Acknowledgements:
Financial supports from FAPESP and a Research Excellence Fund of the State of Michigan are gratefully acknowledged.

Figure Captions:

1. Partition function zeros in the complex $`u`$ plane for $`N=30`$. For the Zimm-Bragg model the zeros would be located solely on the negative real $`u`$ axis.
2. Linear regression for -ln (Im$`u_1^0(N)`$) in the range $`N=1030`$.
3. Scaling behavior of the first $`j`$ complex zeros closest to $`u_c=0.6284`$, for chain lengths $`N=10,15,20`$ and 30.
4. Scaling plot for the specific heat $`C_N(T)`$ as a function of temperature $`T`$, for poly-alanine molecules of chain lengths $`N=10,15,20,`$ and $`30`$.


text:
# Nucleation theory and the phase diagram of the magnetization-reversal transition

## Abstract

The phase diagram of the dynamic magnetization-reversal transition in pure Ising systems under a pulsed field competing with the existing order can be explained satisfactorily using the classical nucleation theory. Indications of single-domain and multi-domain nucleation and of the corresponding changes in the nucleation rates are clearly observed. The nature of the second time scale of relaxation, apart from the field driven nucleation time, and the origin of its unusual large values at the phase boundary are explained from the disappearing tendency of kinks on the domain wall surfaces after the withdrawal of the pulse. The possibility of scaling behaviour in the multi-domain regime is identified and compared with the earlier observations.

Arkajyoti Misra111 e-mail : arko@cmp.saha.ernet.in  and Bikas K. Chakrabarti222 e-mail : bikas@cmp.saha.ernet.in

*Saha Institute of Nuclear Physics, 1/AF Bidhannagar, Calcutta 700 064, India.*

PACS. 64.60 Ht - Dynamic Critical Phenomena.

PACS. 64.60 Kw - Mulicritical Points.

PACS. 64.60.Qb - Nucleation.

The study of dynamical phase transitions in pure Ising systems, with ferromagnetic short range interactions, under the influence of time dependent external magnetic field has recently become one of the significant areas of interest in statistical physics . The effect of external magnetic fields which are periodic in time was first dealt with mean field theory . Subsequently, through extensive Monte Carlo studies, the existence of a dynamic phase transition was established and properly characterized . A relevant investigation in this context was the study of the effect of external pulsed field which is uniform in space but applied for a finite duration. All these pulsed field studies are concerned with the system below its static critical temperature $`T_c^0`$, where the system has a long range order characterized by a non-zero magnetization $`m_0`$. Under the influence of a ‘positive’ field or a field applied along the direction of prevalent order, the system does not show any new phase transition . However, a ‘negative’ pulsed field competes with the existing order and the system may show a transition from an initial equilibrium state with magnetization $`m_0`$ to a final equilibrium state with magnetization $`m_0`$ . Such a transition can be brought about by tuning the duration ($`\mathrm{\Delta }t`$) and strength ($`h_p`$) of the pulse, and the ‘phase diagram’ in the $`h_p\mathrm{\Delta }t`$ plane gives the minimal combination of these two parameters to bring about the transition at any temperature $`T`$ below $`T_c^0`$. The transition is dynamic in nature and both length and time scales are seen to diverge across this transition .

According to classical theory of nucleation, there could be two mechanisms through which the droplet of a particular spin grows in the sea of opposite spins. When the external magnetic field is relatively weak, a single droplet grow to cover the whole system size and this regime is called the single-droplet (SD)  or nucleation regime . On the other hand, under stronger magnetic fields, many small droplets grow simultaneously and eventually they coalesce to span the entire system. This is called the multi-droplet (MD) or coalescence regime. The crossover from the SD to the MD regime takes place at the dynamic spinodal field or $`H_{DSP}`$, which is a function of the system size $`L`$ and of the temperature $`T`$. There are two time scales in the problem : (i) the nucleation time $`\tau _N`$ is the time taken by the system to leave the metastable state under the influence of the external field and (ii) relaxation time $`\tau _R`$ which is defined as the time taken by the system to reach the final equilibrium state after the external field is withdrawn. In this letter, we show that the phase diagram of the magnetization-reversal transition can be explained satisfactorily by employing the classical nucleation theory . The nature of the transition changes form a discontinuous to a continuous one, giving rise to a ‘tricritical’ point as the system goes from the SD regime to the MD regime, depending on the temperature and the strength of the applied pulse. The dimension dependent factor by which the nucleation time differs in the two regimes is confirmed by Monte Carlo simulation of $`2d`$ kinetic Ising model evolving under Glauber dynamics below $`T_c^0`$. Our observations here also support the validity of the finite size scaling of the of the order parameter used earlier to obtain the dynamic critical exponents of the transition .

Monte Carlo simulation using single spin-flip Glauber dynamics have been used to study pure Ising system in 2$`d`$ below $`T_c^0`$ ($`2.27,\text{ }`$ in units of the strength of the cooperative interaction). The system is subjected to a magnetic field $`h(t)`$ of finite strength $`h_p`$ for a finite duration $`\mathrm{\Delta }t`$ :

$$\begin{array}{ccccc}h(t)& =& \hfill h_p& \text{for }t_0tt_0+\mathrm{\Delta }t\hfill & \\ & =& \hfill 0& \text{otherwise}.\hfill & \end{array}$$
(1)
In order to prepare the system in an ordered state corresponding to a temperature $`T<T_c^0`$, $`t_0`$ is taken much larger than the equilibrium relaxation time of the system at that temperature. The negative sign signifies that the applied field competes with the prevalent order characterized by equilibrium magnetization $`m_0`$, which is a function of $`T`$ only. Depending on the values of $`h_p`$ and $`\mathrm{\Delta }t`$, the system can either go back to the original state (with magnetization $`m_0`$) or to the other equivalent ordered state, characterized by equilibrium magnetization $`m_0`$, eventually after the withdrawal of the field. This gives rise to the magnetization-reversal transition which brings the system from one equilibrium ordered state to the other, below $`T_c^0`$.

It appears that the average magnetization at the time of withdrawal of the pulse $`m_wm(t_0+\mathrm{\Delta }t)`$ is a very relevant quantity : depending on the sign of $`m_w`$, the system chooses the final equilibrium state. On an average, the transition to the $`m_0`$ state occurs for negative values of $`m_w`$, whereas for positive $`m_w`$ the system goes back to the original $`+m_0`$ state. There occur of course fluctuations, where the transition to the opposite order takes place even for small positive $`m_w`$ values or vice versa. The phase boundary is therefore identified as the $`m_w=0`$ line. Fig. 1 shows the Monte Carlo phase boundaries at a few different temperatures below $`T_c^0`$ for a 100$`\times `$100 lattice with periodic boundary condition. At very low temperatures (typically for $`T<0.5`$), the transition is discontinuous in nature all along the phase boundary. Crossover to a continuous transition region appears along the phase boundary for higher temperatures, thereby giving rise to a tricritical point. Along a typical phase boundary where both kinds of transition are observed, the continuous transition region appears for smaller values of $`\mathrm{\Delta }t`$ (higher values of $`h_p`$), whereas the discontinuous region appears for higher values of $`\mathrm{\Delta }t`$ (smaller values of $`h_p`$).

It is instructive to look at the droplet picture of the classical nucleation theory to describe the nature of the phase diagram for the magnetization-reversal transition. The typical configuration of a ferromagnet below $`T_c^0`$ consists of clusters or droplets of down spins in a background of up spins or vice versa. According to the classical nucleation theory , the equilibrium number of droplets of size $`l`$ is then given by $`n_l=N\mathrm{exp}(\beta ϵ_l)`$, where $`\beta =1/k_BT`$ and $`ϵ_l`$ is the free energy of formation of a droplet of size $`l`$. Assuming spherical shape of the droplets, one can write $`ϵ_l=2hl+\sigma (T)l^{(d1)/d}`$, where $`\sigma `$ is the temperature dependent surface tension and $`h`$ is the external magnetic field. Droplets of size greater than a critical value $`l_c`$ are then favoured to grow. One obtains $`l_c`$ by maximizing $`ϵ_l`$ : $`l_c=\left[\sigma (d1)/2d\left|h\right|\right]^d`$. The number of supercritical droplets is then given by $`n_l^{}=N\mathrm{exp}\left(\beta K_d\sigma ^d/\left|h\right|^{d1}\right)`$, where $`K_d`$ is a constant depending on dimension only. In the nucleation regime or the SD regime, there is only one supercritical droplet that grows and engulfs the entire system. The nucleation time $`\tau _N^{SD}`$ is inversely proportional to the nucleation rate $`I`$. According to the Becker-D$`\ddot{\text{o}}`$ring theory , $`I`$ in turn is proportional to $`n_l^{}`$, and thus

$$\tau _N^{SD}I^1\mathrm{exp}\left(\frac{\beta K_d\sigma ^d}{\left|h\right|^{d1}}\right).$$

In the coalescence or MD regime, where due to stronger magnetic field many droplets grow simultaneously and eventually coalesce to form a system spanning droplet, the nucleation time is given by

$$\tau _N^{MD}I^{1/(d+1)}\mathrm{exp}\left(\frac{\beta K_d\sigma ^d}{(d+1)\left|h\right|^{d1}}\right).$$

During the time when the field is ‘on’, the only relevant time scale of the problem is the nucleation time. The phase boundary of the magnetization-reversal transition corresponds to the threshold value of the field pulse ($`h_p^c`$), which can bring the system from an equilibrium state with magnetization $`m_0`$ to a state with magnetization $`m_w=0_{}`$ in time $`\mathrm{\Delta }t`$, so that the system eventually evolves to a state with magnetization $`m_0`$. Equating therefore $`\mathrm{\Delta }t`$ with the nucleation time $`\tau _N`$, one gets for the magnetization-reversal phase diagram

$$\begin{array}{cccccc}\mathrm{ln}(\mathrm{\Delta }t)& =& c_1& +& \frac{C}{h_p^{d1}}& \text{in SD regime}\hfill \\ & =& c_2& +& \frac{C}{(d+1)h_p^{d1}}& \text{in MD regime},\hfill \end{array}$$
(2)
where $`C=\beta K_d\sigma ^d`$ and $`c_1`$, $`c_2`$ are constants. Therefore a plot of $`\mathrm{ln}(\mathrm{\Delta }t)`$ against $`\left(h_p^c\right)^{d1}`$ would show different slopes in the two regimes. This is indeed the case as shown in Fig. 2, where two distinct slopes are seen at higher temperatures (typically for $`T>0.5`$). The ratio of the slopes $`R`$ corresponding to the two regimes has a value close to $`3`$ as predicted for $`2d`$ by the classical nucleation theory. The point of intersection of the two straight lines gives the value of $`H_{DSP}`$ at that temperature and system size. At lower temperatures, however, the MD region is absent which is marked by a single slope (Fig. 2(d)).

According to the classical nucleation theory $`l_c\left|h_p\right|^d`$ and at any fixed temperature, therefore, one expects stronger (weaker) fields to bring the system to the MD (SD) regime. Fig. 3 shows snapshots of the spin configuration at a particular temperature, where the dots correspond to down spins. It is clear from the figure that many droplets of down spins are formed for a smaller value of $`\mathrm{\Delta }t`$ or equivalently for large $`h_p`$ whereas only a single down spin droplet is formed for a larger $`\mathrm{\Delta }t`$ or weaker $`h_p`$. The snapshots of the system are taken at the time of withdrawal of the pulse, beyond which the system is left to itself to relax back to either of the equilibrium states. The time taken by the system to reach the final equilibrium after the withdrawal of the pulse is defined as the relaxation time ($`\tau _R`$). It is observed  that $`\tau _R\kappa (T,L)\mathrm{exp}\left[\lambda (T)\left|m_w\right|\right]`$, where $`\kappa (T,L)`$ is a constant depending on temperature and system size ($`\kappa \mathrm{}\text{ as }L\mathrm{}`$) and $`\lambda (T)`$ is a constant depending on temperature. This is in distinct contrast with the normal relaxation of a ferromagnet to the equilibrium state with magnetization $`m_0`$, starting from a random initial state where the average magnetization is close to zero. The effect of the field is to initiate the nucleation process and by the time the pulse is withdrawn, the droplet(s) has (have) very few kinks along the surface(s). Once a droplet or a domain forms a flat boundary, it becomes a rather stable configuration and the domain wall movement practically stops; thereby restricting further nucleation. When the system is trapped in such a metastable state, it is left for large fluctuations to initiate further movement of the domain walls and resume the process of nucleation. The closer to zero is $`m_w`$, more is the chance that the system will be trapped in a metastable state owing to more number of flat domain walls. Thus the effect of the pulse is to create a domain structure which does not favour fast nucleation and even for $`TT_c^0`$, $`\tau _R`$ therefore diverges at the phase boundary.

To determine the order of the phase transition one can look at the probability distribution $`P(m_w)`$ of $`m_w`$, as one approaches a phase boundary. Fig. 4 shows the variation of $`P(m_w)`$ across the phase boundary corresponding to a particular $`T`$ at two different values of $`\mathrm{\Delta }t`$. The existence of a single peak of $`P(m_w)`$ in (a), which shifts its position continuously from $`+1`$ to $`1`$, indicates the continuous nature of the transition; whereas in (b), one obtains two peaks simultaneously, close to $`\pm m_0(T)`$, indicating that the system can simultaneously reside in two different phases with finite probabilities. This is a sure indication of discontinuous phase transition. The crossover from the discontinuous transition to a continuous one along the phase boundary is not very sharp; instead one gets a ‘tricritical region’ where the determination of the nature of the transition is not very conclusive. As is evident from Fig. 2, the data points in this particular region do not fit to either of the two straight lines corresponding to the two regimes. Thus for $`h_p^cH_{DSP}`$ the system belongs to the MD regime and the nature of the transition is continuous, whereas for $`h_p^cH_{DSP}`$ the system is brought to the SD regime where the transition is discontinuous in nature.

In our earlier work , we clearly noticed that finite size scaling of the order parameter was possible only at higher temperatures for moderately low $`\mathrm{\Delta }t`$. The reason behind that becomes clear now, as at higher temperatures the tricritical point shifts towards higher values of $`\mathrm{\Delta }t`$ and hence in the lower $`\mathrm{\Delta }t`$ region the system belongs to the MD regime and the transition is continuous in nature. This supports the scaling. Since at lower temperatures the continuous region gradually shrinks before disappearing altogether, all attempts for a finite size scaling fit failed. This observation also compares with that of Sides et al.  for the scaling behaviour in the case of dynamic transition under periodic fields.

In this letter we have shown that the phase diagram of the magnetization-reversal transition in the Ising model can be explained satisfactorily from the classical nucleation theory. Across $`H_{DSP}`$, the system goes from MD or coalescence regime where the transition is continuous in nature to SD or nucleation regime where the transition is discontinuous. A tricritical point separates these two regimes which moves towards smaller $`\mathrm{\Delta }t`$ (larger $`h_p`$) values as one decreases $`T`$ until it disappears altogether at low temperatures. The dimensional factor in the nucleation time of the two regimes is found to be quite accurately reproduced in the Monte Carlo simulations. There are two time scales in the problem, viz. $`\tau _N`$ and $`\tau _R`$. $`\tau _N`$, which is determined by the pulse duration $`\mathrm{\Delta }t`$, brings the system from one regime to the other depending on the temperature; while $`\tau _R`$ determines the relaxation of the field-released system and it diverges at the phase boundary even for $`TT_c^0`$. Finite size scaling of the order parameter for the transition is also justified in the MD regime where the transition is continuous in nature.

\***

A. M. would like to thank Muktish Acharyya and Burkhard Duenweg for useful discussions.

Figure Captions

Fig. 1. Monte Carlo phase diagram for magnetization-reversal transition in a $`100\times 100`$ lattice: (a) $`T=0.2`$, (b) $`T=0.6`$, (c) $`T=1.0`$, (d) $`T=1.4`$, (e) $`T=2.0`$.

Fig. 2. Plot of $`\mathrm{ln}(\mathrm{\Delta }t)`$ against $`\left[h_p^c(\mathrm{\Delta }t,T)\right]^1`$, obtained from the phase diagram in Fig. 1, at different temperatures. (a) $`T=1.4`$, (b) $`T=1.0`$, (c) $`T=0.7`$, (d) $`T=0.2`$. The observed values of the slope ratio $`R`$ are close to $`3.38`$, $`3.09`$ and $`3.24`$ in (a), (b) and (c) respectively.

Fig. 3. Snapshots of the spin configuration of the system at $`t=t_0+\mathrm{\Delta }t`$ for $`T=1.0`$; (a) $`\mathrm{\Delta }t=2`$, $`h_p=1.91`$ (MD regime) (b) $`\mathrm{\Delta }t=100`$, $`h_p=0.65`$ (SD regime).

Fig. 4. Probability distribution $`P(m_w)`$ of $`m_w`$ at $`T=0.6`$. In (a) where $`\mathrm{\Delta }t=1`$ (SD regime) , the position of the peak changes sharply from $`+1`$ to $`1`$ indicating a discontinuous transition: (i) $`h_p=2.05`$, (ii) $`h_p=2.15`$, (iii) $`h_p=2.25`$. In (b) where $`\mathrm{\Delta }t=500`$ (MD regime) , the position of the peak changes continuously from $`+1`$ to $`1`$ indicating a continuous transition: (i) $`h_p=0.92`$, (ii) $`h_p=1.02`$, (iii) $`h_p=1.12`$.


text:
# Averages of static electric and magnetic fields over a spherical region: a derivation based on the mean-value theorem

preprint: Submitted to Notes and Discussions in the American Journal of Phyiscs

The electromagnetic theory of dielectric and magnetic media deals with macroscopic electric and magnetic fields, because the microscopic details of these fields are usually experimentally irrelevant (with certain notable exceptions such as scanning tunneling microscopy). The macroscopic fields are the average of the microscopic fields over a microscopically large but macroscopically small region. This averaging region is often chosen to be spherical, denoted here as

$$\overline{𝐄}\frac{3}{4\pi R^3}_𝒱𝑑𝐫^{}𝐄(𝐫^{}),$$
(1)
where $`𝒱`$ is a sphere of radius $`R`$ centered at $`𝐫`$ (with a similar definition for $`\overline{𝐁}`$). $`R`$ is a distance which is macroscopically small but nonetheless large enough to enclose many atoms. The macroscopic $`\overline{𝐄}`$ and $`\overline{𝐁}`$ fields obtained by averaging over a sphere exhibit properties which prove useful in certain arguments and derivations. These properties are as follows:

1. if all the sources of the $`𝐄`$-field are outside the sphere, then $`\overline{𝐄}`$ is equal to the electric field at the center of the sphere,
2. if all the sources of the E-field are inside the sphere, $`\overline{𝐄}=𝐩/(4\pi ϵ_0R^3)`$, where $`𝐩`$ is the dipole moment of the sources with respect to the center of the sphere,
3. if all the sources of the $`𝐁`$-field are outside the sphere, then $`\overline{𝐁}`$ is equal to the magnetic field at the center of the sphere,
4. if all the sources of the B-field are inside the sphere, $`\overline{𝐁}=\mu _0𝐦/(2\pi R^3)`$, where $`𝐦`$ is the magnetic dipole moment of the sources.

These results can be derived in a variety of ways. For example, Griffiths derives properties 1 and 2 using a combination of results from Coulomb’s law and Gauss’ law, and properties 3 and 4 by writing down the $`𝐁`$-field in terms of the vector potential in the Coulomb gauge, and explicitly evaluating angular integrals. The purpose of this note is to describe a relatively simple derivation of all four results, based on the well-known mean-value theorem theorem (described in most textbooks on electromagnetic theory): if a scalar potential $`\mathrm{\Phi }(𝐫)`$ satisfies Laplace’s equation in a sphere, then the average of $`\mathrm{\Phi }`$ over the surface of the sphere is equal to $`\mathrm{\Phi }`$ at the center of the sphere; that is, if $`^2\mathrm{\Phi }=0`$ in a spherical region of radius $`r^{}`$ centered at $`𝐫`$, then

$$\mathrm{\Phi }(𝐫)=\frac{1}{4\pi }𝑑\mathrm{\Omega }\mathrm{\Phi }\left(𝐫+r^{}\widehat{𝐧}_\mathrm{\Omega }\right),$$
(2)
where $`\mathrm{\Omega }`$ is the solid angle relative to $`𝐫`$ and $`\widehat{𝐧}_\mathrm{\Omega }`$ is the unit vector pointing in the direction of $`\mathrm{\Omega }`$.

Taking the gradient $``$ of both sides of Eq. (2) with respect to r, we immediately obtain

$$𝐄(𝐫)=\mathrm{\Phi }(𝐫)=\frac{1}{4\pi }𝑑\mathrm{\Omega }\left[\mathrm{\Phi }\left(𝐫+r^{}\widehat{𝐧}_\mathrm{\Omega }\right)\right]=\frac{1}{4\pi }𝑑\mathrm{\Omega }𝐄\left(𝐫+r^{}\widehat{𝐧}_\mathrm{\Omega }\right);$$
(3)
that is, if $`𝐄=0`$ inside a sphere, the average of $`𝐄`$ over the surface of the sphere is equal to the $`𝐄`$ at the center of sphere. Eq. (3) is the basis of the derivation of the four properties listed above. To simplify the notation, henceforth in this note it is assumed that $`𝒱`$ is a sphere of radius $`R`$ and is centered at the origin $`\mathrm{𝟎}`$.

1. Electric field with sources outside sphere

When all the static charge sources are outside $`𝒱`$, $`𝐄=\mathrm{\Phi }`$ and $`𝐄=0`$ inside $`𝒱`$ and hence Eq. (3) is valid. We shall now see that property 1 follows quite trivially as a special case of Eq. (3) \[and hence of Eq. (2)\].

The average of $`𝐄`$ over $`𝒱`$ can be written as weighted integral of the average of $`𝐄`$ over surfaces of spheres with radius $`r^{}<R`$ centered at $`\mathrm{𝟎}`$,

$`\overline{𝐄}{\displaystyle \frac{3}{4\pi R^3}}{\displaystyle _𝒱}𝑑𝐫^{}𝐄(𝐫^{})`$ $`=`$ $`{\displaystyle \frac{3}{R^3}}{\displaystyle _0^R}r^2𝑑r^{}\left[{\displaystyle \frac{1}{4\pi }}{\displaystyle 𝑑\mathrm{\Omega }𝐄(r^{}\widehat{𝐧}_\mathrm{\Omega })}\right].`$ (4)
Using Eq. (3) in Eq. (4) immediately yields property 1,

$$\overline{𝐄}=\frac{3}{R^3}_0^Rr^2𝑑r^{}𝐄(\mathrm{𝟎})=𝐄(\mathrm{𝟎}).$$
(5)
2. Electric field with sources inside sphere

For clarity we first prove property 2 for a single point charge inside the sphere. The general result can be inferred from the single point charge result using the superposition principle, but for completeness we generalize the proof for a continuous charge distribution.

Utilizing the vector identity

$$_𝒱\mathrm{\Phi }d𝐫=_𝒮\mathrm{\Phi }𝑑𝐚,$$
(6)
where $`𝒮`$ is the surface of $`𝒱`$, the average of the electric field over $`𝒱`$ centered at $`\mathrm{𝟎}`$ can be written as

$$\overline{𝐄}=\frac{3}{4\pi R^3}_𝒱\mathrm{\Phi }(𝐫)𝑑𝐫=\frac{3}{4\pi R^3}_S\mathrm{\Phi }𝑑𝐚.$$
(7)
Eq. (7) implies the average $`\overline{𝐄}`$ is determined completely by the potential on the surface $`𝒮`$.

We now use a well-known result from the method of images solution for the potential of a point charge next to a conducting sphere: the potential on the surface $`S`$ for a charge $`q`$ at $`𝐫`$ inside the sphere is reproduced exactly by an image charge $`q^{}=Rq/d`$ at $`𝐫^{}=(R^2/r)\widehat{𝐫}`$ outside the sphere ($`\widehat{𝐫}`$ is a unit vector). Eq. (7) therefore implies that the $`\overline{𝐄}`$ for a point charge $`q`$ at $`𝐫`$ is exactly equal to that of an image point charge $`q^{}`$ at $`(R^2/r)\widehat{𝐫}`$. But since the image charge is outside the sphere, we can use property 1 to determine $`\overline{𝐄}`$,

$$\overline{𝐄}=𝐄_{\mathrm{image}}(\mathrm{𝟎})=\frac{q^{}\widehat{𝐫}}{4\pi ϵ_0|𝐫^{}|^2}=\frac{q𝐫}{4\pi ϵ_0R^3}=\frac{𝐩}{4\pi ϵ_0R^3},$$
(8)
where $`𝐩=q𝐫`$ is the dipole moment for a single point charge.

Generalization to continuous charge distributions – Assume the charge distribution inside the spherical volume $`𝒱`$ is $`\rho (𝐫)`$. A volume element $`dV`$ at $`𝐫`$ inside $`𝒱`$ contains charge $`dq=\rho (𝐫)dV`$. The image charge element outside the sphere which gives the same average electric field as $`dq`$ is $`dq^{}=\rho ^{}(𝐫^{})dV^{}=dqR/r`$ at $`𝐫^{}=(R^2/r)\widehat{𝐫}`$. As in the discrete case, the contribution of $`dq`$ to the average $`𝐄`$-field in $`𝒱`$ is equal to the electric field at the origin due to $`dq^{}`$,

$$d\overline{𝐄}=d𝐄_{\mathrm{image}}(\mathrm{𝟎})=\frac{\rho ^{}(𝐫^{})dV^{}}{4\pi ϵ_0r^2}\widehat{𝐫}=\frac{𝐫\rho (𝐫)dV}{4\pi ϵ_0R^3}.$$
(9)
The averaged electric field due to all charges in the sphere $`𝒱`$ is therefore

$$\overline{𝐄}=_𝒱𝑑\overline{𝐄}=\frac{1}{4\pi ϵ_0R^3}_𝒱𝑑V𝐫\rho (𝐫)=\frac{𝐩}{4\pi ϵ_0R^3},$$
(10)
where $`𝐩=_𝒱𝑑V𝐫\rho (𝐫)`$ is the dipole moment with respect to the origin of all the charges in $`𝒱`$.

3. Magnetic field with sources outside the sphere

When magnetic field sources are absent in $`𝒱`$, both $`𝐁=0`$ and $`\times 𝐁=0`$, and hence $`𝐁=\varphi _M`$ where $`^2\varphi _M=0`$ inside the sphere. Therefore, the same derivation for property 1 holds here; that is, the average of the $`𝐁`$-field over a sphere is equal to its value at the center of the sphere.

4. Magnetic field for sources inside the sphere

Using the vector potential description of the magnetic field, $`𝐁=\times 𝐀`$, the average over the sphere $`𝒱`$ can be written as

$$\overline{𝐁}=\frac{3}{4\pi R^3}_𝒱\times 𝐀=\frac{3}{4\pi R^3}_𝒮𝐀\times 𝑑𝐚.$$
(11)
The second equality in the above equation is a vector identity. This equation shows that, as in the case of the $`𝐄`$-field and the scalar potential, the average $`𝐁`$-field over any volume $`𝒱`$ can be computed from $`𝐀`$ on the surface of $`𝒱`$.

We now consider the contribution to $`\overline{𝐁}`$ of current element $`𝐉(𝐫)dV`$ inside $`𝒱`$. We can do this by determining the image current element outside the sphere that exactly reproduces the vector potential due to $`𝐉(𝐫)dV`$ on the surface of the sphere. We choose the Coulomb gauge

$$A_i(𝐫)=\frac{\mu _0}{4\pi }\frac{J_i(𝐱)}{|𝐫𝐱|}𝑑𝐱,$$
(12)
where $`i`$ denotes spatial component. Since $`A_i`$ is related to $`J_i`$ in the same way as the electric potential $`\mathrm{\Phi }`$ is to the charge density $`\rho `$, the method of electrostatic images can also be used here to determine the image current element.

The proof of property 4 proceeds analogously to that of property 2. For $`𝐉(𝐫)dV`$ inside the sphere $`𝒱`$, $`𝐀`$ on the surface $`𝒮`$ is reproduced by an image current element $`𝐉^{}(𝐫^{})dV^{}=(R/r)𝐉(𝐫)dV`$, where $`𝐫^{}=(R^2/r)\widehat{𝐫}`$. Since the image current is outside the sphere, we can use property 3. Thus, the $`𝐉(𝐫)dV`$ contribution to the average $`𝐁`$-field in $`𝒱`$ is equal to the $`𝐁`$-field at the center due to $`𝐉^{}(𝐫^{})dV^{},`$

$$d\overline{𝐁}=d𝐁_{\mathrm{image}}(\mathrm{𝟎})=\frac{\mu _0}{4\pi }\frac{𝐉^{}(\mathrm{r}^{})\times \widehat{𝐫}}{r^2}dV^{}=\frac{\mu _0}{4\pi R^3}𝐫\times 𝐉(𝐫)dV.$$
(13)
The contribution for the entire current distribution in $`𝒱`$ is therefore

$$\overline{𝐁}=\frac{\mu _0}{4\pi R^3}_𝒱𝑑V𝐫\times 𝐉(𝐫)=\frac{\mu _0𝐦}{2\pi R^3},$$
(14)
where $`𝐦=\frac{1}{2}_𝒱𝑑V𝐫\times 𝐉(𝐫)`$ is magnetic moment of a current distribution in $`𝒱`$.

Finally, note that similar arguments hold for charge distributions which are constant along the $`z`$-direction, since the potentials for these satisfy Laplace’s equation in two dimensions, and the method of images is also applicable for line charges.

Acknowledgement

Useful correspondence with Prof. David Griffiths is gratefully acknowledged.


text:
# Conformal Invariance and Particle Aspects in General Relativity

## 1 Introduction

Conformal invariance has been playing a particularly important role in the investigation of gravitational models ever since the emergence of such theories. In a system which includes matter, it is well known that conformal invariance requires the vanishing of the trace of the stress tensor in the absence of dimensional parameters. In the presence of dimensional parameters, the conformal invariance can be established for a large class of theories  if the dimensional parameters are conformally transformed according to their dimensions. One general feature of conformally invariant theories is, therefore, the presence of varying dimensional coupling constants. In particular, one can say that the introduction of a constant dimensional parameter into a conformally invariant theory breaks the conformal invariance in the sense that a preferred conformal frame is singled out, namely that in which the dimensional parameter has the assumed constant configuration. Thus the breakdown of conformal invariance may be established by introducing a constant dimensional parameter into the theory. The determination of the corresponding preferred conformal frame depends on the nature of the problem at hand. In a conformally invariant gravitational model one usually considers the symmetry breaking as a cosmological effect. This would mean that one breaks the conformal symmetry by defining a preferred conformal frame in terms of the large scale properties of cosmic matter distributed in a finite universe. In this way, the breakdown of conformal symmetry was found to be a framework in which one can look for the origin of the gravitational coupling of matter, both classical  and quantum , at large cosmological scales.

The purpose of this paper is to show that the cosmological breakdown of conformal invariance in a conformally invariant gravitational model together with a local change of the corresponding preferred conformal frame may be used to model a particle concept in general relativity.

The organisation of the paper is as follows: We first study the breakdown of conformal symmetry in a conformally invariant gravitational model and define the resulting preferred conformal frame in terms of some cosmological characteristics in close correspondence to the work in . We then show that, by a local change of the preferred conformal frame, it is possible to derive a Hamilton-Jacobi type equation with adjustable mass. In our presentation there is a substantial dynamical interplay between a particle in the ensemble and the applied conformal factor in a form which is similar to the dynamical effect of the quantum potential on a particle in the context of the causal interpretation of relativistic quantum mechanics . This result seems to be interesting because it suggests that the emergence of quantal behaviour of matter may be a general feature of a conformal invariant gravitational model. We shall work with a metric having the signature (- + + +).

## 2 Breakdown of conformal invariance

In this section we briefly review the work in . Consider the action functional

$$S[\varphi ]=\frac{1}{2}d^4x\sqrt{g}(g^{\mu \nu }_\mu \varphi _\nu \varphi +\frac{1}{6}R\varphi ^2)$$
(1)
which describes a system consisting of a real scalar field $`\varphi `$ non-minimally coupled to gravity, $`R`$ is the scalar curvature. Variations with respect to $`\varphi `$ and $`g_{\mu \nu }`$ lead to the equations

$$(\mathrm{}\frac{1}{6}R)\varphi =0$$
(2)
$$G_{\mu \nu }=6\varphi ^2\tau _{\mu \nu }(\varphi )$$
(3)
where $`G_{\mu \nu }=R_{\mu \nu }\frac{1}{2}g_{\mu \nu }R`$ is the Einstein tensor and

$$\tau _{\mu \nu }(\varphi )=[_\mu \varphi _\nu \varphi \frac{1}{2}g_{\mu \nu }_\alpha \varphi ^\alpha \varphi ]\frac{1}{6}(g_{\mu \nu }\mathrm{}_\mu _\nu )\varphi ^2$$
(4)
with $`_\mu `$ denoting the covariant derivative. Taking the trace of (3) gives

$$\varphi (\mathrm{}\frac{1}{6}R)\varphi =0$$
(5)
which is consistent with equation (2). This is a consequence of the conformal symmetry of action (1) under the conformal transformations

$$\varphi \overline{\varphi }=\mathrm{\Omega }^1(x)\varphi ,g_{\mu \nu }\overline{g}_{\mu \nu }=\mathrm{\Omega }^2(x)g_{\mu \nu }$$
(6)
where the conformal factor $`\mathrm{\Omega }(x)`$ is an arbitrary, positive and smooth function of space-time. Adding a matter source $`S_m`$ independent of $`\varphi `$ to the action (1) in the form

$$S=S[\varphi ]+S_m$$
(7)
yields the field equations

$$(\mathrm{}\frac{1}{6}R)\varphi =0$$
(8)
$$G_{\mu \nu }=6\varphi ^2[\tau _{\mu \nu }(\varphi )+T_{\mu \nu }]$$
(9)
where $`T_{\mu \nu }`$ is the matter energy-momentum tensor. The following algebraic requirement

$$T_\mu ^\mu =0$$
(10)
then emerges as a consequence of comparing the trace of (9) with (8). This implies that only traceless matter can couple consistently to such gravity models.

We may break the conformal symmetry by adding a dimensional mass term $`\frac{1}{2}d^4x\sqrt{g}\mu ^2\varphi ^2`$, with $`\mu `$ being a constant mass parameter, to the action (7). This leads to

$$\mu ^2\varphi ^2=T_\mu ^\mu .$$
(11)
Consequently, the field equations become

$$(\mathrm{}\frac{1}{6}R\mu ^2)\varphi =0$$
(12)
$$G_{\mu \nu }3\mu ^2g_{\mu \nu }=6\varphi ^2[\tau _{\mu \nu }(\varphi )+T_{\mu \nu }].$$
(13)
Now the basic input is to consider the invariance breaking as a cosmological effect. This would mean that one may take $`\mu ^1`$ as the length scale characterizing the typical size of the universe $`R_0`$ and $`T_\mu ^\mu `$ as the average density of the large scale distribution of matter $`MR_0^3`$, where $`M`$ is the mass of the universe. This leads, as a consequence of (11) to the estimation of the constant background value of $`\varphi `$

$$\varphi ^2R_0^2(M/R_0^3)^1R_0/MG$$
(14)
where the well-known empirical cosmological relation $`GM/R_01`$ has been used. Inserting this background value of $`\varphi `$ into the field equations (12), (13) leads to the following set of Einstein equations

$$G_{\mu \nu }3\mu ^2g_{\mu \nu }=6\varphi ^2T_{\mu \nu }GT_{\mu \nu }$$
(15)
with a correct coupling constant $`8\pi G`$, and a term $`3\mu ^2`$ which appears as an effective cosmological constant $`\mathrm{\Lambda }`$ of the order of $`R_0^2`$. The field equation (12) for $`\varphi `$ contains no new information.

We should emphasize that the invariance breaking explained above dealt with a broken conformal invariance in the presence of a dimensional matter source which effectively appeared as a cosmological constant $`\mathrm{\Lambda }`$. By implication, we get a preferred conformal frame $`(\varphi ,g_{\mu \nu },\mu )`$ for the gravitational variables, namely that in which $`\varphi ^2G^1`$, $`\mu ^2\mathrm{\Lambda }`$ and $`g_{\mu \nu }`$ determined by the field equations (15). This preferred confromal frame has the remarkable property that it incorporates a correct coupling of the cosmic matter to gravity. We shall call it the cosmological frame.

## 3 Particle interpretation

A key feature of any fundamental theory consistent with a given symmetry is that its breakdown would lead to effects which can have various manifestations of physical importance. Therefore, in the case of conformal symmetry, one would expect that the corresponding cosmological invariance breaking, considered in the last section, would have an important effect on the local structure of the underlying theory. Here we would like to study one particular effect which illustrates the local particle aspects.

There are two chracteristic scales of squred mass in the cosmological frame $`(\varphi ,g_{\mu \nu },\mu )`$ which are subjected to the scale hierarchy 111We use units in which $`\mathrm{}=c=1`$. In these units $`\varphi `$ has the dimension of mass.

$`\mu ^2\varphi ^2G^1.`$ (16)
In terms if this scale hierarchy, it is possible to single out a congruence of hypersurface orthogonal timelike curves, if we subject the corresponding generating vector field $`_\mu S`$ to the condition

$`_\mu S^\mu S=(\varphi ^2\mu ^2).`$ (17)
We can observe that no known elementary particle can move along the congruence defined by (17) because of the large mass-scale defined by the right hand side of (17). Actually, we know from (16) that

$$\varphi ^2\mu ^2\varphi ^2G^1=m_p^2$$

where $`m_p`$ is the Plank mass . The situation, however, changes if we consider a local change of the cosmological frame. To illustrate this point, let us write equation (17) in a new frame $`(\overline{\varphi },\overline{g}_{\mu \nu },\overline{\mu })`$ which is locally connected to the frame $`(\varphi ,g_{\mu \nu },\mu )`$ by a conformal transformation (6). Taking into account that $`S`$ as a dimensionless quantity does not transform under conformal invariance, we find

$`\overline{}_\mu S\overline{}^\mu S=(\overline{\varphi }^2\overline{\mu }^2)`$ (18)
with $`\overline{\varphi }=\mathrm{\Omega }^1\varphi `$, $`\overline{\mu }=\mathrm{\Omega }^1\mu `$ and the bar quantities refering to the new frame. Now, from the field equation (12) we find

$$\overline{\mu }^2=\frac{\stackrel{}{\mathrm{}}\overline{\varphi }}{\overline{\varphi }}\frac{1}{6}\overline{R}$$

which, if combined with equation (17) leads to

$`\overline{}_\mu S\overline{}^\mu S=\overline{\varphi }^2+{\displaystyle \frac{\stackrel{}{\mathrm{}}\overline{\varphi }}{\overline{\varphi }}}{\displaystyle \frac{1}{6}}\overline{R}.`$ (19)
Now, if the new conformal frame $`(\overline{\varphi },\overline{g}_{\mu \nu },\overline{\mu })`$ is taken to be subjected to

$`\overline{}_\mu S\overline{}^\mu \overline{\varphi }=0`$ (20)
then equation (19) can be interpreted as a Hamilton-Jacobi equation for a particle with adjustable mass-scale $`\overline{\varphi }^2`$. The relation (20) ensures that this mass-scale will not change along the particle trajectory. In addition, we get a dynamical effect on the particle trajectories, reflected in the term $`\frac{\stackrel{}{\mathrm{}}\overline{\varphi }}{\overline{\varphi }}\frac{1}{6}\overline{R}`$, which illustrates a modification of the particle mass due to the scalar curvature and the applied conformal transformation. In summary, we observe that starting from the given preferred cosmological frame $`(\varphi ,g_{\mu \nu },\mu )`$, we may derive a particle concept in terms of a corresponding local change of the conformal frame. In this approach no separate particle action needs to be introduced, that of a conformally invariant gravitational model suffices. The particle aspects emerge from an internal condition connecting the local properties of a time-like congruence of curves associated with a characteristic scale hierarchy in the cosmological frame with particle properties in a new conformal frame. This observation emphasizes that general relativity, if suitably formulated as a conformal invariant field theory, does not ascribe any special significance to a separate particle action. We should, however, note that there is a certain limitation on the applicability of the particle concept. Actually, in a universe in which $`\varphi ^2`$ is smaller than $`\mu ^2`$, the above particle concept does not apply, for the mass scale of (17) becomes tachionic. This limitation, however, does not seem to be a weakness of our particle concept because the condition $`\varphi ^2<\mu ^2`$ describes a universe of trans-Plankian size.

## 4 Concluding remarks

In this paper we have shown that the cosmological breakdown of conformal symmetry in a conformally invariant gravitational model together with a local change of the corresponding preferred conformal frame leads to a picture consistent with a particle concept. This picture may be considered as a manifestation of Mach’s principle in that the particle concept emerges as a local effect emerging from large scale cosmological consideration.
We emphasize the similarity of the term $`\frac{\stackrel{}{\mathrm{}}\overline{\varphi }}{\overline{\varphi }}`$ on the right hand side of (19) with the quantum potential term in the context of the causal interpretation of relativistic quantum mechanics . This similarity merits attention because it provides an indication for a possible geometrization of quantal behaviour of relativistic particles in the framework of a conformal invariant gravitational model.


text:
# Formalizing the gene centered view of evolution

## Abstract

A historical dispute in the conceptual underpinnings of evolution is the validity of the gene centered view of evolution.  We transcend this debate by formalizing the gene centered view as a dynamic version of the mean field approximation. This establishes the conditions under which it is applicable and when it is not. In particular, it breaks down for trait divergence which corresponds to symmetry breaking in evolving populations.

The gene centered view addresses a basic problem in the interplay of selection and heredity in sexually reproducing organisms. Sexual reproduction disrupts the simplest view of evolution because the offspring of an organism are often as different from the parent as organisms that it is competing against. In the gene centered view the genes serve as indivisible units that are preserved from generation to generation. In effect, different versions of the gene, i.e. alleles, compete rather than organisms. It is helpful to explain this using the “rowers analogy” introduced by Dawkins.  In this analogy boats of mixed left- and right-handed rowers are filled from a common rower pool. Boats compete in heats and it is assumed that a speed advantage exists for boats with more same-handed rowers. The successful rowers are then reassigned to the rower pool for the next round. Over time, a predominantly and then totally single handed rower pool will result. Thus, the selection of boats serves, in effect, to select rowers who therefore may be considered to be competing against each other. Other perspectives on evolution distinguish between vehicles of selection (the organisms) and replicators (the genes). However, a direct analysis of the gene centered view to reveal its domain of applicability has not yet been discussed. The analysis provided here, including all the equations, is applicable quite generally, but for simplicity it will be explained in terms of the rowers analogy.111The rowers analogy may be considered a model of a single gene in an $`n`$-ploid organism with $`n`$ the number of rowers, or a model of $`n`$ genes with two alleles per gene and each pair labeled correspondingly. The formal discussion applies to complete genomes i.e. to homolog and non-homolog genes.

The formal question is: Under what conditions (assumptions) can allelic (rower) competition serve as a surrogate for organism (boat) competition in the simple view of evolution. Formalizing this question requires identifying the conditions attributed to two steps in models of evolution, selection and reproduction. In the selection step, organisms (boats) are selected, while in the sexual reproduction step, new organisms are formed from the organisms that were selected. This is not fully discussed in the rowers model, but is implicit in the statement that victorious rowers are returned to the rower pool to be placed into new teams. The two steps of reproduction and selection can be written quite generally as:

$`\{N(s;t)\}`$ $`=`$ $`R[\{N^{}(s;t1)\}]`$ (1)
$`\{N^{}(s;t)\}`$ $`=`$ $`D[\{N(s;t)\}]`$ (2)
The first equation describes reproduction. The number of offspring $`N(s;t)`$ having a particular genome $`s`$ is written as a function of the reproducing organisms $`N^{}(s;t1)`$ from the previous generation. The second equation describes selection. The reproducing population $`N^{}(s;t)`$ is written as a function of the same generation at birth $`N(s;t)`$. The brackets on the left indicate that each equation represents a set of equations for each value of the genome. The brackets within the functions indicate, for example, that each of the offspring populations depends on the entire parent population.

To formalize the gene centered view, we introduce a dynamic form of what is known in physics as the mean field approximation. In the mean field approximation the probability of appearance of a particular state of the system (i.e. a particular genome, $`s`$) is the product of probabilities of the components (i.e. genes, $`s_i`$)

$$P(s_1,\mathrm{},s_N)=P(s_i).$$
(3)
In the usual application of this approximation, it can be shown to be equivalent to allowing each of the components to be placed in an environment which is an average over the possible environments formed by the other components of the system, hence the term “mean field approximation.” The key to applying this in the context of evolution is to consider carefully the effect of the reproduction step, not just the selection step.

In many models of evolution that are discussed in the literature, the offspring are constructed by random selection of surviving alleles (a panmictic population). In the rowers analogy the return of successful rowers to a common pool is the same approximation. This approximation eliminates correlations in the genome that result from the selection step and thus imposes Eq. (3), the mean field approximation, on the reproduction step for the alleles of offspring. Even though it is not imposed on the selection step, inserting this approximation into the two step process allows us to write both of the equations in Eq. (4) together as an effective one-step update

$$P^{}(s_i;t)=\stackrel{~}{D}[\{P^{}(s_i;t1)\}]$$
(4)
which describes the allele population change from one generation to the next of offspring at birth. Since this equation describes the behavior of a single allele it corresponds to the gene centered view.

There is still a difficulty pointed out by Lewontin.  The effective fitness of each allele depends on the distribution of alleles in the population. Thus, the fitness of an allele is coupled to the evolution of other alleles. This is apparent in Eq. (4) which, as indicated by the brackets, is a function of all the allele populations. It corresponds, as in other mean field approximations, to placing an allele in an average environment formed from the other alleles. For example, there is a difference of likelihood of victory (fitness) between a right-handed rower in a predominantly left-handed population, compared to a right-handed rower in a predominantly right-handed population. Since the population changes over time, fitnesses are time dependent and therefore not uniquely defined. This problem with fitness assignment would not be present if each allele separately coded for an organism trait. While this is a partial violation of the simplest conceptual view of evolution, however, the applicability of a gene centered view can still be justified, as long as the contextual assignment of fitness is included. When the fitness of organism phenotype is dependent on the relative frequency of phenotypes in a population of organisms it is known as frequency dependent selection, which is a concept that is being applied to genes in this context.

A more serious breakdown of the mean field approximation arises from what is known in physics as symmetry breaking. This corresponds in evolution to trait divergence of subpopulations. Such trait divergence arises when correlations in reproduction exist so that reproduction does not force complete mixing of alleles. The correlations in reproduction do not have to be trait related. For example, they can be due to spatial separation of organisms causing correlations in reproduction among nearby organisms. Models of spatially distributed organisms are sometimes called models of spatially structured environments. However, this terminology suggests that the environment itself is spatially varying and it is important to emphasize that symmetry breaking / trait divergence can occur in environments that are uniform (hence the terminology “symmetry breaking”). In the rowers model this has direct meaning in terms of the appearance of clusters of mostly left and mostly right handed rowers if they are not completely mixed when reintroduced and taken from the rower pool. Trait related correlations in sexual reproduction (assortive mating) caused by, e.g. sexual selection, would also have similar consequences. In either case, the gene centered view would not apply.

Historically, the gene-centered view of evolution has been part of the discussion of attitudes toward altruism and group selection and related socio-political as well as biological concerns.  Our focus here is on the mathematical applicability of the gene-centered view in different circumstances. While the formal discussion we present here may contribute to the socio-political issues, we have chosen to focus here on mathematical concerns.

The problem of understanding the mean-field approximation in application to biology is, however, also relevant to the problem of group selection. In typical models of group selection asexually (clonally) reproducing organisms have fecundities determined both by individual traits and group composition. The groups are assumed to be well defined, but periodically mixed. Similar to the gene-centered model, an assumption of random mixing is equivalent to a mean field theory. Sober and Wilson (1998) have used the term “the averaging fallacy” to refer to the direct assignment of fitnesses to individuals. This captures the essential concept of the mean-field approximation. However, both the limitations of this approximation in some circumstances and its usefulness in others do not appear to be generally recognized. For example, it is not necessary for well defined groups to exist for a breakdown in the mean-field approximation to occur. Correlations in organism influences are sufficient. Moreover, standard group-selection models rely upon averaging across groups with the same composition. For this case, where well defined groups exist and correlations in mixing satisfy averaging (mean-field) assumptions by group composition, equations developed by Price222See discussion on pp. 73–74 of . separate and identify both the mean field contribution to fitness and corrections due to correlations. These equations do not apply in more general circumstances when correlations exist in a network of interactions and/or groups are not well defined, and/or averaging across groups does not apply. It is also helpful to make a distinction between the kind of objection raised by Lewontin to the use of averaging, and the failure that occurs due to correlations when the mean-field approximation does not apply. In the former case, the assignment of fitnesses can be performed through the effect of the environment influencing the gene, in the latter case, an attempt to assign fitnesses to a gene would correspond to inventing non-causal interactions between genes.

The mean field approximation is widely used in statistical physics as a “zeroth” order approximation to understanding the properties of systems. There are many cases where it provides important insight to some aspects of a system (e.g. the Ising model of magnets) and others where is essentially valid (conventional BCS superconductivity). The application of the mean field approximation to a problem involves assuming an element (or small part of the system) can be treated in the average environment that it finds throughout the system. This is equivalent to assuming that the probability distribution of the states of the elements factor. Systematic strategies for improving the study of systems beyond the mean field approximation both analytically and through simulations allow the inclusions of correlations between element behavior. An introduction to the mean-field approximation and a variety of applications can be found in Bar-Yam (1997).

In conclusion, the gene centered view can be applied directly in populations where sexual reproduction causes complete allelic mixing, and only so long as effective fitnesses are understood to be relative to the prevailing gene pool. However, structured populations (i.e. species with demes—local mating neighborhoods) are unlikely to conform to the mean field approximation / gene centered view. Moreover, it does not apply to considering the consequences of trait divergence, which can occur when such correlations in organism mating occur. These issues are important in understanding problems that lie at scales traditionaly between the problems of population biology and those of evolutionary theory: e.g. the understanding of ecological diversity and sympatric speciation.