text:
# 1 Introduction

## 1 Introduction

The mechanosensitive channel (MscL) is a compelling example of the interaction between a protein and the surrounding bilayer membrane. The channel is gated mechanically by applied tension and is believed to function as an emergency relief valve in bacteria . MscL is a member of a growing class of proteins which have been determined to be mechanosensitive , . The dependence of the conductance on applied tension has been studied extensively in patch clamp experiments , , . In terms of the observed conductance, these studies have revealed that the channel is very nearly a two state system. MscL spends the vast majority of its life in either a closed state (C) or an open state (O) characterized by a discrete conductance. When the bilayer tension is small, the protein is exclusively in the closed configuration. As the tension grows, the open state becomes ever more prevalent, until it dominates at high tension. The simplest structural interpretation of this conductance data is to assume that each discrete conductance corresponds to a well defined channel conformation. This assumption seems to be compatible with the conductance data111The area of the states seems roughly independent of applied tension .. Patch clamp experiments have also revealed that there are at least three additional discrete, intermediate conductance levels  suggesting three additional short lived substates (S1-S3). Rees and coworkers  have solved the structure for one conformation using X-ray crystallography. This state appears to be the closed state , . Perozo et al.  have trapped MscL in the open configuration and used electron paramagnetic resonance spectroscopy (EPR) and site-directed spin labeling (SDSL) to deduce its geometry. Sukharev et al.  have also proposed an open state conformation based on structural considerations.

The conformational landscape of the MscL channel is extremely complex, depending on huge number of microscopic degrees of freedom which are analytically intractable. Even from the standpoint of numerical calculations, this number is still very large . As an alternative to a detailed microscopic picture of MscL, we consider a simplified free energy function where we divide the free energy of the system into two contributions, namely,

$$G=G_P+G_{},$$
(1)
where $`G_P`$ is the free energy associated with the conformation of the protein and $`G_{}`$ is the deformation free energy from the bulk of the bilayer . In general, these two terms are coupled. The conformation of the protein depends on the forces applied by the bilayer. The bilayer deformation is induced by the external geometry of the protein. We denote this external geometry with a state vector, $`X`$, which captures the radius of the channel as well as its orientation relative to the surrounding bilayer as described in more detail below. We calculate the induced bilayer deformation energy, $`G_{}(X)`$, by minimizing the free energy of the bilayer and solving the resulting boundary value problem using an analytic model developed for the study of bilayer mechanics  and protein-bilayer interactions , , , . We then apply asymptotic approximations to the exact solutions of this model for cylindrically symmetric inclusions, permitting all of the results to be expressed, estimated, and understood with simple scaling relations. The advantage of this model is that it permits us to characterize the protein-bilayer system in a way that is at once analytically tractable and predictive. By understanding the consequences of the simplest models, we develop a framework in which to understand the richer dynamics of the real system. There is a wealth of useful, physical intuition to be gleaned from this model relating to both the function of MscL and that of mechanosensitive transmembrane proteins more generally. In a forthcoming paper, we will show that the mechanics of the bilayer must play an integral role in mechanotransduction and channel function. Specifically, we will present detailed analytic estimates of the free energy generated by bilayer deformation induced by the channel and show that these free energies are of the same order as the free energy differences measured by Sukharev et al. . These analytic calculations reveal a compelling and intuitive model for the gating of the MscL channel which is the subject of this current paper. The competition between hydrophobic mismatch and applied tension, in the presence of radial constraints, generates a bistable system that is implicitly a mechanosensitive channel. Furthermore, this simple model provides a picture which is both qualitatively and quantitatively consonant with the measured dependence of the free energy on acyl chain length as observed by Perozo et al. . In addition, these results may also explain the stabilization of the open state by spontaneous curvature inducing lysophospholipids observed by Perozo et al. , although more experiments are required to check the consistency of this proposal.

## 2 The energy landscape of the bilayer

In the calculations considered here, the geometry of the protein, characterized by the conformational state vector $`X`$, is described by three geometrical parameters

$$X=(R,W,H^{}),$$
(2)
where $`R`$ is the radius of the channel, $`W`$ is the hydrophobic thickness, and $`H^{}`$ is the mid-plane slope. See figure 1 for details. Although we have parameterized the conformation space of the protein with these three parameters, in this paper, we will focus on the radial dependence alone, claiming that even in this reduced description, the model provides a rich variety of predictions which are compatible with previous observation and suggest new experiments. The radial dependence of the bilayer deformation energy is particularly important for MscL since the radius undergoes a very large change between the open and closed states . The bilayer deformation energy can be written explicitly in terms of the channel radius as

$$G_{}=G_0+f2\pi R\alpha \pi R^2,$$
(3)
where $`G_0`$ and $`f`$ do not explicitly depend on $`R`$ and $`\alpha `$ is the applied tension which triggers channel gating. $`G_0`$ is a radially independent contribution to the deformation energy which is a function of the other geometrical parameters of the protein. Its importance in gating the channel is most likely secondary since it is independent of $`R`$ and it will be ignored in the remainder of the discussion. The dependence of bilayer deformation energy on applied tension can be explained intuitively . The free energy contribution for a small change in the channel area due to the applied tension can be written $`\alpha dA`$ which is the two dimensional analogue of the $`PdV`$ term for a gas in three dimensions. At high enough applied tension, the state with the largest inclusion area will have the lowest free energy.

The line tension, $`f`$, contributes an energy proportional to the circumference. Line tensions arise very naturally since the interface area of the inclusion is proportional to the circumference. In what follows, we will discuss the two dominant contributions to this line tension: thickness deformation , ,  and spontaneous curvature . The thickness deformation free energy is induced by the mismatch between the the equilibrium thickness of the bilayer and the hydrophobic thickness of the protein. The importance of this hydrophobic mismatch in the function of transmembrane proteins has already been established . The bilayer deforms locally to reduce the mismatch with the protein as shown in figure 1. Symbolically the thickness deformation energy is ,

$$G_U=f_U2\pi R=\frac{1}{2}𝒦U^22\pi R,$$
(4)
where $`𝒦=2\times 10^2kT\mathrm{\AA }^3`$ is an effective elastic modulus defined in the appendix which is roughly independent of acyl chain length222See the appendix for an brief discussion of scaling of the elastic constants. and $`U`$ is half the hydrophobic mismatch as defined in figure 1. Naturally the energetic penalty for this deformation is proportional to the mismatch squared since the minimum energy state corresponds to zero mismatch. The area of that part of the bilayer which is deformed is roughly equal to the circumference of the channel times an elastic decay length. As a result, the contribution of thickness deformation to the total free energy budget scales with the radial dimension of the channel. We also note that the thickness deformation free energy is always positive.

In contrast, the free energy induced by spontaneous curvature can be either negative or positive. Physically, this free energy comes from locally relieving or increasing the curvature stress generated by lipids or surfactants which induce spontaneous curvature , , . Again the radial dependence of this free energy will be linear since the effect is localized around the interface. Since the leaflets of the bilayer can be doped independently , the spontaneous curvatures of the top and bottom leaflets, $`C_\pm `$, can be different. It is convenient to work in terms of the composite spontaneous curvature of the bilayer, $`C\frac{1}{2}(C_+C_{})`$. The contribution to the deformation energy arising from spontaneous curvature is given by

$$G_C=f_C2\pi R=K_BCH^{}2\pi R,$$
(5)
where $`H^{}`$ is the mid-plane slope and $`K_B=20(a/20\mathrm{\AA })^3kT`$, is the bending modulus which roughly scales as the third power of the bilayer thickness333See the appendix for an brief discussion of scaling of the elastic constants.. We will discuss these results in more detail elsewhere. Notice that if $`C`$ and $`H^{}`$ have opposite signs, the deformation energy and the corresponding line tension, $`f_C`$, will be negative. We note that the elastic theory of membrane deformations associated with proteins like MscL permit other terms (such as mid-plane deformation, for example) which can be treated within the same framework and which give rise to the same radial dependence as that described here. However, for the purposes of characterizing the energetics of MscL, these other terms are less important than the two considered here.

Typically, in the absence of large spontaneous curvature, the line tension, $`f`$, will be dominated by the mismatch and will be positive. A potential of the form described by equation 3 is depicted schematically in figure 2. In this figure, we have implied that there are steric constraints for the range of radii accessible to the protein. Assuming that there is a lower bound on the radius of the inclusion is very natural. It can be understood as the radius below which the residues begin to overlap. This steric constraint will generate a hard wall in the protein conformation energy, forbidding lower radii. Similar but slightly more elaborate arguments can made for an upper bound. The bilayer deformation energy generates a barrier between small radius and large radius states. The location of the peak of this barrier is the turning point

$$R_{}\frac{f}{\alpha }.$$
(6)
At small tension, the turning point is very large and is irrelevant since it occurs at a radius not attainable by the channel due to the steric constraints, but as the tension increases the position of the turning point decreases. This behavior is reminiscent of the competition between surface tension and energy density for nucleation processes which give rise to a similar barrier (e.g. ).

### 2.1 Induced tension and bilayer induced stabilization

Although the conformational landscape of the MscL channel is certainly very complicated, there is an intriguing possibility that the channel harnesses the elastic properties of the bilayer which quite naturally provide the properties we desire in a mechanosensitive channel: a stable closed state at low tension and a stable open state at high tension. In effect, we will treat the bilayer deformation energy as an external potential with respect to the conformational energy landscape of the protein. The physical effects of the radial dependence of the bilayer deformation energy on the inclusion conformation can be recast in a more intuitive form by appealing to the induced tension which accounts not only for the applied far field tension, but also for induced tension terms due to bilayer deformation. The applied tension is not the whole story! The generalized forces are obtained by differentiating the bilayer deformation energy with respect to bilayer excursions. The net tension induced by the bilayer on the inclusion interface is

$$\alpha _\mathrm{\Sigma }=\alpha \frac{f}{R},$$
(7)
where we have denoted the net tension $`\alpha _\mathrm{\Sigma }`$ since we have already used $`\alpha `$ to denote the applied tension. For radii smaller than the turning point, $`R_{}`$, the bilayer deformation energy is an increasing function of radius and therefore the net tension is negative and acts to compress the channel. For radii larger than that at the turning point, the bilayer deformation energy is a decreasing function of radius and the net tension is positive and acts to expand the channel. The combination of these constraints and the bilayer deformation energy lead to a bistable system where the closed and open states correspond to the constraint-induced radial minimum and maximum, respectively. Recall that the net tension on the closed state will be compressive as long as its radius is smaller than that at the turning point, namely,

$$R_C<R_{}=\frac{f}{\alpha }.$$
(8)
This inequality defines the range of applied tension over which the closed state is stabilized by the bilayer deformation energy. The net tension on the open state will be expansive as long as its radius is greater than that at the turning point:

$$R_O>R_{}=\frac{f}{\alpha }.$$
(9)
This inequality defines the range of applied tension over which the open state is stabilized by the bilayer deformation energy. There is an intermediate range of tensions for which both states are stabilized by the bilayer,

$$\frac{f}{R_O}<\alpha <\frac{f}{R_C}.$$
(10)
The bilayer deformation energy naturally destabilizes the open state for applied tension below this range while stabilizing the closed state for applied tensions up to the limit of this range. Both effects help to prevent the channel from leaking at low applied tension. This bistability is precisely the desirable behavior for a mechanosensitive channel designed to relieve internal pressure and yet surprisingly little is required from the protein conformational landscape, $`G_P`$, except for steric constraints which arise very naturally. In figure 3, we have depicted a mechanical analogue of this bistable system. In figure 4, we have depicted the way in which MscL mimics this mechanical analogue by using the sum of a schematic protein energy and the bilayer deformation energy to form energy minima corresponding to the open and closed states.

Our general discussion of the role of the bilayer enables us to make some rather general observations about the nature of the substates. In order to generate a substate, we assume that there is more than one gating transition in the protein conformational energy. One gating transition would correspond to a closed to open transition. An additional gating transition allows three conductance states. We will assume that these transitions are themselves bistable in nature since the conductance data would seem to imply the lifetimes of the transition states are very short compared to the conductance states . In other words, the conformational gating transition occurs near a local maximum in the conformational free energy like that depicted in figure 3. If we add two such transitions to $`G_P`$, we generate a substate of intermediate radius between these two transition state radii. A schematic example of this is illustrated in figure 4. Sukharev et al.  have shown that all the substates are short lived and have estimated the areas of each state based on the tension dependence of their free energies444There is now evidence for additional substates .. Specifically, they have shown that the radii of the substates lie between the open and closed state radii. If the bilayer deformation dominates the free energy of the states, the ephemeral nature of the substates is a natural consequence of their intermediate radii. The compressive tension due to the mismatch stabilizes the state of lowest radius at low applied tension. At high applied tension the bilayer stabilizes the state with highest radius. All the states with intermediate radii are never stabilized by the bilayer and are therefore short lived. Our deceptively simple mismatch model quite naturally leads to short lived substates at intermediate radii.

In this section we have shown that some very general assumptions about the form of the bilayer deformation energy lead to channel behavior which is mechanosensitive. Although we have suggested a specific source, thickness deformation, for the line tension, $`f`$, this linear radial dependence arises very naturally. For instance, if we assume the bilayer is too stiff to deform and the mismatch causes hydrophobic residues to come into contact with the solvent, a larger positive line tension would be induced. The competition between line tension and applied tension is therefore a quite general phenomena for mechanosensitive transmembrane proteins. We shall demonstrate below that, for MscL, this mechanism appears to explain the free energy dependence on bilayer acyl chain length surprisingly well.

## 3  Results

The patch clamp experiments of Perozo et al.  go beyond the earlier work of Sukharev et al.  by providing experimental values for the free energy difference between the open and closed states for bilayers of several thicknesses. These results can be compared with our predictions. In order to apply our model, we must determine the geometrical parameters of the state vector $`X`$ for the open and closed states and in particular the open and closed radii. The radius of the closed state is known from x-ray crystallography : $`R_C23\mathrm{\AA }`$. Structural studies  and EPR and SDSL  experiments have suggested an open state radius of roughly $`R_O35\mathrm{\AA }`$. In order to estimate the line tension and free energy generated by hydrophobic mismatch, we must determine the hydrophobic thickness $`W`$. (We ignore the difference in hydrophobic thickness between the open and closed states.) In principle, one might have thought this could be deduced from the atomic-level structure of MscL, but in practice, real structures are complicated, often lacking a clear transition from hydrophobic to hydrophilic residues on the interface. However, this width may be deduced from the EPR and SDSL data of Perozo et al. . EPR and SDSL experiments measure inter-subunit proximity and spin-label mobility, respectively . Compressive tension in the bilayer suppresses the fluctuations of the protein, increasing the subunit proximity and reducing the spin label mobility. In the experiments of Perozo et al. , the applied tension is low implying that the net tension is dominated by the line tension, induced by thickness deformation,

$$\alpha _\mathrm{\Sigma }\frac{𝒦U^2}{2R},$$
(11)
in the absence of spontaneous curvature. This tension is compressive and proportional to the mismatch squared. Therefore, when the mismatch is zero, the tension reaches a minimum, implying that mobility and subunit separation should reach a maximum. The EPR and SDSL data of Perozo et al.  may turn over for PC12 bilayers, implying that the mismatch is zero which would imply in turn that

$$W2a_{n=12},$$
(12)
but due to the quadratic dependence on $`U`$, the slope in the vicinity of the turnover is small. Since PC lipids with acyl chain length shorter than $`n=10`$ do not form stable bilayers , it is difficult to extensively check the quadratic dependence on $`U`$. The predicted turnover would be more pronounced for PC bilayers with $`n<10`$. We shall see that this deduced hydrophobic mismatch is compatible with the patch clamp measurements of Perozo et al. .

### 3.1 Hydrophobic mismatch

Perozo et al.  have measured $`\mathrm{\Delta }G_0`$, the free energy difference between the open and closed state at zero tension555The free energy measured by Perozo et al.  is equivalent to the free energy at zero tension modulo several assumptions . for three acyl chain lengths. Using the value we have deduced for $`W`$, we can now calculate the free energy difference between the open and closed states due to bilayer deformation at zero tension, $`\mathrm{\Delta }G_{0,}`$, which is given by the line tension contribution alone:

$$\mathrm{\Delta }G_{0,}=f2\pi \mathrm{\Delta }Rf_U2\pi \mathrm{\Delta }R,$$
(13)
where $`\mathrm{\Delta }R`$ is the difference between the open and closed radii. The theoretical result, $`\mathrm{\Delta }G_{0,}`$, is plotted with the experimental measurements of $`\mathrm{\Delta }G_0`$ in fig 5. The agreement between experiment and theory is embarrassingly good given the cavalier fashion in which we have chosen the geometrical parameters and that we have neglected the protein conformational energy, $`G_P`$, entirely. There is a very important point to be made about these results. Perozo et al.  have measured three data points and our model is quadratic, implying that we could have chosen the parameters of our model to fit the data points perfectly since any three points lie on a parabola, but our parameters have in fact been deduced independently rather than fit, which is why this correspondence with the data is remarkable. This model corresponds to a channel where the free energy difference between the open and closed states is dominated by the bilayer deformation rather than protein conformation. Our model implies that $`\mathrm{\Delta }G_0`$ for PC10, PC12, and PC14 should be very small. Unfortunately these bilayers have proved too weak for patch-clamp measurements of $`\mathrm{\Delta }G_0`$ . Certainly none of these bilayers trap the channel in the open state .

The opening tension is defined as the tension at which the open and closed state probabilities are equal or, analogously, the tension at which the free energies of the open and closed states are equal. The opening tension is

$$\alpha _{1/2}=\frac{f}{\overline{R}}+\frac{\mathrm{\Delta }G_P}{\mathrm{\Delta }A},$$
(14)
where $`\overline{R}\frac{1}{2}(R_C+R_O)`$ is the mean radius, $`\mathrm{\Delta }G_P`$ is the difference in the open and closed state protein conformation energy, and $`\mathrm{\Delta }A`$ is the difference in open and closed state area. When the bilayer deformation energy dominates, the opening tension is determined by the first term alone. Changing the sensitivity of the channel is straightforward from this perspective. Changes in the length of the hydrophobic region of the protein can increase or decrease the opening tension of the channel. For example, MscL channels might be engineered with an expanded hydrophobic region which matches PC14 bilayers. Our mismatch based theory would predict that the free energy versus acyl chain length curve would simply be translated to higher $`n`$ so that the minimum $`\mathrm{\Delta }G_{0,}`$ is realized for a PC14 bilayer. This shift should be measurable, reducing $`\mathrm{\Delta }G_0`$ for PC16, PC18, and PC20 bilayers. The reduction in mismatch may also allow MscL to be reconstituted into PC22 bilayers, allowing an additional data point. The proposed shift should also be measurable in EPR and SDSL measurements of residue proximity and mobility. The maximum mobility and separation should now be centered around $`n=14`$, perhaps permitting a clear measurement of the rise in induced tension for a PC10 bilayer predicted by the quadratic dependence of the line tension on the mismatch.

### 3.2 Spontaneous curvature

Perozo et al.  have proposed that asymmetric bilayer stresses play a central role in MscL gating. They have proposed this model based on patch clamp, EPR, and SDSL experiments showing that spontaneous curvature can induce MscL channel opening. Specifically, Perozo et al. find that MscL reconstituted into PC vesicles with high enough concentrations of asymmetrically incorporated LPC stabilizes the open state of the channel, while MscL reconstituted into PC vesicles with symmetrically incorporated LPC, does not stabilize the open state. Unlike Keller et al. , Perozo et al. have measured neither the spontaneous curvature for the mixed bilayer nor the free energy difference between open and closed states as a function of LPC concentration. In the absence of these quantitative experimental results, it is difficult to make concrete comparisons between our model and the experimental data. For large spontaneous curvature  but a relatively modest complementary mid-plane slope, the free energy difference between the two states due to spontaneous curvature is

$$\mathrm{\Delta }G_C2\pi (R_OR_C)K_BCH^{}=16\left(\frac{20\mathrm{\AA }}{C^1}\right)\left(\frac{H^{}}{0.2}\right)kT,$$
(15)
an energy typically large enough to stabilize the open state. In section 2, we have made some rather general arguments about the shape of the bilayer deformation energy landscape. We now return to this picture briefly to discuss the consequences of spontaneous curvature. In our discussion, we assumed that the line tension, $`f`$, was typically positive, but we remarked that this need not be the case in the presence of large spontaneous curvature. If $`f`$ is negative, as depicted by the dotted curve in figure 2, the only state stabilized by the bilayer is the open state, which very naturally gives rise to the open state stabilization observed by Perozo et al. . Alternatively, this result can be understood from the predicted opening tension in equation 14. When $`\mathrm{\Delta }G_0`$ is bilayer deformation dominated, a negative $`f`$ implies that the opening tension is itself negative! A compressive force is required to stabilize the closed state. This argument gives a tantalizingly simple explanation for the open state stabilization but in the absence of measured values for the spontaneous curvature induced by LPC, we can only conclude that spontaneous curvature could stabilize the open states for rather generic values of the parameters. An experimental consistency check of these results is fairly simple. Perozo et al. have incorporated LPC asymmetrically. The same experiment might be repeated with $`H_{II}`$ phase inducing lipids which can also be used to generate spontaneous curvature but of the opposite sign. For the DOPC/DOPE system of Keller et al. , the spontaneous curvature is known and tunable as a function of concentration. Our results predict that $`\mathrm{\Delta }G_0`$ should be linear in $`C`$  and compatible with equation 15.

### 3.3 Transition Rates

Although we have not discussed transition states in detail, the proposed model also predicts that the transition rates between states depends on the line tension, $`f`$. If we assume the transition rates to be given by the standard Arrhenius form:

$$k=k_0\mathrm{exp}(G_B/kT),$$
(16)
where $`G_B`$ is the free energy barrier height for the transition and $`k_0`$ is roughly constant, the results can be understood intuitively by consulting figure 4. At fixed tension, forward transition rates (transitions to increasing radii) are reduced by line tension while backward rates are increased. At any fixed tension, we would therefore expect the forward rates for MscL reconstituted into a PC16 bilayer to be faster than for a PC18 bilayer while the backward rates are slower since the PC18 bilayer has a greater line tension. The ratio of these rates can be calculated but depends on the radial location of the transition state. At the opening tension for each bilayer, it can be shown that if the bilayer deformation energy dominates, the predicted rates for the transitions $`CO`$ will be faster for PC16 bilayer than for PC18 bilayer. Again, the ratio of these rates can be calculated in terms of the radii of the states and transition states. These two sets of qualitative rate predictions are straightforward to check experimentally.

In the argument above, we focused on the radial dependence of bilayer deformation free energy, and fixed the other components of the state vector, $`X`$. In principle, there is a potentially important piece of radial dependence we are missing. The internal conformation may effectively couple the radius to the other parameters in the state vector $`X`$, adding additional implicit radial dependence. For example, the thickness of the inclusion, $`W`$, is almost certainly a function of radius. It is also very natural to couple the mid-plane slope to the radius. We have ignored these dependences in order to develop an intuitive and simple one dimensional picture with as few undetermined constants and couplings as possible. Provided that the bilayer deformation energy change is dominated by the radial change, this simplified model is a useful tool for understanding the bilayer-inclusion interaction. More elaborate models might easily be built from the general analytic framework we have constructed. This framework will be described in a forthcoming paper.

## 4 Discussion

We have argued that the bilayer deformation energy is harnessed by MscL to govern channel gating. Indeed, we have shown that a model which attributes the entire free energy difference between the open and closed states to the bilayer deformation energy is compatible with the experimental data. These results are somewhat surprising since it has been shown experimentally that the mutation of a single residue in the vicinity of the channel gate, can significantly effect channel gating , . The protein conformational energy cannot be neglected in general. In fact, we have assumed that the protein conformational energy is large enough to constrain the channel geometry since we have assumed it is the bilayer which deforms rather than the protein. In principle, the closed state could have been stabilized by protein conformation alone, rather than mismatch, but exploiting bilayer deformation provides a robust mechanism for mechanotransduction, a design principle which functions in spite of the enormous number of nearly degenerate microstates endemic to proteins. Even for proteins as simple as myoglobin, Frauenfelder et al.  have shown that the macroscopic conformation corresponds to an enormous number of structurally distinct microstates. These ideas have already been exploited for channel proteins. Goychuk and Hänggi  have used this degenerate landscape to derive the empirical rate law for voltage gated channels. In light of these results, it is very natural to suppose that the protein conformational energy of the MscL protein gives rise to a vast number of nearly degenerate states as well. The bilayer deformation energy naturally breaks this degeneracy and forms a mechanotransducing channel. The ensemble of microstates we observe as the closed states is stabilized at low applied tension by the line tension, while the ensemble of microstates we observe as the open state is stabilized by high applied tension. The importance of bilayer deformation in mechanotransduction may help to explain why there are no obvious sequence motifs associated with mechanosensitivity  since a mismatch requirement does not imply sequence specificity. Harnessing bilayer elasticity does have one noted disadvantage. The gating of a channel will be affected by the membrane environment which surrounds it. This is precisely what the experiments of Perozo et al.  have shown. In realistic cell membranes there is an enormous diversity of proteins and lipids which would imply that the free energy, and therefore opening tension, in these membranes would be heterogeneous. Sukharev and co workers  have evidence for exactly this variability for MscL in giant spheroplasts.

The correspondence between our simple theoretical model for the gating of the MscL channel and experiment is at least strongly suggestive that this mechanism is exploited by the MscL channel. Our model can also naturally explain the stabilization of the open state by LPC  as well as the ephemeral nature of the substates . The tractable nature and simplicity of the model allow extensive analytic calculations to be made, which have in turn lead to numerous experimental predictions, discussed in section 3. Specifically we have predicted (i) a shift in the curve relating the free energy difference and acyl chain length when the hydrophobic thickness of the channel is altered, (ii) the dependence of the free energy on spontaneous curvature and, in particular, on the concentration of spontaneous curvature-inducing molecules, and (iii) the qualitative dependence of the transition rates on acyl chain length and spontaneous curvature.

We have developed an extensive framework for studying bilayer-inclusion interactions in the MscL system. The model we have discussed here is the simplest implementation of these results and a more thorough description of the model will appear in a future paper. There are several very natural extensions to the current work. For example, we have focused here on the radial dependence only, but, as we have briefly alluded to in section 3.3, there are two additional geometric parameters which may also play important roles in the function of the MscL channel. More detailed measurements of the rates and free energies of the various states will no doubt prove our simplified model incomplete and provide motivation and insight into a more detailed model of channel gating. The simplicity and generality of the competition between applied tension and line tension suggests that it may be a quite general phenomena for mechanotransduction. We hope to apply similar ideas to other mechanosensitive systems. More generally, we are also intrigued by the possibility of finding analogous bilayer deformation driven conformational changes for other transmembrane proteins which do not exhibit mechanosensitive function, perhaps illuminating a more general qualitative design principle for the function of transmembrane proteins.

## 5 Methods

### 5.1 Lipids

When discussing the lipids used by other authors, we have used the same naming convention they employed: 10:0 dicaproyl-phosphatidylcholine (PC10), 12:0 dilauroyl-phosphatidylcholine (PC12), 14:1 dimirstoyl-phosphatidylcholine (PC14), 16:1 dipalmitoleoyl-phosphatidylcholine (PC16), 18:1 dioleoyl-phosphatidylcholine (PC18,DOPC), 20:1 Eicossenoyl-phosphatidylcholine (PC20), lysophospholipid (LPL), lysophosphatidylcholine (LPC), dioleoyl-phosphatidylethanolamine (DOPE). We have preformed calculations based on elastic moduli measured by Rawicz et al. . To extrapolate the effects of changes in the acyl chain length, we use the moduli for a typical 18:1 phosopholipid (1,2-dioleoyl-sn-glycero-3-phosphocholine) and the scaling relations listed below

$$\begin{array}{ccccc}2a& =& 37\mathrm{\AA }\hfill & & a^1\\ K_A& =& 0.64kT\mathrm{\AA }^2\hfill & & a^1\\ K_B& =& 21kT\hfill & & a^3\\ 𝒦& =& 2.1\times 10^3kT\mathrm{\AA }^3\hfill & & a^0\end{array}$$
(17)
where $`2a`$ is the equilibrium thickness of the bilayer, $`K_A`$ is the area expansion modulus, $`K_B`$ is the bending modulus, and

$$𝒦\sqrt{2}\left(\frac{K_A^3K_B}{a^6}\right)^{1/4},$$
(18)
is the composite elastic modulus appearing in the line tension, $`f_U`$. The proportionality to bilayer thickness is approximate only and is deduced from treating the bilayer as a thin shell. For more elaborate models and scaling arguments, see Rawicz et al. . By fitting the Peak-to-Peak head group thickness  for saturated and mono-unsaturated lipids, we find that the relation between bilayer thickness and acyl chain length is roughly

$$2a=1.3n+17\mathrm{\AA }.$$
(19)
Again, more elaborate models and scaling arguments are found in Rawicz et al. .

### 5.2 Line tension

In terms of $`h(r)`$, the bilayer mid-plane position, and $`u(r)`$, half the difference of the bilayer thickness and the equilibrium thickness, the mean curvature contributions to the free energy density are

$$𝒢_\mathrm{B}=\frac{K_B}{2}[\underset{}{\underset{}{\left(^2h\right)^2+\left(^2u\right)^2}}\underset{}{\underset{}{C_+^2[h+u]C_{}^2[uh]}}],$$
(20)
where the variation of the $``$ terms contribute to the action in the bulk (bilayer), the $``$ terms are total derivatives and can be evaluated at the interface, the constant terms are dropped, $`K_B`$ is the bending modulus for the bilayer, and $`C_\pm `$ are the spontaneous curvatures of the upper and lower leaflets, respectively. The tension contribution to the free energy density is

$$𝒢_\alpha =\underset{}{\underset{}{\frac{\alpha }{2}\left(h\right)^2}},$$
(21)
where the thickness deformation term can be ignored. The interaction free energy density between the two layers is

$$𝒢_\mathrm{I}=\underset{}{\underset{}{\frac{K_A}{2a^2}u^2}},$$
(22)
where $`K_A`$ is the area expansion modulus and $`a`$ is the equilibrium thickness of a leaflet.

The equilibrium equations that result from the minimization of $`u(r)`$ and $`h(r)`$ are

$`0={\displaystyle \frac{\delta G[u,h]}{\delta u}}`$ $`=`$ $`\left[K_B^4+{\displaystyle \frac{K_A}{a^2}}\right]u(r),`$ (23)
$`0={\displaystyle \frac{\delta G[u,h]}{\delta h}}`$ $`=`$ $`\left[K_B^4\alpha ^2\right]h(r).`$ (24)
In the asymptotic regime, large $`R`$, we can ignore the curvature of the interface and replace the Laplacians with $`d^2/dr^2`$. For $`u(r)`$, the radially decreasing solutions are exponentials with complex wave numbers

$$\beta _\pm =\left(\frac{K_A}{a^2K_B}\right)^{1/4}e^{\pm i\pi /4}.$$
(25)
For $`h(r)`$, there is only a single radially decreasing solution with decay length:

$$\beta _h=\sqrt{\frac{\alpha }{K_B}}.$$
(26)
At the interface, the bilayer deforms to match the hydrophobic thickness of the protein:

$$u(R)=U=\frac{1}{2}Wa.$$
(27)
If the inclusion interface is relatively flat and the top and bottom surfaces of the bilayer are normal to the inclusion interface as depicted in figure 1, the radial slope of $`u(r)`$ is zero at the interface:

$$u^{}(R)=0,$$
(28)
although this too may be used as a geometrical parameter . For $`h(r)`$, we set the radial slope at the interface:

$$h^{}(R)=H^{}.$$
(29)
At infinity, the bilayer is unperturbed, implying

$`u(\mathrm{})`$ $`=`$ $`0,`$ (30)
$`u^{}(\mathrm{})`$ $`=`$ $`0,`$ (31)
$`h(\mathrm{})`$ $`=`$ $`0.`$ (32)
The equilibrium solutions can be found by matching the boundary conditions above.

The line tension is

$$f=_R^{\mathrm{}}𝑑r\left(𝒢_\mathrm{B}+𝒢_\alpha +𝒢_\mathrm{I}\right).$$
(33)
Integrating this equation by parts twice yields a bulk term proportional to the equilibrium equations, which is zero since the equations are satisfied, and surface terms. The total line tension is:

$$f=\underset{f_U}{\underset{}{\frac{1}{2}𝒦U^2}}+\underset{f_H}{\underset{}{\frac{1}{2}\sqrt{K_B\alpha }H^2}}+\underset{f_C}{\underset{}{K_BCH^{}}},$$
(34)
where $`𝒦`$ has been defined in equation 18, $`C\frac{1}{2}(C_+C_{})`$ is the composite spontaneous curvature, and $`f_H`$ is the mid-plane deformation contribution to the line tension which is typically small enough to ignore.

## 6 Acknowledgments

We are grateful to Doug Rees, Tom Powers, Jané Kondev, Klaus Schulten, Evan Evans, Sergei Sukharev, Eduardo Perozo, Olaf Andersen, Sylvio May, Ben Freund, and Mandar Inamdar for useful discussions, suggestions, and corrections. We also acknowledge support from the Keck Foundation, NSF and the NSF funded Center for Integrative Multiscale Modeling and Simulation. PAW acknowledges support through an NSF fellowship.


format: text

text:
# Temporal Bias in the Clustering of Massive Cosmological Objects

## I Introduction

Large-scale structure in the Universe is believed to have originated from a primordial Gaussian random field of matter fluctuations, the product of quantum fluctuations that were shifted to larger spatial scales during cosmological inflation \[e.g.  1\]. Cosmic Microwave Background observations show that when the Universe was 100,000 years old, the gaseous component was extremely smooth, with temperature variations $`10^5`$. These tiny inhomogeneities, in concert with inhomogeneities in the unseen massive dark-matter, were amplified through gravitational instability, eventually forming the galaxies, clusters, and other cosmological objects we see today.

While the precise details of structure formation are highly complex, the simplicity of the initial random field allows us to easily compute the overall distributions of cosmological objects. Galaxies and galaxy clusters represent the peaks in the initial density distribution which accrete matter at the expense of the diffuse regions between them, and thus their number densities can be simply related to the number densities of the peaks of the initial random field. This technique has been applied most cleanly to galaxy clusters, whose densities and evolution provide strong constraints on the overall matter density ek96 .

Similarly, because the peaks in a random field are more clustered than the overall distribution, the spatial clustering of cosmological objects is stronger than the underlying mass distribution. Furthermore, this “geometrical bias” is a systematic function of the mass of these structures, an effect that has been well-studied analytically and numerically \[e.g.  3, 4, 5\].

Yet, this is only half of the story. Here we conduct a detailed numerical simulation that shows that the spatial correlation function of objects that are in the midst of accreting substantial amounts of material is significantly enhanced over that of the general population. This temporal bias causes them to mimic the properties of higher-mass structures, with important astrophysical implications as discussed below.

The structure of this work is as follows: In §2 we describe our numerical simulation, discuss our group-finding algorithms, and develop a robust definition of accreting groups. In §3 we present our results for the spatial correlation functions of these samples, and in §4 we discuss the astrophysical implications of our results. Further details of this study are given in sc03 .

## II Simulations and Group Finding

Our numerical simulation traced the growth of primordial density fluctuations by dynamically evolving a large number of point test particles. The distribution of these particles was then used to determine the nonlinear evolution of the spatial correlation function of massive objects at late times. Driven by measurements of the Cosmic Microwave Background, the number abundance of galaxy clusters, and high redshift supernova estimates \[e.g.  7, 2, 8\] we focused our attention on a Cold Dark Matter cosmological model with parameters $`H=70`$ km/s/Mpc, $`\mathrm{\Omega }_0`$ = 0.3, $`\mathrm{\Omega }_\mathrm{\Lambda }`$ = 0.65, $`\mathrm{\Omega }_b=0.05`$, and $`\sigma _8=0.87`$, where $`H`$ is the Hubble constant today, $`\mathrm{\Omega }_0`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, and $`\mathrm{\Omega }_b`$ are the total matter, vacuum, and baryonic densities in units of the critical density, and $`\sigma _8^2`$ is the present variance of linear fluctuations on the $`8\times (100/70)`$ Mpc scale (where 1 Mpc is $`3.26\times 10^6`$ light years).

Periodic boundary conditions, which approximate large-scale homogeneity and isotropy, were taken, and the mass within our simulation volume was held fixed. Thus the overall box expanded along with the cosmological expansion, such that each side at any given redshift $`z`$ was $`73/(1+z)`$ Mpc across. This box was populated with $`350^3`$ dark matter particles that interact only gravitationally and represent the dominant mass component of the Universe. The mass of each particle, $`4.3\times 10^8`$ solar masses ($`M_{}`$), was chosen to match the observed mass density of the Universe, and the simulation was started at an initial redshift of $`z=49`$. The simulation used a parallel OpenMP-based version of the HYDRA code  with 64-bit precision.

To demonstrate the robustness of our results we have chosen two distinct group-finding approaches; the friends-of-friends approach da85  (FOF) and the HOP algorithm ei98 . FOF works by linking together all pairs of particles within a fixed “linking-length” of each other, and then taking each such group of “friends” to be an identified cosmological object. Although it remains popular, the FOF mass estimates are known to have significant scatter due to a problem that can occur as small strings of particles fall within the linking length.

The HOP algorithm works by using the local density for each particle to trace (‘hop’) along a path of increasing density to the nearest density maxima, at which point the particle is assigned to the group defined by that local density maximum. As this process assigns all particles to groups, a ‘regrouping’ stage is needed in which a merger criterion for groups above a threshold density $`\delta _{outer}`$ is applied. This criterion merges all groups for which the boundary density between them exceeds $`\delta _{saddle}`$, and all groups thus identified must have one particle that exceeds $`\delta _{peak}`$ to be accepted as a group (see ei98  for explicit details).

Beginning from $`z=4.89`$, we saved particle positions every 50 million years up to the final output at $`z=3`$. For the final 5 outputs we found FOF groups using a linking parameter of $`b=0.18`$, and HOP groups using the parameters: $`N_{dens}=48`$, $`N_{hop}=20`$, $`N_{merge}=5`$, $`\delta _{peak}=160`$, $`\delta _{saddle}=140`$, and $`\delta _{outer}=80`$. Visual inspection showed strong similarities between the two populations, with a small amount of unavoidable noise coming from groups around the 80 particle resolution limit (a group found by FOF at this limit may not be found by HOP and vice versa). We compare groups from one output to another by tracing back all particles with a given group index from the later output to the earlier output. Particles that show no membership to a group at the earlier time are regarded as ‘smooth infall’ while other non-null indices describe the merger history of the object.

To give a rough estimate of the accuracy of the group finding methods in Fig. 1 we plot the mass of the most massive progenitor at $`t_1(z=3.059)`$, versus the mass at $`t_2(z=3)`$, such that $`\mathrm{\Delta }t=5\times 10^7`$ years. As compared to the FOF groups, the smaller fraction of HOP groups lying above the equal mass line shows that the HOP algorithm identifies groups that are more likely to be massive at later outputs. The effect of this difference is significant.

Our definition of accreting groups is similar to that of pe03 , except that we select the subset that grew by 20% from $`t_1`$ to $`t_2`$, which implicitly includes mass accretion via smooth infall and results in 545(980) HOP(FOF) groups if $`\mathrm{\Delta }t=5\times 10^7`$ years. Note that the mass of each group is that at the end of each time interval, such that we tag all groups that experienced appreciable infall. The 20% value is arbitrary, but we selected it primarily because it appears to lie outside the central ‘noise’ band in the FOF data (see Fig 1). The groups corresponding to this cut appear in Fig 1, as points to the right of the dashed lines.

## III Temporal Bias

In Figure 2 we show the spatial correlation functions of the groups selected by both the HOP and FOF algorithms and compare them with $`\xi (r)`$ of the accreting groups. This function measures the excess probability of finding a pair of groups at a given separation $`r`$ relative to a random distribution, and is calculated for a separation bin $`r_l`$ as

$$1+\xi (r_l)=N(r_l)/N_{\mathrm{random}}(r_l),$$
(1)
where $`N(r_l)`$ is the number of pairs separated by distances between $`r_l`$ and $`r_{l+1}`$, and $`N_{\mathrm{random}}(r_l)=\frac{1}{2}N^2\frac{4\pi }{3}(r_{l+1}^3r_l^3)/V`$, with $`N`$ the total number of groups and $`V`$ the volume of the simulation. In the accreting case we co-added the correlation functions calculated from the differences from the last four $`\mathrm{\Delta }t=5\times 10^7`$ year intervals and the last two $`\mathrm{\Delta }t=10\times 10^7`$ year intervals. Radial bins of 1/80 the simulation size, corresponding to $`0.92`$ comoving Mpc, were taken throughout. For comparison, in each panel of Fig. 2 we also show the correlation function of all the groups in the next largest mass bin. The amplitudes of the correlation functions obtained using the full set of HOP and FOF groups agree with each other to within statistical uncertainties, as well as with analytical estimates.

The upper row demonstrates a clear enhancement of the clustering of accreting groups at both the $`10^{10.5}M_{}`$ and $`10^{11.0}M_{}`$ mass scales, with their correlation functions roughly matching those of objects three times greater in mass (no conclusion can be drawn from the high mass bin as the sample is too small). This “temporal biasing” arises from the fact that both objects accreting substructure as well as those experiencing considerable smooth infall tend to be found in the densest regions of space, which are themselves highly clustered. This conclusion is supported by the fact that the average local overdensity of groups in the $`10^{10.5}(10^{11})M_{}`$ mass bin is 0.82(0.87) (measured in 4 Mpc (comoving) spheres, corresponding to a mass scale of $`1.2\times 10^{13}`$ $`M_{}`$), whereas the same mass bin for the entire population exhibits an overdensity of 0.60(0.73).

In the second row of Fig. 2 we take a longer interval of $`\mathrm{\Delta }t=10\times 10^7`$ yr. This has a only slight dampening effect on temporal bias, which can not be definitively distinguished from statistical noise in our measurements. In the $`\mathrm{\Delta }t=20\times 10^7`$ yr case, however, only a very weak enhancement of $`\xi (r)`$ was measured.

In the lower two rows of this figure, we repeat our analyses using the FOF group finder. Although this approach is more susceptible to statistical noise, the same trends are apparent as in the HOP case. If $`\mathrm{\Delta }t=5\times 10^7`$ yr, this temporal bias is roughly equal to the geometrical bias of the groups three times more massive, while if $`\mathrm{\Delta }t=10\times 10^7`$ yr, $`\xi (r)`$ is boosted to a slightly lesser degree.

Finally, to quantify our results, we have computed the effective temporal bias in each mass bin, $`\mathrm{\Delta }t`$, and group finder. We define $`b_t^2`$ as the ratio of the correlation function of the accreting groups to the overall correlation function, weighted by the number of points in each bin in the overall function; $`b_t^2_{i=0}^{20}\frac{\xi _{\mathrm{accreting},i}N_{\mathrm{all},i}}{\xi _{\mathrm{all},i}N_{\mathrm{all},i}}`$, where the sum is carried out over all bins within $`r20`$ comoving Mpc. These values are labeled in each panel, and in the $`\mathrm{\Delta }t=20\times 10^7`$ yr case, $`b_t^2=`$ 1.1(1.0) in the $`10^{10.5}(10^{11.0})M_{}`$ HOP bins and 1.1(1.3) in the respective FOF bins.

## IV Astrophysical Implications

While temporal bias is a general property of the peaks of a gravitationally amplified Gaussian random field, our results have specific implications for the large sample of $`z3`$ galaxies made available by the Lyman-break color-selection technique st98 . Lyman break galaxies (LBGs) are observed to have enormous star-formation rates on the order of $`50M_{}`$ per year ad00 , implying that these objects are likely to be accreting large amounts of material. Furthermore, although the clustering of LBGs brighter than $`_{AB}25.5`$ is roughly that expected from the geometrical bias of $`10^{12}M_{}`$ objects \[e.g. 16, 17\], the linewidths measured from a spectroscopic sample these galaxies correspond to total masses $`10^{11}M_{}`$ pe01 .

To relate our result to LBGs we plot the spatial correlation function of $`_{AB}25.5`$ LBGs, as derived in we01 , in the center column of Fig. 2. Although there are significant uncertainties involved in computing this quantity, since comparisons are more naturally conducted in angular coordinates, the shaded regions provide a guide to the range of $`\xi (r)`$ values consistent with observations. In these panels, we see that if $`\mathrm{\Delta }t=5\times 10^7`$ yr is chosen, then temporal bias boosts the correlation function of $`10^{11}M_{}`$ groups into reasonable agreement with observations.

This mass is marginally consistent with the upper mass bound inferred from the rotation curves of a somewhat bright ($`_{AB}24`$) spectroscopic subset of LBGs pe01 . Furthermore, only $`4\%`$ of all groups exhibit appreciable accretion in each $`\mathrm{\Delta }t=5\times 10^7`$ year time interval and the density of $`10^{11}M_{}`$ groups is $`2\times 10^2`$ Mpc3, at $`z=3`$ in our assumed cosmology. Thus associating such objects with $`5\times 10^7`$ year starbursts results in a density $`5\times 10^4`$ Mpc3, comparable with that observed.

While quite suggestive, these comparisons are not meant as a complete model, and may not prove to be the final explanation of the discrepant mass estimates of LBGs. Kinematic models have been explored, for example, in which the observed velocity dispersions of LBGs are much less than the circular velocities of the groups in which they are contained mo99 . What is clear however, is that this bias can not be ignored and must be carefully considered when interpreting the clustering of these objects. While perhaps only part of the story, temporal biasing represents an important factor that must be taken into account when studying the properties of Lyman break galaxies.

###### Acknowledgements.

ES would like to express his sincere thanks for the hospitality shown to him by Jon Weisheit and the T-6 group at Los Alamos National Laboratory, where this work was initiated. We are grateful to Marc Davis for fruitful suggestions. ES was supported in part by an NSF MPS-DRF fellowship. RJT acknowledges funding from the Canadian Computational Cosmology Consortium and use of the CITA computing facilities. This work was supported by the National Science Foundation under grant PHY99-07949.


text:
# High-Temperature Atomic Superfluidity in Lattice Boson-Fermion Mixtures

(June 8, 2004)

## Abstract

We consider atomic Bose-Fermi mixtures in optical lattices and study the superfluidity of fermionic atoms due to $`s`$-wave pairing induced by boson-fermion interactions. We prove that the induced fermion-fermion coupling is always attractive if the boson-boson on site interaction is repulsive, and predict the existence of an enhanced BEC–BCS crossover as the strength of the lattice potential is varied. We show that for direct on-site fermion-fermion repulsion, the induced attraction can give rise to superfluidity via $`s`$-wave pairing, at striking variance with the case of pure systems of fermionic atoms with direct repulsive interactions.

The recent spectacular acceleration in the experimental manipulation Bloch1 ; Kasevich ; Bloch2 ; Modugno  and theoretical understanding Jaksch ; Zwerger  of systems of neutral atoms in optical lattices is leading to new and far reaching possibilities in the study of complex systems of condensed matter physics. Besides opening the way to the controlled simulation and experimental testing of models of strongly correlated systems Jaksch ; Hofstetter , such as, e.g., high-$`T_c`$ superconductors and Hall systems, atomic physics in optical lattices hints at the possibility of discovering and probing new quantum phases of matter. Among others, very interesting recent studies in this direction are concerned with lattice boson-fermion mixtures Albus ; Buechler ; Lewenstein ; Burnett ; Cramer .

Optical lattices are stable periodic arrays of microscopic potentials created by the interference patterns of intersecting laser beams. Atoms can be confined to different lattice sites, and by varying the amplitude of the periodic potential, it is possible to tune the interatomic interactions at will. Therefore, optical lattices provide an ideal tool to reach strong coupling regimes even in the dilute limit Zwerger . The theory of neutral bosonic atoms in optical lattices has been developed Jaksch  by assuming that the atoms are confined to the lowest Wannier band of the periodic potential. It can then be shown Jaksch  that the system is effectively described by a single–band Bose–Hubbard model Hamiltonian Fisher . In such a model the superfluid–insulator transition is predicted to occur when the on-site boson–boson interaction energy becomes comparable to some multiple of the hopping energy between adjacent lattice sites. This situation can be experimentally achieved by increasing the strength of the lattice potential, which results in a strong suppression of the kinetic (hopping) energy term. In this way, the superfluid–Mott-insulator quantum phase transition has been realized in a series of beautiful experiments by loading an ultracold atomic Bose–Einstein condensate in an optical lattice Bloch2 . The theory of dilute mixtures of interacting bosonic and fermionic neutral atoms subject to an optical lattice has been recently developed by deriving an effective single–band Bose–Fermi Hubbard (BFH) Hamiltonian from the underlying microscopic many-body dynamics Albus . The zero-temperature phase diagram of the BFH model has been analyzed both in the homogeneous and inhomogeneous cases Albus ; Buechler ; Lewenstein ; Burnett ; Cramer . These recent studies have shown that Bose-Fermi mixtures are very fundamental systems of condensed matter exhibiting complex phase diagrams and a potentially very rich physics that is just beginning to be unravelled.

The quest for superfluidity of fermionic atoms in dilute, degenerate quantum gases is the current object of widespread theoretical and experimental efforts Hofstetter ; Stringari ; Thomas . In the present work we propose a new mechanism for fermion pairing in atomic systems. We consider mixtures of bosonic and fermionic atoms in optical lattices. The fermions are assumed to belong to the same atomic species, but trapped in different spin states, so that they can differ by angular momentum or generalized spin $`\sigma `$ ($`\sigma =\{,\}`$). This situation can be easily achieved, for instance, by all–optical trapping. The mixture is then characterized by independent bosonic and fermionic hoppings and by direct on-site boson-boson, boson-fermion, and fermion-fermion couplings. For such a system we determine the structure of the induced fermion-fermion coupling due to the presence of the bosons, in analogy with recent and important studies on the effects of induced interactions in boson-fermion mixtures in free space Viveritall . We show that the regime of static linear response is an excellent approximation when the bosons are fully superfluid (weak optical lattice potential for the bosons) and, provided that the bosons are endowed with repulsive self-interactions, we find that the induced on-site fermion-fermion interaction is always attractive, thus leading to local pairing, irrespective of the sign of the boson-fermion interaction. In the case of a direct, repulsive on-site fermion-fermion interaction we show that the addition of the induced, attractive fermion-fermion interaction leads to an enhanced crossover regime between extended BCS pairing and Bose-Einstein condensation of local molecules, and allows for the occurrence of $`s`$-wave superfluidity of fermionic atoms. This is at striking variance with the case of pure Fermi systems with repulsive interactions, for which $`s`$-wave pairing is impossible and superfluidity might be achieved, at much lower temperatures, via higher-order mechanisms such as anisotropic $`d`$-wave pairing.

The single-band BFH system in an optical lattice with nearest-neighbor hopping and on-site direct boson-boson, boson-fermion, and fermion-fermion interaction is described by the Hamiltonian Albus

$`\widehat{H}`$ $`=`$ $`J_B{\displaystyle \underset{\{i,j\}}{}}\left(\widehat{a}_i^{}\widehat{a}_j+\text{h.c.}\right)J_F{\displaystyle \underset{\{i,j\},\sigma }{}}\left(\widehat{c}_{i,\sigma }^{}\widehat{c}_{j,\sigma }+\text{h.c.}\right)`$ (1)
$`+`$ $`{\displaystyle \frac{U_{BB}}{2}}{\displaystyle \underset{i}{}}\widehat{n}_i^B(\widehat{n}_i^B1)+U_{BF}{\displaystyle \underset{i}{}}\widehat{n}_i^B\widehat{n}_i^F`$
$`+`$ $`U_{FF}^{dir}{\displaystyle \underset{i}{}}\widehat{n}_{i,}\widehat{n}_{i,},`$
where $`J_B`$ and $`J_F`$ are, respectively the boson and fermion hopping amplitudes between adjacent lattice sites, and $`U_{BB}`$, $`U_{BF}`$, and $`U_{FF}^{dir}`$ are, respectively the direct boson-boson, boson-fermion, and fermion-fermion on-site interaction strengths. The operators $`\widehat{a}_i`$ are the on-site bosonic annihilation operators, and $`\widehat{c}_{i,\sigma }`$ are the on-site fermionic annihilation operators for states with spin $`\sigma `$. The bosonic on-site occupation number is $`\widehat{n}_i^B=\widehat{a}_i^{}\widehat{a}_i`$, while $`\widehat{n}_{i,\sigma }=\widehat{c}_{i,\sigma }^{}\widehat{c}_{i,\sigma }`$ denotes the fermionic on-site occupation number for states with spin $`\sigma `$, and, finally, $`\widehat{n}_i^F=\widehat{n}_{i,}+\widehat{n}_{i,}`$ denotes the total fermionic on-site occupation number. The bosonic and fermionic lattice potentials $`V_B`$ and $`V_F`$ are produced by a standing laser wave: $`V_{B,F}(x,y,z)=V_{B,F}(\mathrm{sin}^2kx+\mathrm{sin}^2ky+\mathrm{sin}^2kz)`$, where $`k=\pi /a`$ is the wave vector of the laser field, $`a`$ denotes the lattice spacing, and the lattice stregths $`V_{B,F}`$ that are felt by the bosons and the fermions depend on the intensity of the laser light and the detuning of the laser fequency against the respective atomic frequencies. Scaling the tunable lattice amplitudes $`V_{B,F}`$ with the bosonic and fermionic recoil energies, respectively $`E_B^R=(\pi \mathrm{})^2/2a^2m_B`$ and $`E_F^R=(\pi \mathrm{})^2/2a^2m_F`$, one can define the dimensionless bosonic and fermionic lattice amplitudes $`\eta _{B,F}=V_{B,F}/E_{B,F}^R`$. Given the wavelength $`\lambda _L`$ of the laser field, and the internal-state wavelengths $`\lambda _B`$, $`\lambda _F`$ for bosons and fermions, one has $`\eta _F/\eta _B=((\lambda _F/\lambda _B)^4(\mathrm{\Gamma }_F\mathrm{\Delta }\lambda _BE_B^R))/(\mathrm{\Gamma }_B\mathrm{\Delta }\lambda _FE_F^R)`$, where $`\mathrm{\Delta }\lambda _{B,F}=\lambda _{B,F}\lambda _L`$ are the boson and fermion detunings from the laser field, and the ratio of the natural line widths $`\mathrm{\Gamma }_F/\mathrm{\Gamma }_B`$ is usually of order one. By properly adjusting the detunings (blue or red) with the gauging of the laser field wavelength $`\lambda _L`$, one can realize lattices whose depths can either be very different or nearly equal for the bosons and the fermions. The hopping and on site energy amplitudes exhibit a characteristic behavior as a function of the parameters of the optical lattice: $`J_{B,F}=E_{B,F}^R(2/\sqrt{\pi })\eta _{B,F}^{3/4}\mathrm{exp}(2\eta _{B,F}^{1/2})`$; $`U_{BB,FF}=E_{B,F}^Ra_{BB,FF}k\sqrt{8/\pi }\eta _{B,F}^{3/4}`$, where $`a_{BB}`$ and $`a_{FF}`$ are the direct boson-boson and fermion-fermion $`s`$-wave scattering lengths, and $`U_{BF}=E_F^Ra_{BF}k(4/\sqrt{\pi })(1+m_F/m_B)(\eta _F^{1/2}+\eta _B^{1/2})^{3/2}`$, where $`a_{BF}`$ is the direct boson-fermion $`s`$-wave scattering length. The ratio of the hopping amplitudes reads $`J_B/J_F=[(E_B^R\eta _B^{3/4})/(E_F^R\eta _F^{3/4})]\mathrm{exp}[2(\sqrt{\eta _B}\sqrt{\eta _F})]`$.

The bosonic and fermionic dispersion relations in a $`3`$-D optical lattice with nearest-neighbor hopping take the following form:

$`ϵ_{B,F}(\text{q})`$ $`=`$ $`4J_{B,F}\left[\mathrm{sin}^2{\displaystyle \frac{q_xa}{2}}+\mathrm{sin}^2{\displaystyle \frac{q_ya}{2}}+\mathrm{sin}^2{\displaystyle \frac{q_za}{2}}\right],`$ (2)
where $`J_B`$ and $`J_F`$ are, respectively, the bosonic and fermionic hopping amplitudes. The bosonic energy spectrum in the Bogoliubov approximation reads

$$\omega _B(\text{q})=\sqrt{ϵ_B(\text{q})\left[ϵ_B(\text{q})+2n_BU_{BB}\right]}.$$
(3)
In the low momentum limit for the bosons, we recover the standard quadratic single particle spectrum $`ϵ_B(\text{q})3J_Ba^2q^2`$, where $`q^2|\text{q}|^2=q_x^2+q_y^2+q_z^2`$. Correspondingly, the phononic Bogoliubov excitation spectrum in the long wavelength regime reads $`\omega _B(\text{q})\sqrt{3J_B(aq)^22n_BU_{BB}}=c_Bq`$, where $`c_B=\sqrt{6n_BJ_BU_{BB}}a`$ is the velocity of sound.

Let us next consider the momentum and frequency dependent bosonic dynamic response function $`\chi _B(\omega ,\text{q})`$:

$$\chi _B(\omega ,\text{q})=\frac{2n_Bϵ_B(q)}{\omega ^2\omega _B^2(q)}.$$
(4)
Regarding the effects induced on the fermions, it is possible to show that the bosonic dynamic response function can be replaced by its static limit approximation. Let us focus on the fermion-fermion interaction $`U_{FF}^{ind}(\omega ,\text{q})`$ induced by the bosons. If the frequency contribution is not neglected, the gap function $`\mathrm{\Delta }(\omega ,\text{q})`$ for the fermions depends on frequency averages over the off-diagonal fermionic propagator $`G_{12}(\omega ,\text{q})`$.

$$\mathrm{\Delta }(\omega ,\text{q})=U_{FF}^{ind}(\omega ^{}\omega ,\text{q}^{}\text{q})G_{12}(\omega ^{},\text{q}^{})𝑑\omega 𝑑\text{q}^{}.$$
(5)
The maximum scale of integration for the fermions is fixed by the Fermi energy: $`\mathrm{}\omega _{max}E_F`$. The condition $`\omega _{max}^2\omega _B^2(q)|_{q=q_F}`$ can be satisfied for widely different scenarios of fermion filling factors. We will consider half-band fermion filling $`n_F1`$. This amounts to $`q_Fa\pi /2`$ and $`E_F6J_F`$. We then have $`\mathrm{}^2\omega _{max}^236J_F^2`$, and $`\omega _B^2(q)|_{q=q_F}36J_B^2+12J_Bn_BU_{BB}`$, so that the static limit condition reads

$$3J_F^23J_B^2+J_Bn_BU_{BB}.$$
(6)
This condition can be satisfied in several ways, depending on the type of boson-fermion mixture that one actually considers, and on the lattice strenghts and hopping amplitude ratios that can be realized in experiments. For boson-fermion mixtures with almost equal atomic masses, like 6Li-7Li, and large red detuning of the laser field from the typical $`D`$ lines of the Lithium isotopes, one has $`J_FJ_B`$, which implies a static limit condition of the form $`6(J_FJ_B)n_BU_{BB}`$. This condition is well satisfied at the onset of bosonic superfluidity $`U_{BB}10J_B`$ and bosonic filling factors $`n_B1`$. Alternatively, the static limit condition can be satisfied by realizing small detunings such that $`\eta _F>2\eta _B`$; this yields in general $`J_BJ_F`$, which again allows inequality (6) to hold. In this second instance, the bosonic lattice is usually some times weaker than the fermion lattice, and the bosons are in the deep superfluid regime. This setting is suited to treat mixtures irrespective of the atomic mass differences, and is thus the one we will consider in the following when studying the degenerate mixture of 40K and 87Rb. We will study quantitatively the case of half-band fermion filling $`n_F1`$, because this is a situation that can be easily realized experimentally, it provides the optimal scenario for the highest attainable critical temperatures, and finally allows a direct comparison with the work of Hofstetter et al. on the high-temperature superfluidity of pure fermionic systems in optical lattices Hofstetter . Taking the static bosonic response function $`\chi _B(0,\text{q})=2n_Bϵ_B(q)/\omega _B^2(q)`$, the resulting induced fermion-fermion interaction (due to the bosonic density-density fluctuations caused by the boson-fermion coupling $`U_{BF}`$) reads $`U_{FF}^{ind}(0,\text{q})=\chi _B(0,\text{q})U_{BF}^2`$. To discuss $`s`$-wave fermion pairing and evaluate the critical temperature, one must compute the average induced interaction $`\overline{U}_{FF}^{ind}`$ over the Fermi surface Viveritall . Provided that the bosonic lattice is not extremely weak (i.e. $`\eta _B`$ not $`1`$), contributions from nearest-neighbor induced interactions are negligible in first approximation, and one finds $`\overline{U}_{FF}^{ind}=U_{BF}^2/U_{BB}`$. We see that, irrespective of the sign of the boson-fermion $`s`$-wave scattering length $`a_{BF}`$, and thus of the boson-fermion on-site coupling $`U_{BF}`$, the induced fermion-fermion interaction is always attractive as long as the bosons are endowed with a positive $`s`$-wave scattering length $`a_{BB}>0`$ and thus with a repulsive on-site interaction $`U_{BB}>0`$. This is the case we will always consider in the following. For distinguishable fermions, such as ensembles of spin unpolarized neutral atoms, the total on-site fermion-fermion coupling $`U_{FF}^{tot}`$ is then

$$U_{FF}^{tot}U_{FF}^{dir}+\overline{U}_{FF}^{ind}=U_{FF}^{dir}\frac{U_{BF}^2}{U_{BB}}.$$
(7)
Thus, regarding interatomic pairing, fermions in a lattice Bose-Fermi mixture behave as interacting particles in a pure lattice fermionic system, but for the crucial difference of acquiring a dressed attraction that modifies the direct fermion-fermion interaction through the boson-fermion coupling. This has relevant consequences on the possible pairing mechanisms in the system, in view of the relative sign between the direct and the induced part of the coupling strength. Because the sign of $`U_{FF}^{dir}`$ is determined solely by the sign of the fermion-fermion $`s`$-wave scattering length $`a_{FF}`$, we need to treat separately the two possible instances $`a_{FF}<0`$ of direct on-site attraction, and $`a_{FF}>0`$ of direct on-site repulsion.

I) Fermions with direct on-site attraction – In this case the fermions can undergo a transition to $`s`$-wave superfluidity at a higher critical temperature $`T_c`$ than the analogous pure Fermi case, as the direct and induced interactions are both attractive and add together (See Eq. (7)). All the considerations valid for pure fermionic systems Hofstetter ; Ranninger  hold as well in the case of a Bose-Fermi mixture, with the crucial difference that the thermodynamic properties will now depend on the total fermion-fermion interaction $`U_{FF}^{tot}`$ rather than the direct one $`U_{FF}^{dir}`$. Starting with a lattice of low or intermediate depth, the BCS picture for $`s`$-wave pairing holds, predicting a critical temperature for the transition with a scaling behavior $`k_BT_c6J_F\mathrm{exp}(7J_F/|U_{FF}^{tot}|)`$. As the depth of the lattice is increased, the on-site interaction becomes at first comparable and then finally dominates over the tunneling amplitude; then the fermionic atoms form localized on-site bosonic molecules that can undergo a Bose-Einstein condensation (BEC) into the superfluid state at a transition temperature $`k_BT_c6J_F^2/|U_{FF}^{tot}|`$. In this situation the reduced ability of the pairs to move around in the lattice leads to a net decrease of the critical temperature. The maximum of the critical temperature $`T_c^{max}`$ is then achieved at the region of crossover between the BCS and BEC regimes, $`|U_{FF}^{tot}|10J_F`$ Scalettar ; Ranninger . In Fig. 1 we show the behavior of the critical temperature for the superfluid (SF) transition of 6Li atoms as a function of the strength of the optical lattice, at $`n_F1`$. When trapped in the state $`|,()`$, the fermionic Lithium atoms are endowed with a very large and negative scattering length $`a_{FF}2500a_0`$. This is the most favorable case of direct fermion-fermion attraction known so far, and considered by Hofstetter et al. Hofstetter . Their result is compared in Fig. 1 with the case of a 6Li-7Li mixture, for which $`a_{BF}38a_0`$ Schreck . We see that the presence of the bosons further enhances the maximum critical temperature attainable for the SF transition of the fermionic atoms. For a lattice spacing $`a=10^4a_0`$ the maximum critical temperature $`T_c^{max}=0.11\mu K`$ for pure 6Li; for the mixture 6Li-7Li we have $`T_c^{max}=0.14\mu K`$.

II) Fermions with direct on-site repulsion – In the case of direct on-site repulsive fermion-fermion interactions, $`U_{FF}^{dir}>0`$, the picture departs more radically from the analogous pure Fermi system. The direct repulsion and the induced attraction that determine the total interaction strength are now in competition, and the net total sign depends on their relative importance. In fact, because $`U_{FF}^{dir}`$ and $`U_{BB}`$ are both positive, we can evaluate the total fermion-fermion interaction as $`U_{FF}^{tot}=(U_{BB}U_{FF}^{dir}U_{BF}^2)/U_{BB}`$. In several important instances of mixtures of alkali atoms of current experimental interest, with direct fermion-fermion repulsion, the measured value of the boson-fermion $`s`$-wave scattering length $`a_{BF}`$ turns out to be much larger than both the boson-boson and fermion-fermion $`s`$-wave scattering lengths $`a_{BB}`$ and $`a_{FF}`$, with the largest difference attained in degenerate mixtures of fermionic 40K and bosonic 87Rb. Namely, in this case one has that $`a_{FF}=104.8\pm 0.4a_0`$ ($`a_0`$ being the Bohr radius) in the singlet state, and $`a_{FF}=174\pm 7a_0`$ in the triplet state Jin , while $`a_{BB}=100.2\pm 1.8a_0`$, and $`a_{BF}=410\pm 80a_0`$ Modugno2 . A similar situation occurs for a mixture of fermionic 6Li with direct repulsion and bosonic 7Li, where one has $`a_{BB}a_{FF}5a_0`$, and $`a_{BF}38a_0`$. In these important instances the direct fermion-fermion repulsion is thus dominated by the induced fermion-fermion attraction, and high-temperature superfluidity of the fermionic atoms occurs by $`s`$-wave pairing. This is at striking variance with the case of systems of pure fermionic atoms with direct on-site repulsion: for such systems $`s`$-wave pairing is obviously always forbidden, and superfluidity can be achieved via different mechanisms of higher order, at lower critical temperatures, such as anisotropic $`d`$-wave pairing. Therefore, mechanism II) realizes the first atomic analogue of electron paring in crystal lattices. A further dramatic increase in the transition temperature $`T_c`$ can be obtained by combining the boson-induced interaction with the experimental manipulation of the $`s`$-wave boson-fermion scattering length via Feshbach resonances. In fact, the existence has been recently predicted Ferlaino  of a Feshbach resonance at an applied magnetic field $`B725`$ G in the degenerate mixture of 40K and 87Rb, leading to an enhancement of the boson-fermion $`s`$-wave scattering length to the very high value $`a_{BF}687a_0`$. In Fig. 2 we show the behavior of the critical temperature for the superfluid transition of 40K atoms induced by the presence of the 87Rb atoms. Assuming the isotopes trapped at the $`D`$-line values of current experiments ($`795`$nm for rubidium, $`767`$nm for potassium), and a blue-detuned laser at $`764`$nm, we have $`\eta _F4\eta _B`$, so that in the range of values of $`\eta _F7`$ around the maximum critical temperature, we have $`\eta _B2`$ (fully superfluid bosons) and $`J_B^220J_F^2`$ (static regime holds). We compare results for the case with (continuous curve) and without the presence of a Feshbach resonance (dotted curve). For a lattice spacing $`a=10^4a_0`$ the first case yields $`T_c^{max}=0.05\mu K`$, while without resonance we have $`T_c^{max}=0.04\mu K`$.

In conclusion, we have shown that for mixtures of bosonic and fermionic neutral atoms in optical lattices, the fermionic atoms acquire an induced on-site attraction mediated by the boson-fermion interaction, irrespective of the sign of the latter. This fact allows the fermions to undergo a high-temperature superfluid transition via total $`s`$-wave pairing in both instances of direct attractive and repulsive on-site fermion-fermion interactions. Recent progresses in the manipulation of mixtures of bosonic and fermionic atoms in optical lattices Modugno  indicate that these predictions may be subject to experimental verification in the near future. We warmly thank Dr. K. Bongs for providing us with up-to-date values on forthcoming experiments with degenerate mixtures of 40K and 87Rb. This work has been supported by the ESF under project BEC2000+, the INFM, the INFN, and the Italian Ministery for Scientific Research, under project PRIN-COFIN 2002. We dedicate this work to Noam Chomsky on his 75th birthday, for his tireless struggle in favor of free research, human dignity and justice throughout the world.


text:
# The self-induced charge current in the experiments detecting photon rest mass

## Abstract

This note considers two related topics regarding the self-induced charge current in the experiments detecting photon rest mass: (i) the effective rest mass, $`(0.31)\times 10^{54}`$ Kg, of photons arising from the ions in the secondary cosmic ray is compared to Luo et al.’s new upper limit on photon rest mass ($`1.2\times 10^{54}`$ Kg); (ii) there are two main classes of experiments to detect photon rest mass based respectively on the Ampère-Maxwell-Proca equation and the Lagrangian (or Hamiltonian) density of electromagnetically interacting system.

PACS number(s): 03.50.De

During the last three decades, photon rest mass problem captured special attention of many investigators who have reported several experimental upper limits on the photon mass $`m_\gamma `$ by using various methods. Since the radiation from pulsars and/or certain explosive astrophysical events at high redshift can be used to place the severe limits on both the fractional variation ($`\mathrm{}c/c`$) in the speed of light with frequency and the photon rest mass, in 1969 Feinberg analyzed the data obtained in the detection of the sharply defined radio and optical pulses from pulsars and indicated that the rest mass of a real photon is less than $`10^{47}`$ Kg. Goldhaber et al. discussed the method used to set a limit on the photon rest mass (e.g., sensitive techniques which test Coulomb’s law and its analog in magnetostatics) and showed that the laboratory tests of Coulomb’s law gave an upper limit on photon rest mass of $`2\times 10^{50}`$ Kg and the constant “external” magnetic field experiment at the Earth’s surface obtained a limit $`m_\gamma 4\times 10^{51}`$ Kg. The experiments of Pioneer-10 measurement of Jovian magnetic field and hydromagnetic wave in the Earth’s magnetosphere and solar wind derived an upper limit on $`m_\gamma `$ of about $`10^{52}`$ Kg. In 1992 Chernokov et al. suggested a new method to set an upper limit on the photon mass at low temperature ($`1.24`$ K) based on a null test of Ampère’s law and their experiment resulted in an upper limit about $`10^{48}`$ Kg. Fischbach et al. investigated a new geometric limit on the photon rest mass which was derived from an analysis of satellite measurements of the Earth’s magnetic field. The order of magnitude of the upper limit on $`m_\gamma `$ obtained by Fischbach et al. was $`10^{51}`$ Kg. In 1998, Lakes reported that their experimental approach based on a toroid Cavendish balance which was used to evaluate the product of photon rest mass squared and the ambient cosmic magnetic vector potential gave an upper limit on $`m_\gamma `$ of about $`4\times 10^{52}`$ Kg. More recently, Luo et al. obtained the most new upper limit on photon rest mass of $`1.2\times 10^{54}`$ Kg by means of rotating torsion balance experiment.

Although both Luo et al.’s experimental scheme and their obtained results are excellent and impressive, here we have a supplement to their results, since our evaluation shows that the photon effective rest mass due to the self-induced charge currents in the environmental dilute plasma (e.g., the muon ($`\mu `$) component and alpha-particles in the secondary cosmic ray flux) is just the same order of magnitude as Luo’s obtained upper limit on photon mass.

In addition, it is worthwhile to point out that there are two main classes of experiments to detect photon rest mass based on the Ampère-Maxwell-Proca equation and the Lagrangian (or Hamiltonian) density of electromagnetically interacting system, respectively. The torsion balance experiments considered here fall into the latter class (i.e., based on the electromagnetic Lagrangian and Hamiltonian density).

In what follows we will discuss these two problems.

We consider the potential effects of self-induced charge currents arising in the torsion balance experiments. The Lagrangian density of electrodynamics reads

$$=\frac{1}{4}F_{\mu \nu }F_{\mu \nu }\frac{1}{2}\mu _\gamma ^2A_\mu A_\mu +\mu _0J_\mu A_\mu ,$$
(1)
where $`\mu _\gamma ^2=\left(\frac{m_\gamma c}{\mathrm{}}\right)^2`$ with $`m_\gamma `$, $`\mathrm{}`$ and $`c`$ being the photon rest mass, Planck’s constant and speed of light in a free vacuum, respectively, and $`\mu _0`$ denotes the magnetic permeability in a vacuum. Here the summation over the repeated indices is implied. It follows from (1) that the canonical momentum density of electromagnetic field reads $`\pi _\mu =\frac{}{\dot{A}_\mu }`$ or $`\stackrel{}{\pi }(x)=𝐄(x)`$. In the ions plasma, due to the conservation law of canonical momentum density of interacting electromagnetic system, i.e., $`\frac{\mathrm{d}}{\mathrm{d}t}(m𝐯+e𝐀)=0`$, we have $`𝐯=\frac{e}{m}𝐀+𝐂`$ with $`𝐂`$ being a constant velocity. By using the formula $`=𝐄\dot{𝐀}`$ for the system Hamiltonian and the electric current density $`𝐉=Ne𝐯=\frac{Ne^2}{m}𝐀+Ne𝐂`$ with $`N`$ standing for the volume density of ions , one can arrive at (in SI)

$$=\frac{1}{2\mu _0}\left[\frac{𝐄^2}{c^2}+𝐁^2+\left(\mu _\gamma ^2+2\frac{\mu _0Ne^2}{m}\right)𝐀^2+\frac{1}{\mu _\gamma ^2c^2}\left(𝐄\frac{𝐉_0}{ϵ_0}\right)^2\right]+𝐉_0𝐀_0Ne𝐂𝐀,$$
(2)
where $`𝐉_0=\rho `$ and $`𝐀_0=\varphi `$ are respectively the electric charge density and electric scalar potential, and $`ϵ_0`$ represents the electric permittivity in a vacuum and $`c=\frac{1}{\sqrt{ϵ_0\mu _0}}`$. Note that here $`\frac{1}{\mu _\gamma ^2c^2}\left(𝐄\frac{𝐉_0}{ϵ_0}\right)^2`$ results from both the mass term $`\frac{\mu _\gamma ^2}{c^2}𝐀_0^2`$ and the Gauss’s law $`𝐄=\mu _\gamma ^2𝐀_0+\frac{𝐉_0}{ϵ_0}`$.

It follows from (2) that the total effective rest mass squared of electromagnetic fields in ions plasma is given as follows

$$\mu _{\mathrm{tot}}^2=\mu _\gamma ^2+2\frac{\mu _0Ne^2}{m}.$$
(3)
It is known that at the sea level, the current density of muon component in secondary cosmic rays is about $`1\times 10^2\mathrm{cm}^2\mathrm{s}^1`$. Assuming that the muon velocity approaches the speed of light, the volume density of muon can be derived and the result is $`N_\mu =0.3\times 10^6`$ $`\mathrm{m}^3`$. So, according to the formula (3), the photon effective rest mass due to the self-induced muon charge current takes the form

$$m_{\mathrm{eff}}=\frac{\mathrm{}}{c^2}\sqrt{\frac{2N_\mu e^2}{ϵ_0m_\mu }}$$
(4)
with $`m_\mu `$ denoting the muon mass, electromagnetic wave with wavelength $`\lambda 100`$ m at the sea level acquires an effective rest mass about $`0.3\times 10^{53}`$ Kg. Hence, likewise, in the rotating torsion balance experiment the ambient cosmic magnetic vector potentials (interstellar magnetic fields) with low or zero frequencies may also acquire such effective rest mass.

Note that the above-mentioned muon volume density ($`0.3\times 10^6`$ $`\mathrm{m}^3`$) is the datum obtained only at the sea level. In order to evaluate the muon volume density in the environment where Luo’s experiment was performed, readers may be referred to the following handbook data of muon current density in underground cosmic rays: the muon current densities are respectively $`10^4`$ $`\mathrm{cm}^2\mathrm{s}^1`$ and $`10^6`$ $`\mathrm{cm}^2\mathrm{s}^1`$ at the equivalent water depths of $`100`$ m and $`1000`$ m under the ground. Since in Luo’s experiment the total apparatus is located in a cave laboratory, on which the least thickness of the cover is more than $`40`$ m, it is reasonably believed that the muon density of secondary cosmic rays in the vacuum chamber of Luo’s experiment may be one or two orders of magnitude less than that at the sea level. This, therefore, means that by using the formula (4) the effective rest mass acquired by photons in Luo’s vacuum chamber may be about $`(0.31)\times 10^{54}`$ Kg which can be compared to Luo’s result ($`1.2\times 10^{54}`$ Kg).

It is believed that the air molecules in the low-pressure vacuum chamber of torsion balance experiment cannot be easily ionized by the alpha-particles, electrons and protons of cosmic rays, because, for example, it is readily verified that the mean free path of alpha-particle moving at about $`10^7\mathrm{m}\mathrm{s}^1`$ in the dilute air with the pressure being only $`10^2`$ Pa is too large (i.e., more than $`10^6`$ m) for the dilute air molecules to be ionized. So, the air medium in the low-pressure vacuum chamber has nothing to contribute to the effective rest mass of photons. However, as far as the alpha-particle, $`\pi `$ mesons, electrons etc. in secondary cosmic rays are concerned, these ions will also contribute about $`10^{54}`$ Kg to photon rest mass, since, for instance, the volume density of alpha-particles and other ions at the sea level in secondary cosmic rays is about one order of magnitude less than that of muons.

Thus the photon effective rest mass resulting from both muons and other ions in secondary cosmic rays can be compared to the newly obtained upper limit ($`1.2\times 10^{54}`$ Kg) by Luo et al. in their recent rotating torsion balance experiment. For this reason, Luo’s experimental upper limit on photon mass is said to be just a critical value. It is reasonably believed that Luo’s recent experimental result is of physical interest. However, the photon mass problem still deserves further experimental investigation so as to improve the present upper limit on photon mass. It is claimed that, in the near future, only the contribution of ions in environmental dilute plasma, i.e., the secondary cosmic rays is ruled out, can we deal better with the experimental schemes, which can place severe limit on the photon rest (intrinsic) mass and other related experimental results. In addition, it is apparently seen that not only the photon intrinsic rest mass but also the effective mass due to media dispersion will play important roles in both fundamental and applied physics. For example, taking into account the influence of the media dispersion and the consequent photon effective mass on the geometric phases of photons propagating inside a noncoplanarly curved fiber is physically interesting, which is under consideration and would be published elsewhere.

Additionally, the difference between the two categories of $`\mu _{\mathrm{eff}}^2`$ based on different physical mechanisms should be pointed out: the theoretical mechanism of torsion balance experiment lies in the Lagrangian or Hamiltonian (2) of interacting electromagnetic system, so, the effective rest mass squared of electromagnetic fields is $`\mu _{\mathrm{eff}}^2=2\frac{\mu _0Ne^2}{m}`$. However, for those experimental schemes of photon rest mass stemming from the Ampère-Maxwell-Proca equation $`\times 𝐁=\mu _0ϵ_0\frac{}{t}𝐄+\mu _0𝐉\mu _\gamma ^2𝐀`$ and the consequent

$$\times 𝐁=\mu _0ϵ_0\frac{}{t}𝐄[\mu _\gamma ^2+\frac{\mu _0Ne^2}{m}]𝐀+\mu _0Ne𝐂,$$
(5)
the total mass squared $`\mu _{\mathrm{tot}}^2`$ should be $`\mu _{\mathrm{tot}}^2=\mu _\gamma ^2+\frac{\mu _0Ne^2}{m}`$ rather than that in (3). This minor difference between these two $`\mu _{\mathrm{tot}}^2`$ results from the derivative procedure applied to the Lagrangian density by using the Euler-Lagrange equation. Thus the effective rest mass squared of electromagnetic fields in the experiments based on the Ampère-Maxwell-Proca equation (5) is $`\mu _{\mathrm{eff}}^2=\frac{\mu _0Ne^2}{m}`$ rather than $`\mu _{\mathrm{eff}}^2=2\frac{\mu _0Ne^2}{m}`$. The so-called experimental realizations stemming from the Ampère-Maxwell-Proca equation are as follows: pulsar test of a variation of the speed of light with frequency; geomagnetic limit on photon mass based on the analysis of satellite measurements of the Earth’s field; experimental test of Ampère’s law at low temperature, etc.. In all these experimental schemes, the effective mass squared is $`\mu _{\mathrm{eff}}^2=\frac{\mu _0Ne^2}{m}`$.

It is shown that in the rotating torsion balance experiments the extra torque $`2\frac{\mu _0Ne^2}{m}𝐚_d\times 𝐀`$ produced by the interaction of magnetic dipole vector potential moment $`𝐚_d`$ with the ambient cosmic magnetic vector potential will also acts upon the toroid in the torsion balance experiments whose apparatus is immersed in the dilute plasma such as the cosmic rays and radon gas. Here the photon effective rest mass squared is $`\mu _{\mathrm{eff}}^2=2\frac{\mu _0Ne^2}{m}`$ rather than $`\mu _{\mathrm{eff}}^2=\frac{\mu _0Ne^2}{m}`$, as stated above.

Acknowledgements This project was supported partially by the National Natural Science Foundation of China under the project No. $`90101024`$.