text:
# Untitled Document

An Analysis of the Thermodynamics of Hydrophobic Solvation Based on Scaled Particle Theory

Hong Qian

Department of Applied Mathematics and Bioengineering University of Washington, Seattle, WA 98195

Synopsis

A comprehensive, semi-quantitative model for the thermodynamics of hydrophobic solvation is presented. The model is based on a very simple premise suggested by the scaled particle theory and treats both solute and solvent molecules as hard spheres. A connection between the peculiarly large heat-capacity change for hydrophobic solvation and the large temperature dependence of the thermal expansivity of water is found. Analysis reveals a possible physical origin for the converging behavior of solvation entropies for a series of homologous hydrophobic compounds. The model suggests that the low solubility and the large heat-capacity change of hydrophobic solvation stem from two distinct aspects of water molecules: the static geometry of the molecule and the dynamic hydrogen bonding network, respectively.

Running title: Hydrophobic Solvation Thermodynamics

Introduction

The scaled particle theory (SPT), a statistical mechanical model for liquids and solutions made of hard spherical molecules, is a powerful conceptual framework and a useful computational tool for molecular solvation.$`^{14}`$ Even though the SPT was originally developed only for hard spherical molecules,$`^5`$ recent studies have shown that it also provides insight into the solvation of hydrophobic solutes in associative liquid like water. There are several reasons for this wide applicability of SPT:

1) The SPT is not a statistical mechanical theory based on the first principle. Rather, an independently obtained solvent density has to be provided in addition to temperature, pressure, and molecular parameters. This is in contrast to a genuine statistical mechanical theory of liquid, in which density is derived from temperature, pressure, and molecular parameters. Some essential properties of an associative liquid are contained in the density data.$`^{5,6}`$

2) Recent studies have shown that the SPT provides accurate numerical calculations for solvation free energy of hydrocarbons in water, but does so with less success for calculating the entropy and enthalpy.$`^{2,6}`$ This result is consistent with the concept of entropy-enthalpy compensation, which occurs with the solvent reorganization when introducing a solute molecule into a solvent.$`^{69}`$

3) To apply the SPT to an associative liquid like water, one practical difficulty is how to determine the hard sphere radii for water and the solute molecules.$`^2`$ While these difficulties affect quantitative numerical calculation, it does not alter the qualitative physical insight one obtains from the SPT. In this paper, our primary interest is in the qualitative relationship between different physical quantities, the numerical values are not essential for our analysis.

Based on these arguments, we try to use SPT as a semi-quantitative model for understanding various experimental results which are essential to the solvation of small organic molecules in aqueous solutions. We show that much of the characteristics of hydrophobic solvation can be understood in terms of this over-simplified model. Our approach is very much in the same spirit of the earlier work of Grunwald.$`^8`$ Treating the hydrophobic solutes as inert hard spheres, of course, completely neglects the solute-solvent interaction.$`^{3,6}`$ This should be kept in mind when applying the result of the present analysis to the solvation of hydrocarbons. On the other hand, the dissolution of inert gases like xenon exhibits all the characteristics of hydrophobic solvation. Hence, neglecting solute-solvent interaction is justified in our present study.

Some Basic Facts

The essential features in the data for dissolution of hydrophobic compounds in aqueous solutions are summarized as follows:

1) Dissolution entropy changes, $`\mathrm{\Delta }S(T)`$’s, associated with transferring a series of hydrophobic solutes from solid $``$ water, liquid $``$ water, and gas $``$ water form three distinct, respective groups. Within each group, $`\mathrm{\Delta }S(T)`$’s for different solute species converge to a common $`\mathrm{\Delta }S^{}`$ at a temperature $`T_s^{}`$. While $`\mathrm{\Delta }S^{}`$’s are different for the three groups, $`T_s^{}`$’s for all three groups are approximately equal to $`110^{}`$C.$`^{10}`$ The two groups of data on the transfer from liquid $``$ water and gas $``$ water are in agreement with the Trouton’s rule.

In order to compare free energies calculated from SPT with experimentally obtained solubility, one has to calculate the solvation free energy from a measurement based on concentration scale in molarity rather than in mole fraction.$`^{5,11,12}`$ The numerical difference between the two calculations is $`RTln(v_1/v_2)`$ where $`v_1`$ and $`v_2`$ are the molar volumes of the solvent systems before and after the transfer, respectively. For the gas phase, $`v=k_BT/p`$ where $`p`$ is the gas pressure (= 1 atm); for water, $`v`$ =18.08 $`cm^3/mol`$.$`^6`$ This changing of concentration scale increases the free energy values given by Murphy et al.$`^{10}`$ by about 60 J per degree per mole$`^{13}`$ (e.g., shift the lines in their Figure 1 upward).

2) The heat-capacity changes, $`\mathrm{\Delta }C_p`$’s, for all the dissolution reactions are approximately temperature independent.$`^{14}`$

3) There also appears to have a convergence temperature for enthalpy changes in dissolution of organic compounds from gas $``$ water (see Appendix for more discussion). The temperature is not known, but the corresponding $`\mathrm{\Delta }H^{}`$ is about a few KJ/mol.

There have been many different ways of presenting experimental data on the thermodynamics of dissolution. Here we have summarized several key results from the literature which are useful in organizing and relating various data representations. It is important to note that many presentations are in fact concerning same experimental measurements. They are related by the basic thermodynamic formulae, as shown in the following theorems.

Baldwin’s Theorem

Based on the observations of Sturtevant15 and Privalov14, R.L. Baldwin gave the following theorem in 1986.$`^{16}`$ For a series of homologous reactions, designated by $`i`$, with temperature-independent $`\mathrm{\Delta }C_{p_i}`$, the following two statements are equivalent: 1) At one temperature $`T^{}`$, there are unique constants $`a`$ and $`b`$ such that:

$$\mathrm{\Delta }S_i=a+b\mathrm{\Delta }C_{p_i}$$
(1)
2) there is a unique temperature $`T_s^{}`$ at which $`\mathrm{\Delta }S(T_s^{})`$ is unique. That is, $`T_s^{}`$ and $`\mathrm{\Delta }S_i(T_s^{})`$ are independent of $`i`$. Hence:

$$\mathrm{\Delta }S_i(T)=\mathrm{\Delta }S(T_s^{})+\mathrm{\Delta }C_{p_i}ln(T/T_s^{})$$
(2)
Compare Eq. (1) with (2), we have $`a`$ = $`\mathrm{\Delta }S(T_s^{})`$ and $`b`$ = $`ln(T^{}/T_s^{})`$. From now on, we will follow Lee$`^{13}`$ and refer to the $`\mathrm{\Delta }S_i`$ versus $`\mathrm{\Delta }C_{p_i}`$ relation in Eq. (1) as a SMPG plot named after Sturtevant, Murphy, Privalov, and Gill.$`^{10,15}`$ Baldwin’s theorem states that a linear SMPG plot is equivalent to having a convergence temperature for $`\mathrm{\Delta }S_i`$, and the intersection $`a`$ in the SMPG plot is the unique $`\mathrm{\Delta }S^{}`$. Conversely, the existence of a convergence temperature indicates a linear SMPG plot at every temperature.

Lee’s Theorem

Lee’s theorem$`^{13}`$ is an application of the mathematical property of a bilinear function: if a function $`f(x,y)`$ is linear as a function of either $`x`$ or $`y`$, as well as a function of their product, $`xy`$, then there exists a $`y^{}`$ at which $`f(x,y^{})`$ is independent of $`x`$. To express this by equation:

$$f(x,y)=a+bx+cy+dxy=d(x+\frac{c}{d})(y+\frac{b}{d})+\frac{adbc}{d}$$

so when $`y=y^{}=b/d`$, $`f(x,y^{})`$ is independent of $`x`$.

The bilinear function has to be a linear function of either variables when the other one is at a fixed value. To apply this result to the thermodynamics of dissolution, we identify the logarithmic temperature and a certain molecular parameter as the variables $`x`$ and $`y`$, and the entropy (or enthalpy or Gibbs free energy) of dissolution as the function $`f(x,y)`$. Hence, at a given temperature, the entropy is a linear function of the molecular parameter $`X`$. Lee proposed the parameter to be a size measure of the solute molecule. In fact, Lee’s proposal of

$$\mathrm{\Delta }S_i(T)=a_s+b_sX_i$$
(3)
is a structural interpretation of the linear SMPG plot given in (1). The subscript “s” here stands for entropy. (There is a similar relation for enthahlpy, with its respective $`a_h`$ and $`b_h`$.$`^{13}`$) It is well known that $`\mathrm{\Delta }C_p`$’s are proportional to molecular size of a hydrophobic solute. It is also worth noting that a linear relationship due to substituents of a series of functional groups is widely observed in organic reactions.$`^{17}`$

Therefore, the bilinear argument immediately leads to the existence of a convergence temperature from the linear entropy relation in Eq. (3). Combining Lee’s and Baldwin’s theorems, it can be shown that the linear entropy relation is sufficient but not necessary for generating a convergence temperature.

The BMDW Theorem

This theorem, which was given by Baldwin and Muller,$`^{18}`$ and also independently by Doig and Williams,$`^{19}`$ establishes an intrinsic relationship between the three convergence temperatures for entropy, enthalpy, and Gibbs free energy. Consider a series of homologous reactions designated by $`i`$ and each has a temperature-independent $`\mathrm{\Delta }C_{p_i}`$ of its own. If two out of the three thermodynamics quantities (entropy, enthalpy, and Gibbs free energy) have convergence temperatures, then there is a convergence temperature for the third quantity. If we denote by $`T_s^{}`$, $`T_h^{}`$, and $`T_g^{}`$ these convergence temperatures, we have:

$$T_g^{}T_h^{}T_g^{}ln(T_g^{}/T_s^{})=0$$
(4)
and also a relation for the corresponding thermodynamic quantities:

$$\mathrm{\Delta }G(T_g^{})=\mathrm{\Delta }H(T_h^{})T_g^{}\mathrm{\Delta }S(T_s^{})$$
(5)
This theorem was used by Baldwin and Muller$`^{18}`$ to explain the intriguing fact of approximately equal convergence temperatures $`T_s^{}`$ and $`T_h^{}`$ for a set of protein folding reactions. It is discovered that this set of proteins have approximately equal melting temperatures, i.e., $`T_g^{}`$ = 331 $`\pm `$ 9 K and $`\mathrm{\Delta }G(T_g^{})`$ $``$ 0. Using the same theorem, Doig and Williams$`^{19}`$ have reached a similar conclusion. They pointed out that most proteins of the same set have approximately equal $`\mathrm{\Delta }G`$ per residue at room temperature. These two arguments are consistent with each other if we note that all the proteins in the set have approximately same $`\mathrm{\Delta }C_p`$ per residue.

A Physical Model Based on SPT

The basic results of SPT are summarized here. It is not necessary for the readers to know the technical details of SPT;$`^5`$ rather, our model starts with Eqs. (6) and (7) below. SPT is a rigorous statistical mechanical theory for liquids made of hard spheres. The theory also provides an approximated formula for calculating the free energy of dissolving a hard spherical solute from gas phase into a hard spherical solvent, which is equivalent to introducing a spherical cavity in the hard sphere liquid. The free energy, which is a function of the radius of solute as well as the density and the radius of pure solvent, has three dominant terms. The first term is independent of solute size; it is associated with the reduction of conformational space of the liquid upon introducing a solute molecule of zero physical size into the solvent. With the physical point being present in the middle of the solution, no solvent molecule can occupy the same point in space. Hence such conformations are no longer accessible. This term is a function of the density and the radius of pure solvent only. The second term is proportional to the square of the radius of the solute, i.e., its surface area. And the third term is proportional to the cubic power of the radius, i.e., its volume, and external pressure. It has been repeatedly demonstrated, both theoretically and empirically, that the volume term is negligible for any molecular size cavity, and the dominant effect is from the surface term which defines a proportional coefficient called surface tension. For hard spherical solvent and solute, SPT gives an expression for the surface tension, and also a minor correction term on the free energy due to different curvature of the solute molecule.

The free energy is given as:$`^5`$

$$\frac{\mathrm{\Delta }G}{k_BT}=ln(1\xi )+\frac{r^2}{a^2}\left[\frac{3\xi }{1\xi }+\frac{9}{2}\left(\frac{\xi }{1\xi }\right)^2\right]$$
(6)
where $`k_B`$ is the Boltzmann constant, $`T`$ is temperature in Kelvin, $`a`$ is the radius of solvent molecule, $`\xi `$ is the packing density for pure solvent (the volume fraction occupied by the hard spheres), and $`r`$ is the radius of solute molecule.

Our model starts with a very simple premise:

$$\mathrm{\Delta }G(T)=\delta (T,\xi )+\sigma (T,\xi )r^2$$
(7)
where $`\delta `$ and $`\sigma `$ are solvation energies for a point and surface tension, respectively. Thermodynamics for solvation will depend on the functional form of $`\delta `$, $`\sigma `$, and implicit temperature dependence of $`\xi `$. Eq. (7) offers the possibility for connecting thermodynamics of solvation with the thermal expansion coefficient of pure solvent (see below).

From Eq. (7), it is straightforward to obtain:

$$\mathrm{\Delta }S(T)=\delta _T\delta _\xi \xi _T(\sigma _T+\sigma _\xi \xi _T)r^2$$
(8)
$$\mathrm{\Delta }H(T)=\delta T\delta _TT\delta _\xi \xi _T+(\sigma T\sigma _TT\sigma _\xi \xi _T)r^2$$
(9)
where subscripts “$`T`$” and “$`\xi `$” represent partial derivatives with respect to these variables.

Let’s now estimate the magnitudes of the various terms using values for water$`^3`$ (SPT parameters for water at room temperature are $`\xi `$ = 0.363 and $`a`$ = $`r_w`$ = 1.38 $`\AA `$.) $`\xi _T`$ is proportional to the thermal expansion coefficient of the pure solvent: $`\xi _T=\alpha \xi `$, where $`\alpha `$ for water is very small. At 1 atm, the values of $`\alpha `$ for water range from $`0.064\times 10^3`$ at $`0^{}`$C to $`0.7\times 10^3`$ at $`100^{}`$C, and $`\alpha `$ = $`0.257\times 10^3`$ at room temperature ($`25^{}`$C). The thermal expansion coefficients for non-associative liquids are around $`1\times 10^3`$ at 1 atm and between $`0^{}`$C to $`100^{}`$C, so these values are not very different from that of water.

Quantitatively, compare Eqs. (6) and (7), we have:

$$|\delta _\xi \xi _T|=\left|\frac{k_BT\alpha \xi }{1\xi }\right|k_B\left(\frac{300\times 0.257\times 10^3\times 0.363}{10.363}\right)=0.04k_B$$

while

$$\delta _T=k_Bln(1\xi )=0.5k_B>>|\delta _\xi \xi _T|$$

where we have used $`T`$ = $`300K`$. Similarly by simple differentiation,

$$\sigma _\xi =k_BT\frac{3(1+2\xi )}{r_w^2(1\xi )^3}=20k_BT/r_w^2.$$

Hence $`|\sigma _\xi \xi _T|`$ $``$ $`0.56k_B/r_w^2`$, while $`\sigma _T`$ = $`3.17k_B/r_w^2`$ $`>>`$ $`|\sigma _\xi \xi _T|`$.

The Heat Capacity $`\mathrm{\Delta }C_p`$

The heat capacity can be obtained from either Eq. (8) or Eq. (9):

$$\mathrm{\Delta }C_p=T(\delta _{TT}+2\delta _{T\xi }\xi _T+\delta _{\xi \xi }\xi _T^2+\delta _\xi \xi _{TT})T(\sigma _{TT}+2\sigma _{T\xi }\xi _T+\sigma _{\xi \xi }\xi _T^2+\sigma _\xi \xi _{TT})r^2$$
(10)
While the values for $`\alpha `$ are not very different between water and non-associative liquids, there is a dramatic difference between $`\alpha `$’s dependence on temperature. This difference contributes a large term to the $`\mathrm{\Delta }C_p`$ of solvation in water, in contrast to non-associative solvent. In other words, even though $`\alpha `$ for water is very small, its temperature dependence is quite large, in contrary to most organic solvents. Hence we will neglect contribution from $`\xi _T`$ but shall keep the terms with $`\xi _{TT}`$. Thus we have:

$$\mathrm{\Delta }C_p=T\xi _{TT}(\delta _\xi +\sigma _\xi r^2)$$
(11)
where, according to SPT (i.e., Eqs (6) and (7)), both $`\delta _\xi `$ and $`\sigma _\xi `$ are explicitly proportional to T. Therefore if $`\xi _{TT}`$ is proportional to $`T^2`$, then $`\mathrm{\Delta }C_p`$ will be approximately temperature independent.

It is clear from the above argument that, for a non-associative solvent, the solvation of a inert solute should have very small $`\mathrm{\Delta }C_p`$ since non-associative liquid has almost zero $`\xi _{TT}`$. For water, $`\delta _\xi `$ = 1.57$`k_BT`$ and $`\sigma _\xi r^2`$ = 20.0$`k_BT(r/r_w)^2`$. Therefore, Eq. (11) gives a linear relationship between $`\mathrm{\Delta }C_p`$ and molecular surface area $`r^2`$ with almost zero intersection when $`r>r_w`$. Note, however, that for hydrocarbon solutes like propane and isobutane, there are significant solute-solvent interactions which contribute to the overall $`\mathrm{\Delta }C_p`$.$`^6`$ Therefore, our present result is only semi-quantitative and has to be augmented with such interactions when applied to real experimental data.

$`\mathrm{\Delta }S^{}`$ and $`\mathrm{\Delta }H^{}`$

We now return to Eq. (8). According to Lee’s theorem:$`^{13}`$

$$\mathrm{\Delta }S^{}=\delta _T\delta _\xi \xi _T$$
(12)
Note that in order to compare calculation from SPT with experimental measurements, we have to obtain the entropy from experimental measurements according to molarity concentration scale.$`^{5,11}`$ When this was done,$`^{13}`$ Lee found that the calculation given by SPT compares favorably with experimental results.

Similarly, we have:

$$\mathrm{\Delta }H^{}=\delta T\delta _TT\delta _\xi \xi _T=T\delta _\xi \xi _T$$
(13)
the second equality is because $`\delta `$ is a linear function of $`T`$ (see Eqs. (6) and (7)). Numerically, $`T\delta _\xi \xi _T`$ is about 0.034 $`k_BT`$, that is, 2.5 kJ/mol. It should be noted that Eq. (13) neglects contribution from solute-solvent interaction. It is known that different hydrocarbons, for example aromatics and aliphatics, have different $`\mathrm{\Delta }H^{}`$. On the other hand, inert gases could be used as a test for the present model.

Thus, according to our analysis, the converging values for entropy and enthalpy are the consequence of point solvation energy. This is an interesting conjecture. A rigorous statistical mechanical treatment of this problem seems possible, but has never been developed. In general, point solvation energy is dependent upon whether it is solid, or liquid, or gas, from which the solutes are transferred. A quantitative theory might be able to explain the small differences between the three SMPG plots.$`^{10}`$ (Note: after correction according to molarity scale, the intersection, $`\mathrm{\Delta }S^{}`$, for gas $``$ water dissolution is increased by about 60 J per degree per mole, see ref. 10).

$`T_s^{}`$ and $`T_h^{}`$

It is natural to suspect that the unique convergence temperature, $`T_s^{}`$, for entropy of all three different groups of dissolution transferred from either solid, liquid, or gas, is due to some intrinsic properties of water. $`T_s^{}`$ is the temperature at which $`\mathrm{\Delta }S(T)`$ in Eq. (8) equals $`\mathrm{\Delta }S^{}`$ given in Eq. (12). That is:

$$\sigma _T+\sigma _\xi \xi _T=0$$
(14)
i.e.,

$$\xi _T(T_s^{})=\alpha (T_s^{})\xi (T_s^{})=\sigma _T/\sigma _\xi .$$

By an approximated calculation, $`\sigma _T/\sigma _\xi `$ = 0.158/T. Hence:

$$T_s^{}\alpha (T_s^{})\xi (T_s^{})=0.158.$$

This is consistent with the laboratory measurements $`T_s^{}`$ = 383 K, $`\xi `$ = 0.363 and $`\alpha `$ = $`1.1\times 10^3`$.

Similarly for $`T_h^{}`$, from Eqs. (9) and (13) we have:

$$\sigma T\sigma _TT\sigma _\xi \xi _T=0.$$
(15)
According to SPT, the first two terms cancel each other. The third term, as we have indicated, is indeed quite small. Thus unfortunately, Eq. (15) is buried in our various approximations and fails to provide an estimation for $`T_h^{}`$. A more accurate estimation is required to obtain the convergence temperature for $`\mathrm{\Delta }H`$.

Discussion

With the physical insight provided by SPT, we now attempt to answer some key questions concerning the thermodynamics of hydrophobic solvation and the hydrophobic effect.

1) What is the Hydrophobic Effect?

Ever since Kauzmann’s seminal paper on hydrophobic effect,$`^{20}`$ people have believed that hydrophobic effect is mostly due to the reorganization of hydrogen bonds among the solvent molecules around the solute, and the contribution of direct interaction between solvent and solute is rather minimal. So what is the relation between the hydrophobic effect and the solvation thermodynamics for hydrophobic solute in water? To address the question, one has to be precise about the meaning of “hydrophobic effect”. There is an experimental (thermodynamic) side and there is a structural (theoretical) side of conventional wisdom on “hydrophobicity”. The experimental side is that non-polar solutes in water have very low solubility and the dissolution has large heat-capacity changes, in contrast to the dissolution in organic solvent. The structural side is that hydrogen bonding arrangement has been altered when an non-polar solute is dissolved in water.

From our analysis, it seems that low solubility and large $`\mathrm{\Delta }C_p`$ in fact stem from two distinct sources.$`^{21}`$ While the large $`\mathrm{\Delta }C_p`$ is associated with rearrangement of water molecules, the low solubility is primarily due to the geometric properties of water molecules. This suggestion is consistent with our understanding of entropy-enthalpy compensation, which says that the ability of rearrangement of solvent should only have minor effect on solubility.$`^7`$ The more dynamic aspect of water will be reflected only through quantities like entropy, enthalpy, and heat capacity. The reorganization process within solvent contributes to entropy change through heat capacity $`\mathrm{\Delta }C_p`$, which in turn is related to the peculiar large temperature dependence of $`\alpha `$, the thermal expansivity, of water. Hence it seems legitimate to identify $`\mathrm{\Delta }C_p`$ with the hydrophobic effect.

However, could one simply identify the $`\mathrm{\Delta }C_p`$ term as hydrophobic free energy? This indeed is the central issue behind the work by Murphy et al.$`^{10}`$ It was, of course, well recognized that the $`\mathrm{\Delta }C_p`$ term is not uniquely determined until an appropriate reference temperature(s) is chosen.$`^{22}`$ Murphy et al. proposed the using of $`T_s^{}`$ and $`T_h^{}`$ as reference temperatures,$`^{10}`$ and thereafter a $`\mathrm{\Delta }C_p`$ term was uniquely defined. However, as we have seen, the existence of convergence temperatures and their values are not the hallmark of the reorganization of associative solvent. We suspect that many other solvation processes might also have such convergent properties. The basis for the existence of convergence temperatures is Eq. (7), and it is clear that this equation is not unique for hydrophobic solvation (more discussion later).

2) What is the Role of Hydrogen Bond?

An inevitable objection to our approach from many readers will be the complete neglect of hydrogen bond which has central importance in Kauzmann’s structural model for hydrophobicity.$`^{20}`$ We would like to emphasize that we do accept the hydrogen bond in water as the structural base for hydrophobicity, but we want to seek the specific thermodynamic aspect or aspects of the hydrogen-bond structure which are responsible for the thermodynamics of hydrophobic solvation. The large temperature dependence of $`\alpha `$ no doubt is a manifestation of hydrogen bonding reorganization in water, the ultimate source of hydrophobicity. However, other aspects of the water molecules might also be relevant or even crucial; for example, the tetrahedral chemistry of hydrogen bonding,$`^{20}`$ or more generally the non-isotropic pair-wise interaction between two water molecules,$`^{23}`$ and the high ratio between physical volume and thermodynamic volume of water.$`^3`$ In our model, SPT indeed uses all these properties of water, though not explicitly. The fact is that water molecule in the SPT model has large thermodynamic volume, i.e., low packing density, but at the same time a small radius leads to high solvation number around a cavity.$`^3`$ This indicates that the solvent molecule is not isotropic, and there are preferences for these molecules to surround a cavity, which is exactly the Kauzamnn’s argument! A crude analogy will be a wedge-shaped molecule, and that is quite consistent with water molecules.

3) What is the basis of convergence temperature? As pointed out by Lee,$`^{13}`$ the presence of convergence temperature is due to some kind of linear free energy dependence on molecular substituents, i.e., Eq. (7).

4) What determine the magnitudes of $`\mathrm{\Delta }S^{}`$ and $`\mathrm{\Delta }H^{}`$? They are determined primarily by the thermodynamics of solvation of zero size point solutes.

5) Why is there a large $`\mathrm{\Delta }C_p`$? The $`\mathrm{\Delta }C_p`$ of hydrophobic solvation stems from the peculiar large temperature dependence of thermal expansion coefficient, $`\alpha /T`$.

6) How to Obtain Molecular Interaction Energy From the Thermodynamic Data?

This is an age-old question. Twenty-five years ago, T.H. Benzinger proposed a new definition for enthalpy of chemical reaction.$`^{24}`$ His argument was that for a chemical reaction with non-zero $`\mathrm{\Delta }C_p`$, there would be no unique heat of formation for the reaction. Since the true heat of formation is mechanical (athermal), Benzinger suggested to use $`\mathrm{\Delta }H(0)`$ at zero Kelvin as the “true” heat of formation, and argued that the remains of free energy should be lumped into one term:

$$\mathrm{\Delta }W(T)=\mathrm{\Delta }H(0)\mathrm{\Delta }G(T).$$

When taking into account the fact $`\mathrm{\Delta }S(0)`$ = 0, he obtained:

$$\mathrm{\Delta }W(T)=_0^T\mathrm{\Delta }C_p(X)\left(\frac{T}{X}1\right)𝑑X$$

In some sense, what Benzinger did was similar to what Murphy et al. did. They were both trying to divide the total free energy into a part with direct (mechanical) interaction and the rest part with surrounding effect due to thermodynamics. They both realized that $`\mathrm{\Delta }C_p`$ term was related to the latter since it characterizes the fluctuation in enthalpy due to thermal agitation. The crucial question is of course whether it is possible to find an appropriate reference temperature(s) based on purely thermodynamic analysis without any molecular model. Chan and Dill recently have extensively discussed this issue, and they concluded that without a molecular model, purely thermodynamic analysis would not provide much meaningful result.$`^{12}`$

7. The Validity of Using SPT to Model Hydrophobic Solvation.

Let’s now reiterate the rationales for using SPT to model solvation in associative solvents. In SPT, the reorganization of solvent is considered implicitly through experimental data on $`\xi `$ as function of temperature. This approach is consistent with the assertion that the ultimate reason for reorganization in a solvent is its temperature dependence as a pure liquid.$`^7`$ The presence of hydrogen bonds between solvent molecules is manifested in the experimental data on water density and its temperature dependence. SPT, of course, completely neglects the soft interaction between solute and solvent.

Our second defense is based on a recent analysis of how thermodynamic systems respond to small perturbations.$`^7`$ It has been shown that if we classify thermodynamic quantities by the orders of derivative of free energy, there is a relationship between thermodynamic of perturbation and the thermodynamics of unperturbed system. Since one can treat solvation as a perturbation, one only needs the thermodynamics of one order higher for pure solvent in order to calculate the thermodynamics of solvation. For example, $`\alpha _T`$ of solvent gives $`\mathrm{\Delta }C_p`$ of solvation. All the structural changes in hydrogen bonding will be captured in these thermodynamic quantity of pure water, and our analysis made use of them.

Finally, we would like to emphasize that we are not attempting to use SPT to model the properties of water; rather, we are merely using SPT to relate the thermodynamics of dissolution of an inert solute in water to that of pure water. In the past, many models for water which are based on multi-state of water conformation have been successful in providing calculations for thermal expansion coefficient, $`\alpha `$, but have failed to deal with its temperature dependence.$`^{25}`$ This situation is completely in accord with our analysis.

Acknowledgments

I thank Buzz Baldwin, B.K. Lee, and John Schellman for many detailed and helpful discussion in the past several years. I wish to specially dedicate this work to Professor Baldwin on the occasion of his 70th birthday. His pioneer work on helical peptide, protein hydrogen exchange, DNA flexibility, and hydrophobic effect have been guiding lights for my research.

References

1. Stillinger, F. H. (1973) J. Solution Chem., 2, 141-158. 2. Postma, J.P.M., Berendsen, H.J.C., & Haak, J.R. (1982) Faraday Symp. Chem. Soc., 17, 55-67. 3. Lee, B. (1985) Biopolymers, 24, 813-825. 4. Morton, A., Baase, W.A., & Matthews, B.W. (1995) Biochemistry, 34, 8564-8575. 5. Reiss, H. (1965) Adv. Chem. Phys. 9, 1-84. 6. Lee, B. (1991) Biopolymers, 31, 993-1008. 7. Qian, H. & Hopfield, J.J. (1996) J. Chem. Phys., 105, 9292-9298. 8. Grunwald, E. (1984) J. Am. Chem. Soc., 196, 5414-5420. 9. Ben-Naim, A. (1975) Biopolymers, 14, 1337-1355. 10. Murphy, K., Privalov, P.L., & Gill, S.L. (1990) Science, 247, 559-561. 11. Ben-Naim, A. (1978) J. Phys. Chem., 82, 792-803. 12. Chan, H.S. & Dill, K.A. (1997) Ann. Rev. Biophys. Biomol. Struct., 26, 423-457. 13. Lee, B. (1991) Proc. Natl. Acad. Sci. USA, 88, 5154-5158. 14. Privalov, P. (1979) Adv. Protein Chem. 33, 167-241. 15. Sturtevant, J.M. (1977) Proc. Natl. Acad. Sci. USA, 74, 2236-2240. 16. Baldwin, R.L. (1986) Proc. Natl. Acad. Sci. USA, 83, 8069-8072. 17. Leffler, L. & Grunwald, E. (1963) Rates and Equilibria of Organic Reactions, John-Wiley & Sons, New York. 18. Baldwin, R.L. & Muller, N. (1992) Proc. Natl. Acad. Sci. USA, 89, 7110-7113. 19. Doig, A.J. & Williams, D.H. (1992) Biochemistry, 31, 9371-9375. 20. Kauzmann, W. (1959) Adv. Protein Chem., 14, 1-63. 21. Lee, B. (1994) Biophys. Chem., 51, 271-278. 22. Becktel, W.J. & Schellman, J.A. (1987) Biopolymers, 26, 1859-1877. 23. Ben-Naim, A. (1992) Statistical Thermodynamics for Chemists and Biochemists, Plenum Press, New York. 24. Benzinger, T.H. (1971) Nature, 229, 100-102. 25. Kauzmann, W. (1975) Colloq. Inter. C.N.R.S., 246, 63-71.

Appendix: A possible $`\mathrm{\Delta }H^{}`$ in gas to water dissolution?

Let’s proceed with the dissolution of a gas into water (g $``$ w) by first liquefying the gas (g $``$ l) and then transferring the liquid into water (l $``$ w):

$$\mathrm{\Delta }H_{gw}=\mathrm{\Delta }H_{gl}+\mathrm{\Delta }H_{lw}.$$

There seems to be a convergence temperature for dissolution enthalpy change from liquid $``$ water at about 2$`0^{}`$C, and corresponding $`\mathrm{\Delta }H^{}`$ is about few KJ/mol.$`^{16}`$ According to Trouton’s rule, we can write

$$\mathrm{\Delta }H_{gl}=88T_bJ/mol$$

It is well known that the boiling temperature $`T_b`$ is scaled with molecular size, hence $`\mathrm{\Delta }H_{gw}`$ and $`\mathrm{\Delta }H_{lg}`$ should have similar converging $`\mathrm{\Delta }H^{}`$ but at different temperatures. Mathematically, if we have:

$$\mathrm{\Delta }H_{lw}=\mathrm{\Delta }H^{}+b_1(T)r^2$$

where $`b_1`$ = 0 when T = 2$`0^{}`$C, and

$$\mathrm{\Delta }H_{gl}=b_2(T)r^2$$

then:

$$\mathrm{\Delta }H_{gw}=\mathrm{\Delta }H^{}+[b_1(T)+b_2(T)]r^2$$

the convergence temperature is the value of T at which $`b_1+b_2`$ = 0.


format: text

text:
# REFERENCES

Nonlinear Thermomagnetic Waves in Type II Superconductors N.A. Taylanov
Theoretical Physics Department and Institute of Applied Physics,
National University of Uzbekistan,
E-mail: taylanov@iaph.silk.org

bstract

This work is devoted to analysis of the nonlinear dynamics of interrelated thermal and electromagnetic perturbations in a superconductor caused by dissipative effects involved in the viscous motion of the magnetic flux. The structure and evolution of the stationary thermomagnetic waves propagating in the superconductor is studied based on a system of nonlinear equations describing dynamics of the instability development in the system studied. The apperance of these structures is qualititively described and the wave propagation velocity is estimated. The affect a thermoactivated flux creep to the structures of the nonlinear thermomagnetic wave propagation is studied. It is shown that taking into account a nonlinear dependency of the current density $`j`$ on the electric field strength $`E`$ would not qualitatively change the main results, since the character of the equilibrium state on the phase plane remains esentially the same.

At present superconducting systems with high critical field strengs and current densites are widely implemented in various advenced technologies. Hovever, successful operation of the superconducting materials is only possible provided that special measures are taken to prevent a system from the thermal or magnetic breakage of superconductivity and transition to a normal (resistive)state . For this reason, one of the main problems in the investigation of properties of superconductors is that of predicting the superconducting state breacage caused by dissipative and nonlinear effects related to viscous motions of the magnetic flux. This explains a considerable interest in the study of dissipative and nonlinear effects in superconductors that has apose in recent years.

This work is devoted to analysis of the nonlinear dynamics of interrelated thermal and electromagnetic perturbations in a superconductor caused by dissipative effects involved in the viscous motion of the magnetic flux. The structure and evolution of the stationary thermomagnetic waves propagating in the superconductor is studied based on a system of nonlinear equations describing dynamics of the instability development in the system is studied.

The evolution of thermal and electromagnetic perturbations in a superconductor is described by a set involving the thermal conductivity equation \[2-5\]

$$\nu \frac{dT}{dt}=(\kappa T)+\stackrel{}{j}\stackrel{}{E},$$
(1)
the Maxwell equation

$$\frac{4\pi }{c^2}\frac{d\stackrel{}{j}}{dt}=\mathrm{\Delta }\stackrel{}{E},$$
(2)
and the related equation of the critical state

$$\stackrel{}{j}=\stackrel{}{j}(T,\stackrel{}{H})+\stackrel{}{j}_r(\stackrel{}{E}).$$
(3)
Here $`\nu =\nu (T)`$ is the heat capacity, $`\kappa =\kappa (T)`$ is the thermal conductivity coefficient, $`\stackrel{}{j}_c`$ is the critical current density, and $`\stackrel{}{j}_r`$ is the current density in the resistive state.

The above system is essentially nonlinear because the right-hand part of Eq.(1) contains a term describing the Joule heat evolution in the region of a resistive phase. Such a set (1)-(3) of nonlinear parabolic differential equations in partial derivatives has no exact analytical solution.

Let us consider a planar semi-infinite sample $`(x>0)`$ placed in an external magnetic field $`\stackrel{}{H}=(0,0,H_e)`$ growing at a constant rate $`{\displaystyle \frac{d\stackrel{}{H}}{dt}}=const`$. According to Maxwel’s equation (2), there is a vortex electric field $`\stackrel{}{E}=(0,E_e,0)`$ in the sample, directed parallel to the current density $`\stackrel{}{j}`$: $`\stackrel{}{E}\stackrel{}{j}`$; where $`H_e`$ is the amplitude of the external magnetic field and $`E_e`$ is the amplitude of the external electric field.

For the automodel propagating waves of the type $`\xi (x,t)=xvt`$ the initial differential Eqs. (1)-(3) acquire the following form

$$\nu v\frac{dT}{d\xi }=\frac{d}{d\xi }\left[\kappa \frac{dT}{d\xi }\right]+jE,$$
(4)
$$\frac{dE}{d\xi }=\frac{4\pi v}{c^2}j,$$
(5)
$$E=\frac{v}{c}H.$$
(6)
A set of boundary conditions corresponding to Eqs. (4)-(6) is as follows:

$$\begin{array}{c}T(\xi +\mathrm{})=T_0,\frac{dT}{d\xi }(\xi \mathrm{})=0,\hfill \\ \\ E(\xi +\mathrm{})=0,E(\xi \mathrm{})=E_e,\hfill \end{array}$$
(7)
Let us use the Bean-London model of the critical state for the dependence $`j_c(T,H)`$; $`{\displaystyle \frac{dj_c}{dH}}=0`$ . According to this model the dependence $`j_c(T)`$ can be represented as $`j_c(T)=j_0[1a(TT_0)]`$: where $`j_0`$ is the equilibtium current density, $`a`$ is the constant parameter and $`T_0`$ is the temperature of the cooling medium. A characteristic field dependence of $`j_r(E)`$ in the region of sufficiently strong electric fields $`(E>E_f)`$ can be aproximated by a piecewise linear function $`j_r\sigma _fE`$, where $`\sigma _f`$ is the effictive conductivity. In the field of small field strengths $`(E<E_f)`$ we assume a relationship $`j_r(E)j_1\mathrm{ln}{\displaystyle \frac{E}{E_0}}`$,to be valid with $`j_1`$ being a characteristic local current density scatter (related to a pinning force inhomogeneity) on the order of $`j_10,01j_c`$ and $`E_0`$ being a constant. This relationship between $`j_r`$ and $`E`$ is due to a thermoactivated early experiments (see,e.g.). The thermoactivated flux creep is principal mechanism of the energy dissipation during the magnetic flux penetration deep into the sample. This factor must affect the character of the nonlinear thermomagnetic wave propagation in the system studied.

Excluding the variables $`T(\xi )`$ and $`H(\xi )`$, using Eqs.(4) and (6) , and taking into account the boundary conditions (7), we obtain an equation describing the electric field $`E(\xi )`$ distribution:

$$\frac{d^2E}{dz^2}+\beta \left[1+\frac{j_1}{\sigma _dE}\tau \right]\frac{dE}{dz}+\beta ^2\tau \left[\frac{j_0}{\sigma _d}+\frac{j_1}{\sigma _d}\mathrm{ln}\frac{E}{E_0}\right]=\frac{E^2}{2E_\kappa }.$$
(8)
The dimensionless variables used in Eq. (8) are defined as follows:

$$z=\frac{\xi }{L},\beta =\frac{vt_\kappa }{L},t_\kappa =\frac{\nu L^2}{\kappa },\tau =\frac{4\pi \sigma _d\kappa }{c^2\nu },E_\kappa =\frac{\kappa }{aL^2},L=\frac{cH_e}{4\pi j_0},$$

where $`L`$ is the depth of the magnetic flux penetration into the sample, $`\sigma _d`$ is the differential conductivity. The corresponding equation of state is obtained using a relationship

$$\mathrm{\Omega }E^2=X(E)=1+\frac{j_1}{j_0}\mathrm{ln}\frac{E}{E_0},$$
(9)
where $`\mathrm{\Omega }={\displaystyle \frac{\sigma _d}{2\beta ^2\tau j_0E_\kappa }}`$.

As seen from the plots of $`X(E)`$ versus $`E`$ (see figure), there exists a single point of intersection of the curves $`y=\mathrm{\Omega }E^2`$ and $`y=1+(j_1/j_0)\mathrm{ln}(E/E_0)`$, which corresponds to a single stable equilibrium state. The stability of this state is determined by the sign of derivative $`d^2E/dz^2`$ in the vicinity of the equilibrium point. The wave velocity $`v_E`$ can be determined from stationary $`({\displaystyle \frac{dE}{dz}}=0)`$ Eq.(8) with an allowance for the boundary conditions :

$$v_E=\frac{L}{t_\kappa }E_e\left[2\tau \frac{j_0E_\kappa }{\sigma _d}\left(1+\frac{j_1}{j_0}\mathrm{ln}\frac{E_e}{E_0}\right)\right]^{1/2}.$$
(10)
Now we can use Eq. (6) and readily derive an equation describing the field distribution for a nonlinear $`H`$ \- wave. The wave velocity $`v_H`$ of this $`H`$\- wave is given by the formula

$$v_H=\frac{cE}{H_e}\mathrm{exp}\left[\frac{\sigma _d}{2\tau j_1}E_\kappa \left(\frac{LH_e}{ct_\kappa }\right)^2\frac{j_0}{j_1}\right].$$
(11)
As seen, the wave velocity $`v_H`$ exponentially increases with the field amplitude $`H_e`$. For a suffuciently small amplitude $`H_e<H_a`$, where

$$H_a=\frac{ct_\kappa }{L}\left[\frac{2\tau j_0E_\kappa }{\sigma _d}\right]^{1/2}$$
(12)
the wave velocity is negigible small, which corresponds to the case of the thermoactivated flux creep. For $`H_e=H_a`$, the wave propagates at a finite constant velocity. In this case, the maximum heating of a superconductor in the region immiadiately in front of the wave is described by the relationship

$$\frac{TT_0}{T_0}=\frac{H_e^2}{8\pi \nu T_0}.$$
(13)
In conluding, we must note that taking into account a nonlinear dependency of the current density $`j`$ on the electric field strength $`E`$ would not qualitatively change the main results, since the character of the equilibrium state on the phase plane remains esentially the same (see ).


text:
# How Many Genes Are Needed for a Discriminant Microarray Data Analysis ?

## Introduction

There are two types of microarray experiments. In the first type, samples are not labeled, but gene expression levels are followed with time. The goal is to find genes whose expression levels move together. In the second type, samples are labeled as either normal or diseased tissues. The goal is to find genes whose expression levels can distinguish different labels. Using terminology from machine learning, the data analysis of the first type of experiment is “unsupervised”, whereas that of the second type is “supervised” (see, e.g., \[Haghighi, Banerjee, Li, 1999\]). The prototype of the first type of analysis is cluster analysis, whereas that of the second type is discriminant analysis.

We focus on the data analysis of the second type, using the leukemia data from Whitehead/MIT’s group \[Golub et al, 1999\]. The question we are addressing is how many gene expressions are needed for a discriminant analysis. In this data set, 7129 gene expressions are measured on 38 sample points. In the original analysis \[Golub et al, 1999\], 50 out of more than 7000 gene expressions are used for discriminant analysis. We ask whether 50 genes are still too large a number, considering that the sample size is only 38. This seemingly simple question may receive two answers from two different perspectives: model selection (see, e.g., \[Burnham, Anderson, 1998\]) and model averaging (see, e.g., \[Geisser, 1974\]). We will address the two separately below.

## Model Selection

The goal of model selection is to pick the model among many possible models that achieves the best balance between data-fitting and model complexity. A perfect data-fitting performance by a complicated model with many parameters can well be an example of overfitting. On the other hand, a simplistic model with few parameters that fits the data poorly is an example of underfitting. There are two proposals for achieving the balance between data-fitting and model complexity: Akaike information criterion (AIC) \[Akaike, 1974; Parzen, Tanabe, Kitagawa, 1998; Burnham, Anderson, 1998\] and Bayesian information criterion (BIC) \[Schwarz, 1976; Raftery, 1995\]. These two quantities are defined as:

$$\text{AIC}=2\mathrm{log}(\widehat{L})+2K,\text{ and}\text{BIC}=2\mathrm{log}(\widehat{L})+\mathrm{log}(N)K,$$
(1)
where $`\widehat{L}`$ is the maximized likelihood, $`K`$ the number of free parameters in the model, and $`N`$ the sample size. High-order terms of $`O(1/N)`$ (for AIC) and $`O(1)+O(1/\sqrt{N})+O(1/N)`$ (for BIC) are ignored here for simplicity. The model with the lowest AIC or BIC is considered to be the preferred model.

Selecting genes relevant to a discriminant microarray data analysis may become an issue of model selection. But we have to be clear in what context this is the case. Model selection adjusts the number of parameters in the model, not the number of variables. Nevertheless, if we plan to combine variables (additively or any other functional form), each variable has one or perhaps more coefficient. Removing a variable removes the corresponding coefficient. In this context, variable selection is a special case of model selection. We illustrate this by the logistic regression/discrimination for our data set:

$$\text{Prob(AML)}=\frac{1}{1+e^{a_0_{jtopgenes}a_jx_j}}$$
(2)
where $`x_j`$ is the (log, normalized) gene expression level of gene $`j`$ among the top performing genes, and AML (acute myeloid leukemia) is one type of leukemia (ALL, the acute lymphoblastic leukemia, is another type). Using a linear combination of all 7129 genes in logistic discrimination requires 7130 parameters, whereas using one gene requires only 2 parameters. AIC/BIC defined in Eq.(1) will then compare multiple-gene models with single-gene models, and determine which scenario is better.

Our results are summarized in Table 1, where we list the type of model (which variables are additively combined in the logistic regression), number of parameters in the model, -2log of the maximum likelihood, AIC/BIC (absolute and relative), prediction rate in training set, and that in the testing set. The most striking result from this data set is that many genes are strong predictors for the leukemia class, consistent with the observation in \[Golub et al, 1999\]. This observation is confirmed by the following facts in Table 1: logistic regressions using the top 2, 5, 10, 22,and 37 genes all fit the data almost perfectly (as measured by the $`2\mathrm{log}(\widehat{L})`$ value) (note that the model is saturated when the number of variables used in the logistic regression is 37, since the number of parameters is then equal to the number of sample points); stepwise variable selection leads to only two genes (note that stepwise variable selection fails to find the best model, a single-gene model, because it is a local minimization/maximization procedure); the No.10 best-performing gene is only slightly worse than the No.1 performing gene: 35 vs. 36 correct predictions on the training set (though the No.100 and No.200 best-performing genes predict only 30 and 29 correctly on the training set); etc. This situation of strong prediction or easy classification is in contrast with the epidemiology and pedigree data used in human complex diseases, where strong predictors are rather rare \[Li, Sherriff, Liu, 2000; Li, Nyholt, 2001\].

Thest two rows in Table 1 are predictors/classifiers that do not use any gene expression levels. These are random guesses or null models. The first null model (#1) uses the proportion of AML in all samples as Prob(AML) (11/38 in our training set) as the guessing probability for AML. The second null model (#0) uses half-half probabilities. It is interesting to note that to beat both null, the number of genes in logistic regression can not exceed some upper limits: for null model 1, since we require AIC/BIC to be smaller (where $`p`$ is the number of genes used in logistic regression, and $`2\mathrm{log}(\widehat{L})`$ is assumed to be zero for the best-case scenario):

$`0+2(p+1)\text{AIC}`$ $`<`$ $`\text{AIC (fixed)}=47.728`$
$`0+3.637568(p+1)\text{BIC}`$ $`<`$ $`\text{BIC (fixed)}=49.365,`$ (3)
it sets $`p<22.8622`$, and $`p<12.5712`$, respectively. Similarly, to beat the random-guess model, we require:

$`0+2(p+1)\text{AIC}`$ $`<`$ $`\text{AIC (random)}=52.679`$
$`0+3.637568(p+1)\text{BIC}`$ $`<`$ $`\text{BIC (random)}=52.679`$ (4)
which are $`p<25.3425`$, and $`p<13.4813`$, respectively. These specific upper limits are directly related to the fact that the sample size is 38. They will move up and down with the sample size.

## Model Averaging

In the model selection framework, we can not use a logistic regression with too many genes because it may not improve the data-fitting performance enough to compensate for the increase of model complexity. This conclusion is correct when one model (e.g. a logistic regression that uses one gene) is compared to other alternative models (e.g. a logistic regression that uses, say, 10 genes). Nevertheless, it is possible to average/combine many different models each involving one gene. The restriction on model complexity during the model selection process does not apply to model averaging. Model averaging has also been discussed under names such as “committee machines”, “boosting weak predictors”, “mixture of experts”, “stacked generalization” (see, e.g., \[Ripley, 1996\]).

Without guidance from the model selection framework on the number of models (number of genes) to be included, we have tried several empirical approaches. We first examine whether there is a gap in data-fitting performance among top genes. Genes that do not fit the data should not be considered in model averaging. For this purpose, Fig.1 shows the $`2\mathrm{log}(\widehat{L})`$ for the top 1000 genes on the training set, as a function of the (log) rank. A linear fitting of $`2\mathrm{log}(\widehat{L})`$ on $`\mathrm{log}`$(rank) seems to fit the points well, and it is hard to see “better-than-average” genes except the first gene. This is not surprising since g4847 discriminates the 38 training sample points perfectly. Note that the linear trend in Fig.1 is equivalent to a power-law function in likelihood vs. rank plot ( $`\widehat{L}1/r^{2.56}`$). Such power-law is similar to the power-law rank-frequency plots observed in many social and natural data, also known as Zipf’s law \[Zipf, 1949; Li, 1997-2001\]. In any case, there is no discernible gap in Fig.1 that separates relevant and irrelevant genes.

We then check how the prediction rate on the training set correlates with that on the testing set. Fig.2 shows the error rates on both training (x-axis) and testing sets (y-axis) for the top 500 performing genes. The left plot in Fig.2 shows the mean square error, and the right plot shows the prediction errors. The left plot contains both information on the success rate of prediction and that of confidence of prediction, whereas the right plot contains only information on the success rate of prediction. Points along the diagonal line in Fig.2 exhibit similar error rates in the training and testing set, and thus are reasonable predictors. On the other hand, points well above the diagonal line indicate overtraining. The most reasonable predictor based on both the training and testing set is g1882 (on the other hand, the best predictor based on training set is g4847).

Finally we examine the model averaging performance with various numbers of models included, each being a single-gene logistic regression:

$$\text{Prob(AML)}=\underset{j\text{top genes}}{}w_j\left(\frac{1}{1+e^{a_jb_jx_j}}\right).$$
(5)
We have chosen three weighting schemes: the first is proportional to the prediction rate on the training set ($`w_jp_{j,train}`$), the second is the equal weight ($`w_j1`$), and the last one is proportional to the maximum likelihood as obtained from the training set ($`w_j\widehat{L}_j`$). Fig.3-4 shows the behavior of all three weighting schemes on both the training and the testing sets, using either the mean square error or the prediction error, up to 200 genes.

The maximum likelihood weight is equivalent to the Akaike weight ($`w_jexp(\text{AIC}_j/2)`$) \[Parzen, Tanabe, Kitagawa, 1998\] and Bayesian weight ($`w_jexp(\text{BIC}_j/2)`$) \[Raftery, 1995\], since all models being averaged in Eq.(5) have the same number of parameters and same sample size (assuming high-order terms are ignored). Using Bayesian weight in model averaging is essentially a derivation of the posterior predictive distribution \[Gelman, et al, 1995\]: “posterior” because data in a training set is used, and “predictive” because unknown new data in the testing set is to be predicted.

For our data set, this weighting scheme is nevertheless uninteresting: the mean square error only decreases slightly with the number of models used, while the error in prediction rate is unchanged. The reason for this is very simple: the best model (the single-gene logistic regression using the g4847) discriminates the training set perfectly. Its likelihood is much higher than any other models. As a result, the weight of other models is negligible, and the model averaging essentially remains as one model. The equal weight can be potentially incorrect since models that do not fit the data should not contribute as equally as models that fit the data better. In our data set, however, the equal weight scheme is actually similar to the weighting scheme that uses the prediction rate, because many genes (up to 200) exhibit similar prediction rates on the training set.

It is difficult from Fig.3-4 to determine a cutoff on the number of models (number of genes) to be included. Fig.3-4, however, clearly shows that a lower number of models (genes) is typically better (except the maximum likelihood weight, which is insensitive to the number of genes included). Fig.3, in particular, indicates that it is possible to perform better (in term of mean square error) on the testing set than using just one model if the number of models is less than around 25. Of course, this result is obtained from the specific testing set we have at hand. A more conclusive result may require a combining of the training and testing set, or requires more sample points than are currently available.

Due to the weight in a model averaging, the apparent number of terms (models, genes) included does not reflect on the true number of genes involved. For this reason, we may introduce a quantity called “effective number of genes (terms, models)”. The leading term in a model averaging contributes a number of 1 to this quantity, but the contribution from the term j is equal to $`w_j/w_1`$. For example, with the maximum likelihood weight, the relative weights of the first ten terms are 1, 0.031, 0.0043, 0.0034, 0.0032, 0.0024, 0.0022, 0.0014, 0.0014, and 0.0006. The effective number of terms when the top ten genes are included is 1.05, a far less number than 10.

The discrimination used in \[Golub et al, 1999\] is a model averaging instead of a model selection. The number of genes used is 50. To quote from \[Golub, et al, 1999\], “the number was well within the total number of genes strongly correlated with the class distinction, seemed likely to be large enough to be robust against noise, and was small enough to be readily applied in a clinical setting.” Our Fig.3-4 shows that although the prediction rate on the training set stays close to 100% even as the number of models (genes) is increased, the prediction rate on the testing set decreases with more genes. The only exception is the maximum likelihood weight, where the prediction rate is almost unchanged due to the dominance of the best gene. It seems that we should not increase the number of models (genes) in model averaging arbitrarily.

In conclusion, due to the small sample size and the presence of strong predictors, we believe the number of genes used in a discriminant analysis in this data set can be much smaller than 50. Although we can not give a definite answer as to the exact number of genes to be used, one proposal is to use only one or two genes, and other exploratory data analyses indicate an upper limit of 10-20 genes. Similar analysis of other data sets for cancer classifications using microarray will be discussed in \[Li, et al. 2001\] and Zipf’s law in these data sets will be discussed in \[Li, 2001\].

## Acknowledgements

W.Li’s work was supported by NIH grant K01HG00024 and Y. Yang’s work was supported by the grant MH44292.

## References

H Akaike (1974), “A new look at the statistical model identification”, IEEE Transactions on Automatic Control, 19:716-723.

KP Burnham, DR Anderson (1998), Model Selection and Inference (Springer).

A Gelman, JB Carlin, HS Stern, DB Rubin (1995), Bayesian Data Analysis (Chapman & Hall).

S Geisser (1993), Predictive Inference: An Introduction (Chapman & Hall).

TR Golub, DK Slonim, P Tamayo, C Huard, M Gaasenbeek, JP Mesirov, H Coller, ML Loh, JR Downing, MA Caligiuri, CD Bloomfield, ES Lander (1999), “Molecular classification of cancer: class discovery and class prediction by gene expression monitoring”, Science, 286:531-537.

F Haghighi, P Banerjee, W Li (1999), “Application of artificial neural networks in whole-genome analysis of complex diseases” (meeting abstract), Cold Spring Harbor Meeting on Genome Sequencing & Biology, page 75.

W Li (1997-2001), An online resource on Zipf’s law (URL: http://linkage.rockefeller.edu/wli/zipf/).

W Li (2001), “Zipf’s law in importance of genes for cancer classification using microarray data”, submitted.

W Li, A Sherriff, X Liu (2000), “Assessing risk factors of complex diseases by Akaike information criterion and Bayesian information criterion” (meeting abstract), American Journal of Human Genetics, 67 (supp 2), page 222.

W Li, D Nyholt (2001), “Marker selection by AIC/BIC”, Genetic Epidemiology, in press.

E Parzen, K Tanabe, G Kitagawa (1998), Selected Papers of Hirotugu Akaike (Springer).

AE Raftery (1995), “Bayesian model selection in social research”, in Sociological Methodology, ed. PV Marsden (Blackwells), pages 185-195.

BD Ripley (1996), Pattern Recognition and Neural Networks (Cambridge University Press).

G Schwarz (1976), “Estimating the dimension of a model”, Annals of Statistics, 6:461-464.

GF Zipf (1949), Human Behavior and the Principle of Least Effect (Addison-Wesley).


text:
# Nonlocal effects in the shot noise of diffusive superconductor - normal-metal systems

## Abstract

A cross-shaped diffusive system with two superconducting and two normal electrodes is considered. A voltage $`eV<\mathrm{\Delta }`$ is applied between the normal leads. Even in the absence of average current through the superconducting electrodes their presence increases the shot noise at the normal electrodes and doubles it in the case of a strong coupling to the superconductors. The nonequilibrium noise at the superconducting electrodes remains finite even in the case of a vanishingly small transport current due to the absence of energy transfer into the superconductors. This noise is suppressed by electron-electron scattering at sufficiently high voltages.

PACS numbers: 72.70.+m, 74.40+k, 74.50+r

Recently, the noise properties of hybrid systems involving superconducting (S) and normal (N) metals became a subject of intensive studies. A key effect in these properties is the Andreev reflection, in which electrons incident from the normal metal are reflected from the NS interface as holes. In particular, it was found that in the zero-voltage limit, the shot noise in diffusive NS contacts with phase-coherent transport is doubled with respect to its value in a normal contact with the same resistance. This doubling was interpreted as an effective doubling of electron charge and has been experimentally confirmed in a number of papers.

Quite recently, it was shown that the doubled shot noise in diffusive NS contacts survives at finite voltages of the order of the energy gap. Moreover, this noise does not require a phase coherence and may be described in terms of a semiclassical Boltzmann - Langevin equation. In this approach, the increased noise in NS systems is due to an excess heating of electron gas rather than to the doubling of the effective charge.

Along with studying the shot noise in two-terminal systems, a considerable attention received the noise in multiterminal structures. It was found that the noise of normal-metal diffusive structures with more than two electrodes may exhibit exchange effects and nonlocal effects. The latter imply that the noise in the system is affected by a presence of open contacts even in the absence of average current flow through them.

Several authors also calculated current correlations in single-channel multiterminal NS systems. In particular, it was found that the noise measured at the two normal electrodes in a four-terminal system depends on the phase difference between the two superconducting electrodes. A current noise was also calculated at a tip placed on a multimode quantum-coherent NS structure.

In this paper, we consider nonlocal effects in multiterminal mesoscopic diffusive SN systems using the semiclassical description of noise. In particular, we report on a semiclassical ”proximity” effect in the noise where the current noise in a normal conductor is increased by its contact with a superconductor. Another finding is that under certain conditions, a large noise may be induced in an SNS structure even by a small transverse transport current.

Consider a cross-shaped contact with two superconducting and two normal-metal electrodes (see Fig. 1). The normal-metal electrodes are kept at potentials $`\pm V/2`$, and the superconducting electrodes are kept at zero potential. It is also assumed that the cross is symmetric, i.e. it consists of two identical resistances $`R`$ connecting the crossing point with the normal electrodes and two identical resistances $`r`$ connecting it with the superconductors. Each resistor presents a long and narrow diffusive wire with the Thouless energy much smaller than the energy gap of the superconductor. The resistance of the crossing point is assumed to be negligible. The applied voltage $`eV`$

is assumed to be larger than the Thouless energy but smaller than the energy gap of the superconductors. Note that in this geometry, there is now electrical current flow through the superconducting ends of the cross.

Under the above conditions, the kinetics of fluctuations may be described using a semiclassical Langevin equation

$$\delta 𝐣=D\frac{}{𝐫}\delta \rho \sigma \frac{}{𝐫}\delta \varphi +\delta 𝐣^{ext},$$
(1)
where $`D`$ is the diffusion coefficient, $`\sigma `$ is the electric conductivity, $`\delta \rho (𝐫)`$ is the charge-density fluctuation, $`\delta \varphi (𝐫)`$ is the local fluctuation of the electric potential, and the correlator of extraneous currents $`\delta 𝐣^{ext}`$ is given by

$$\delta j_\alpha ^{ext}(𝐫_1)\delta j_\beta ^{ext}(𝐫_2)_\omega =4\sigma \delta _{\alpha \beta }\delta (𝐫_1𝐫_2)T_N(𝐫_1),$$

$$T_N(𝐫)=𝑑\epsilon f(\epsilon ,𝐫)[1f(\epsilon ,𝐫)].$$
(2)
Equation (1) may be integrated along the length of each arms of the cross, which gives for the fluctuations of currents flowing into each of the four electrodes

$$\delta I_i=\frac{1}{R_i}\left(\delta \varphi ^{}+\frac{1}{e^2N_F}\delta \rho ^{}\right)+\delta I_i^{ext},$$
(3)
where $`R_i`$ is the resistance of the corresponding arm, $`\delta \varphi ^{}`$ and $`\delta \rho ^{}`$ are the fluctuations of electrical potential and charge density at the crossing point, and the correlator of extraneous currents $`\delta I_i^{ext}`$ equals

$$\delta I_i^{ext}(t_1)\delta I_j^{ext}(t_2)_\omega =\delta _{ij}\frac{4T_i}{R_i},$$
(4)
where $`T_i`$ is obtained by averaging (2) over the length of the corresponding arm. The system of equations (3) combined with the current-conservation condition at the crossing point

$$\underset{i}{}\delta I_i=0$$

is easily solved giving

$$\delta I_i=\delta I_i^{ext}\frac{1}{R_i}\underset{j}{}\delta I_j^{ext}/\underset{j}{}\frac{1}{R_j},$$
(5)
and hence the cross-correlated spectral density is

$$S_{ij}\delta I_i\delta I_j_\omega =\delta _{ij}\frac{4T_i}{R_i}\frac{4}{R_iR_j}\frac{T_i+T_j}{\underset{k}{}1/R_k}$$

$$+\frac{4}{R_iR_j}\underset{k}{}\frac{T_k}{R_k}/\left(\underset{k}{}\frac{1}{R_k}\right)^2.$$
(6)
In what follows, we restrict ourselves to relatively small voltages $`eV<2\mathrm{\Delta }`$ and zero temperatures. In each arm of the cross, the average distribution function $`f(\epsilon ,x_i)`$ obeys the standard diffusion equation $`^2f=0`$. Introducing the distribution function $`f^{}(\epsilon )`$ at the crossing point, one can write down the distribution in each arm in the form

$$f(\epsilon ,x_i)=\left(1\frac{x_i}{L_i}\right)f^{}(\epsilon )+\frac{x_i}{L_i}f_i(\epsilon ),$$
(7)
where $`L_i`$ is the length of the arm and $`f_i`$ is the distribution at the end of it. The distribution function $`f^{}`$ should be determined from the balance of the diffusion fluxes at the crossing point for any energy

$$\underset{i}{}\frac{f(\epsilon ,x_i)}{x_i}|_{x_i=0}=0.$$
(8)
The distribution functions at the normal ends of the cross $`f_1=f_0(\epsilon eV/2)`$ and $`f_3=f_0(\epsilon +eV/2)`$ are just the equilibrium Fermi functions shifted in energy by $`\pm eV/2`$. Owing to the Andreev reflections from the interfaces with the superconductors, the distribution functions at these interfaces are related to the distribution function at the crossing point by a very simple expression

$$f_2(\epsilon )=f_4(\epsilon )=\frac{1}{2}[1+f^{}(\epsilon )f^{}(\epsilon )].$$
(9)
In the symmetric case where $`R_1=R_3=R`$ and $`R_2=R_4=r`$, it is easily obtained that

$$f^{}(\epsilon )=f_2=f_4=\frac{1}{2}(f_1+f_3).$$
(10)
Substituting Eq. (10) into Eq. (7), calculating the corresponding $`T_N(x_i)`$ by means of (2) and averaging them over the corresponding segments, one easily obtains that $`T_1=T_3=eV/6`$ and $`T_2=T_4=eV/4`$. From this, one readily obtains the cross-correlated spectral densities by means of (6). Taking into account that the average current flowing between the normal electrodes is $`I=eV/2R`$, the expression for the noise at the normal ends may be written in the form

$$S_{11}=S_{33}=\frac{eI}{3}\frac{4R^2+7Rr+2r^2}{(R+r)^2},$$
(11)
The cross-correlated noise at the normal ends is

$$S_{13}=\frac{eI}{3}\frac{r(R+2r)}{(R+r)^2}.$$
(12)
The spectral densities of noise at the superconducting ends of the cross are given by

$$S_{22}=S_{44}=\frac{eV}{6r}\frac{3R^2+8Rr+6r^2}{(R+r)^2}$$
(13)
and

$$S_{24}=\frac{eV}{6r}\frac{R(3R+4r)}{(R+r)^2}.$$
(14)
To complete the description, we also present the cross-correlated spectral density at the normal and superconducting ends

$$S_{12}=\frac{eV}{6}\frac{3R+2r}{(R+r)^2}.$$
(15)
Note that all the cross-correlated spectral densities (12), (14), and (15) are always negative.

Consider now the most intresting case of a strong coupling to the superconductors where $`rR`$. It is seen from Eqs. (11) and (12) that the noise measured at one of the normal ends of the cross is doubled with respect to a two-terminal normal-metal system, as it takes place in diffusive contacts where the transport current flows through an SN interface. However the physics of this effect is different because the cross-correlated spectral density of noise at the normal ends (12) tends to zero.

The reason for the increased shot noise is that the ”superconducting” arms, which do not carry any transport current, are open yet for the current fluctuations and act as noise generators that supply additional electric fluctuations into the system. In the case of normal-metal electrodes 2 and 4, the noise should be also slightly increased, but the increment would reach its 25% maximum ar $`r=R`$ (see Fig. 2). At strong coupling, such cross would just break down into two independent contacts, each with the resistances $`R`$ and a voltage drop $`V/2`$. Hence the noise at electrodes 1 and 3 would be just $`2eI/3`$. However the NS interfaces, while being transparent for the fluctuations of electric current, forbid the energy transport into the superconductors and hence do not allow the distribution function at the crossing point to assume the equilibrium shape. Clearly, the doubling of shot noise should take place also in a three-terminal structure with only one superconducting electrode attached in the middle of the normal conductor.

It is noteworthy that the noise in the ”superconducting” arms remains finite even if $`R\mathrm{}`$ and the transport current through the system vanishes. In this case, $`S_{22}=S_{44}=S_{24}=eV/2r`$. This is in contrast to a purely normal system, where $`S_{22}=(2/3)eV/(R+r)`$ and tends to zero as either $`r`$ or $`R`$ becomes infinitely large. The reason for this is that the electron gas in this case is confined between two interfaces with superconductors, which hinder heat transfer from it. Hence it may be strongly heated even with a small current, much like as in the case of two-terminal diffusive SNS contacts it is strongly heated even by a small voltage.

As $`R0`$, the spectral densities $`S_{22}`$ and $`S_{44}`$ increase to $`eV/r`$, while the cross-correlated spectral density $`S_{24}`$ tends to $`2eV/3r`$, which implies that the current fluctuations at the opposite superconducting ends become only partially correlated.

Consider now the effects of electron-electron scattering on the spectral density $`S_{22}`$ in the case of a strong coupling. In the case of two-terminal SNS contacts, this type of scattering was shown to suppress the nonequilibrium noise at low voltages because the quasiparticles confined between the two NS interfaces may be outscattered from the subgap energy range and escape into the superconducting electrodes thus transferring energy and effectively cooling the electron gas.

Consider the case where the length of a ”normal” arm $`L_R`$ is much larger than the electron-electron scattering length. In this case, the electron gas may be described by a local effective temperature, which depends only on the coordinate $`x`$ in the direction of average current flow and is constant in the ”superconducting” arms at $`x=0`$. This effective temperature obeys a heat-balance equation

$$\frac{\pi ^2}{6}\frac{d^2}{dx^2}\left(T_e^2\right)=\left(\frac{eV}{2L_R}\right)^2+\delta (x)\frac{2}{L_R\epsilon _T}\frac{R}{r}J,$$
(16)
where $`\epsilon _T`$ is the Thouless energy of a ”superconducting” arm and

$$J=\alpha _{ee}\frac{\mathrm{\Delta }^3T^{}}{\epsilon _F}\mathrm{exp}\left(\frac{\mathrm{\Delta }}{T^{}}\right)$$

is the density of flux of energy carried by electrons and holes outscattered from the subgap region by electron-electron collisions, which depends on the effective temperature in the ”superconducting” arms $`T^{}=T_e(0)`$. Making use of the boundary conditions $`T_e(L_R)=T_e(L_R)=0`$, one easily obtains a closed equation for $`T^{}`$ in the form

$$(T^{})^2=\frac{3}{4\pi ^2}(eV)^2\frac{6}{\pi ^2}\alpha _{ee}\frac{R}{r}\frac{\mathrm{\Delta }^3T^{}}{\epsilon _F\epsilon _T}\mathrm{exp}\left(\frac{\mathrm{\Delta }}{T^{}}\right).$$
(17)
From this equation, it is readily seen that in the case of large values of the prefactor $`\lambda =\alpha _{ee}(R/r)(\mathrm{\Delta }^2/\epsilon _F\epsilon _T)`$ the effective temperature $`T^{}`$ and the spectral density $`S_{22}=2T^{}/r`$ become suppressed at voltages $`V>\mathrm{\Delta }/(e\mathrm{ln}\lambda )`$. Unlike the case of two-terminal SNS contacts, the suppression takes place at high voltages rather than at small ones.

In summary, we semiclassically considered nonequilibrium noise in diffusive multiterminal NS structures and found that nonlocal effects in them are by far more pronounced than in purely normal systems with the same geometry. In particular, they allow an observation of a doubled longitudinal noise and of a giant transverse nonequilibrium noise with respect to the direction of transport current.

The cross-shaped SN structures considered above are easily fabricated and hence the theoretical conclusions about the noise may be easily tested by experimentalists. An advantage of this system is that there is no voltage drop between the superconducting electrodes and the effects of electron heating are not obscured by an onset of ac Josephson effect, as it takes place for two-terminal SNS structures.

I am grateful to W. Belzig for attracting my attention to multiterminal geometries, to C. Strunk for a discussion of experimental results, and to M. Büttiker for a discussion.

This work was supported by the Swiss National Science Foundation and Russian Foundation for Basic Research, Grant No. 01-02-17220.


text:
# Is the active gravitational mass of a charged body distance-dependent?

(12 April 2001)

## Abstract

It appears to follow from the Reissner-Nordstrøm solution of Einstein’s equations that the charge of a body reduces its gravitational field. In a recent note Hushwater offered an explanation of this apparent paradox. His explanation, however, raises more questions than solves since it implies that the active gravitational mass of a charged body is distance-dependent and therefore is not equal to its inertial mass.

As discussed by Hushwater  the Reissner-Nortstrøm solution of Einstein’s equations for a charged body of mass $`M`$ and charge $`Q`$

$$ds^2=\left(1\frac{2GM}{c^2r}+\frac{GQ^2}{c^4r^2}\right)c^2dt^2+\left(1\frac{2GM}{c^2r}+\frac{GQ^2}{c^4r^2}\right)^1dr^2+r^2\left(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2\right)$$
(1)
may lead one to the conclusion that the body’s gravitational filed is reduced by its charge. This becomes obvious if the Newtonian limit of general relativity is considered. In that limit the metric tensor of curved spacetime $`g_{\alpha \beta }`$ can be represented by the metric tensor of flat spacetime $`\eta _{\alpha \beta }`$ and another ”perturbation” tensor $`h_{\alpha \beta }`$ whose components are much less than unity (since they are proportional to $`c^n`$ where $`n2`$)

$$g_{\alpha \beta }=\eta _{\alpha \beta }+h_{\alpha \beta }.$$

In the Newtonian limit of Reissner-Nortstrøm metric

$$g_{00}=1+h_{00}$$

where

$$h_{00}=\frac{2GM}{c^2r}+\frac{GQ^2}{c^4r^2}.$$

The equation of motion of a non-relativistically moving test particle in terms of $`h_{00}`$ is then

$$\frac{d^2𝐫}{dt^2}=\frac{c^2}{2}h_{00}$$

or

$$\frac{d^2𝐫}{dt^2}=\left(\frac{GM}{r}\frac{GQ^2}{2c^2r^2}\right).$$
(2)
As seen from (2) it appears that the charge $`Q`$ of the body reduces its gravitational field since the second term in the parentheses is subtracted from the gravitational potential $`GM/r`$. The problem with such a conclusion is that it contradicts the very foundations of relativity according to which the electric field of the body must, in fact, increase its gravitational field since the electric field possesses energy and therefore mass.

Hushwater claims to have resolved this apparent paradox by making use of the concept ”total mass inside a radius $`r`$” felt by a test particle at a distance $`r`$ from the body’s center

$$M\left(r\right)=M\frac{1}{8\pi c^2}|𝐄|^2d^3x$$
(3)
where $`M`$ is the total mass of the charged body consisting of its ordinary mass and the *whole* electromagnetic mass that corresponds to the energy of the body’s electric field $`𝐄`$. The integration in the second term in (3) is taken over the space outside a sphere $`S(r)`$ of radius $`r`$ and therefore that term is the *part* of the electromagnetic mass that is stored in the body’s electric field occupying the space outside $`S(r)`$. In such a way the mass $`M\left(r\right)`$ comprises the ordinary mass of the body and only that part of its electromagnetic mass that corresponds to the energy of the electric field inside the sphere $`S(r)`$. As

$$𝐄=\frac{Q}{r}=\frac{Q}{r^2}𝐧$$

where $`𝐧=𝐫/r`$, the second term in (3) is equal to $`Q^2/2c^2r`$. Therefore for (3) one obtains

$$M\left(r\right)=M\frac{Q^2}{2c^2r}.$$
(4)
The classical equation of motion of a particle in a gravitational potential $`GM\left(r\right)/r`$ caused by the mass $`M\left(r\right)`$ is

$$\frac{d^2𝐫}{dt^2}=\frac{GM\left(r\right)}{r}.$$
(5)
The substitution of (4) in (5) gives (2). This result is regarded by Hushwater as a proof that the apparent paradox disappears if both the ordinary and electromagnetic mass of a charged body are taken into account and the mass $`M\left(r\right)`$ inside a sphere $`S(r)`$ is used.

This resolution of the paradox, however, comes at too high a price since it is based on two implicit assumptions none of which seems to be correct.

1. The mass $`M\left(r\right)`$ is implicitly regarded as the *active gravitational mass* of the charged body, which is felt by a test particle placed at a distance $`r`$ from the body’s center. And indeed, as seen from (5) it is the mass $`M\left(r\right)`$ that gives rise to the gravitational potential $`GM\left(r\right)/r`$. This means, however, that the active gravitational mass of a charged body is *distance-dependent*. Leaving aside the question of what that may mean, it is obvious that a distance-dependent active gravitational mass is not equal to the body’s inertial mass, which consists of its ordinary mass and its *entire* electromagnetic mass; $`M\left(r\right)`$ coincides with the inertial mass of the charged body only when $`r\mathrm{}`$. An assumption that the active gravitational mass of a charged particle is not equal to its inertial mass is not justified since there exists no evidence that the equivalence principle is violated in the case of charged particles .

2. Regarding the mass $`M\left(r\right)`$ as the source of the gravitational potential $`GM\left(r\right)/r`$ in (5) implies that $`M\left(r\right)`$ defines some metric. To find that metric one can substitute the expression for $`M\left(r\right)`$

$$M\frac{Q^2}{2c^2r}=M\left(r\right)$$

from (4) in (1). The result is

$$ds^2=\left(1\frac{2GM\left(r\right)}{rc^2}\right)c^2dt^2+\left(1\frac{2GM\left(r\right)}{rc^2}\right)^1dr^2+r^2\left(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2\right)$$
(6)
which is the Schwarzschild metric in the case of a body of active gravitational mass $`M\left(r\right)`$. In other words, the second implicit assumption is that the Reissner-Nortstrøm metric (1) can be directly obtained from the Schwarzschild metric (6) if the expression for the mass $`M\left(r\right)`$ is substituted in (6). However, it is obviously incorrect to use $`M\left(r\right)`$ in the Schwarzschild solution since it is a vacuum solution $`\left(T_{\alpha \beta }=0\right)`$ whereas in the case of a charged body $`T_{\alpha \beta }0`$.

If the equivalence principle strictly holds for charged particles, then the paradox that the charge of a body reduces its gravitational field according to the Reissner-Nortstrøm solution remains. And if no other explanation of that paradox is found the only way out of this situation seems to be to assume that the active gravitational mass of a charged body is indeed distance-dependent and therefore is not equal to its inertial mass . This would mean that the Reissner-Nortstrøm solution of Einstein’s equations does follow from the Schwarzschild solution if the expression for the mass $`M(r)`$ is substituted in (6). Then a justification for such an illegal at least at first glance operation might be the following. A test particle at a distance $`r`$ from the charged body’s center does not feel the effect of the electromagnetic mass corresponding to the body’s electric field outside the sphere $`S(r)`$ since it cancels out exactly. In this sense the space outside the sphere $`S(r)`$ will appear ”empty” to the test particle.