text:
# An analysis of two-layer models for circumstellar disks

## 1 Introduction

Models of disks around pre-main–sequence stars are currently used to make predictions of a number of observable quantities and to compare them to observations. The complexity of such models varies from so-called power-law disks, where all the relevant disk properties are described by power-law of radius and changed independently (see, for example, Beckwith et al. 1990), to models where the thermal and geometrical properties are computed self-consistently, under the assumption of hydrostatic equilibrium (D’Alessio et al. 1998; Malbet et al. 1991,2001; Bell et al. 1997, 1999; Dullemond et al. 2002). A large effort is currently being made in improving the treatment of radiation transfer, by developing 2D codes that deal efficiently with the very large optical depths that characterize such disks (e.g. Bjorkman & Wood 2001, Dullemond 2002). At the same time, it is clear that the need of simple models, that can be easily used in interpreting the observations, remains. This explain the success of the two-layers disk schematization proposed some years ago by Chiang & Goldreich (1997; henceforth CG97).

Because the CG97 model and its derivatives (e.g. Chiang et al. 2001; Dullemond, Dominik & Natta 2001; Lachaume et al. 2003) are often used in interpreting disk observations, it is important to estimate their reliability and accuracy by comparing them to more detailed and self-consistent model calculations. At present, self-consistent 2-D and 3-D models are still too cumbersome for such a comparison. However, there are a number of 1D vertical structure models based on detailed radiative transfer, which can be used as templates to assess the reliability of 2-layer models.

In this paper, we will make use of the model of Dullemond et al. 2002 (henceforth DZN02). In this model the vertical temperature profile was solved using a detailed plane-parallel 1-D vertical radiative transfer procedure, which treats properly the slanted penetration of the stellar radiation (often referred to as 1+1D approximation). Using the resulting temperature profile, the density structure is re-computed after each iteration step by integration of the equation of vertical hydrostatic equilibrium. In a follow-up paper (Dullemond & Natta 2003) scattering of the stellar light by the dust particles was included in the model. For the case of zero albedo this model is in fact identical to the DZN02 model, and we will use DN03 as our “standard 1+1D disk structure model” against which we will compare the two-layer models.

The goal of this paper is three-fold. First, and most importantly, we compare an improved version of the CG97 model to the the 1+1D vertical structure model of DN03. This initial comparison is done for dust opacities in which scattering is switched off. Then we investigate how big the effects of scattering are in comparison to the typical errors of the improved CG model. Finally we discuss the possibility of using a simple recipe to include scattering into the CG model.

## 2 Description of the models

In this section we give a short description of the improved CG97 model we adopt, and of the 1+1D vertical structure model. Both families of models neglect viscosity, and treat the radial dependence of the surface density as a free parameter. Consequently, they also do not include viscous heating which, depending on the value of the viscosity parameter $`\alpha `$, could be of importance in the inner parts of the disk.

### 2.1 Improved CG97 model (”CG+”)

The central element of the two-layer Chiang and Goldreich model is that, due to the very shallow incident angle $`\alpha 1`$, the impinging stellar radiation is absorbed entirely in the very tenuous upper layers of the disk, high above the disk’s photosphere. While this layer has an optical depth of unity for stellar radiation at that grazing angle (typically around $`\alpha =0.05`$), it is very optically thin in the vertical direction ($`\tau _{\mathrm{vertical}}\alpha `$). Energy conservation requires that this absorbed radiation is re-emitted as infrared radiation, half of which is emitted away from the disk, while the other half is emitted towards the disk midplane. This downwards emitted flux is absorbed by the disk interior and re-emitted once more in the infrared, though this time at longer wavelengths. This typically gives two components in the SED: a warm optically thin component with dust emission features in the near- to mid infrared and a cool thermal blackbody component from the disk interior at longer wavelengths. If the disk is optically thick to the radiation from the midplane, then both components are equally strong, as the cool component is a reprocessed version of the downward radiated surface emission. The temperature $`T_s`$ of the surface layer is fixed: it is the optically thin dust temperature. The emissivity of the surface layer is determined by the surface density of the surface layer, which is calculated from energy conservation: the total emission should equal the total absorbed flux. The emission from the interior is determined by the temperature of the disk interior and the disk’s surface density. In the simple case of a grey opacity and high optical depth this so-called midplane temperature is determined by equating the thermally emitted flux $`F=\sigma T_i^4`$ to the downwards directed surface flux. In general, $`T_i`$ is computed using mean opacities both for the heating and for the cooling radiation.

The original paper of Chiang & Goldreich (1997) describes the two-layer flaring disk model in stages. At first, a grey opacity constant surface density model is presented for which the flaring index $`d\mathrm{log}(H_s)/d\mathrm{log}(R)=9/7`$, where $`H_s`$ is the height of the surface. Then a set of improvements are listed in order to account for effects of low optical depth at various wavelength regimes. In a later paper (Chiang et al. 2001) the model was refined to include more realistic opacities and, more importantly, the self-consistent (and numerically stable) determination of the flaring angle, which is crucial for energy conservation. Dullemond, Dominik & Natta (2001, henceforth DDN01) subsequently added the emission and structure of the disk’s inner rim, and included effects of self-irradiation and self-shadowing. The complete set of model equations is listed in DDN01. In the rest of this paper this improved Chiang & Goldreich model will be referred to as the CG+ model. We will not address here issues related to the inner rim.

### 2.2 Vertical structure models

The reference model against which we compare the CG+ model is the 1+1D vertical structure model of DZN02 and its improved version DN03. The radiative transfer is done in a two-stage procedure: first the primary stellar photons are followed as they penetrate the surface layers of the disk and get absorbed. Then the re-emitted infrared radiation is solved using a 1-D vertical frequency and angle-dependent radiative transfer code applied at each radial annulus. In a third stage the vertical pressure balance is solved, under the assumption that the gas temperature equals the dust temperature and that the gas-to-dust ratio is constant. This three-stage procedure is then iterated until a converged solution is reached for the complete temperature and density structure of the disk. During each iteration, the flaring index has to be recomputed in order to compute the flaring angle (“grazing angle” in the language of Chiang & Goldreich) self-consistently. This is done by estimating the surface height of each annulus by finding the height $`H_s`$ above the midplane where the direct stellar radiation is extincted by $`\mathrm{exp}(1)`$. The double-logarithmic derivative of this $`H_s`$ with cylindric radius R is then the flaring index from which the flaring angle can be computed (see DZN02).

Once the iteration procedure has converged, we have obtained a solution for the dust temperature $`T_{\mathrm{dust}}(R,Z)`$ and the disk’s gas+dust density $`\rho (R,Z)`$.

DN03 included scattering in the following way. Stage 1 of the DZN02 model (the irradiation of the disk by primary stellar photons) was replaced with a 1-D Monte-Carlo code. Primary stellar photons are followed as they penetrate into the surface of the disk and scatter around until they either get absorbed, or escape to infinity. If a photon is absorbed, it leaves behind its energy, which is then used as a source term for the radiative transfer of the second stage. In the second stage, which takes care of the re-emitted infrared radiation from the disk, scattering is not included. This is a valid approximation if the dust grains are small enough (typically not larger than 0.5 $`\mu `$m).

## 3 Comparing models without scattering

As a first step we compare the CG+ model and the 1+1D vertical structure for a test case without scattering. In this way the physical assumptions of both models are the same, so that in principle we expect both models to produce very similar SEDs. The comparison is first done for a single annulus at various radii and optical depths, irradiated by a T Tauri star or by a Herbig Ae/Be star. We then compare the SEDs of models of entire disks, as well as intensity profiles at some selected wavelengths.

### 3.1 Single annulus

As a first example we consider an annulus of a disk at 1 AU from the central star, with width 0.01 AU (i.e. an annulus between 1.00 and 1.01 AU). We take the central star to be a Herbig Ae/Be star ($`M_{}=2M_{}`$, $`R_{}=2R_{}`$, $`T_{\mathrm{eff}}=10000`$K). To simplify matters, we assume the spectrum of the star to be a blackbody spectrum with $`T=T_{\mathrm{eff}}`$. The grazing angle (i.e. the angle at which the stellar radiation enters the disk’s atmosphere) is computed self-consistently in both CG+ and 1+1D models, but for the purpose of the comparison described in this section we assume it to be fixed to $`\alpha =0.05`$. We also assume, for these annulus tests, that half of the surface of the star as seen from the top of the annulus is covered by the very inner parts of the disk (see CG97). This effectively reduces the stellar flux by a factor 0.5. The absorption cross section is that of a silicate grain of 0.1 $`\mu `$m (Draine & Lee 1984), with zero albedo. Different choices of the opacity will not change the results of this paper in any significant way, unless explicitely noted. The surface density of the disk is such that the total vertical absorption optical depth (from $`z=\mathrm{}`$ to $`z=\mathrm{}`$) at 550 nm is $`\tau _{550nm}=10`$. The mean molecular weight is assumed to be $`\mu =2.3`$.

In Fig. 1 the SED of this annulus is shown for both the vertical structure model and the CG+ model. The difference between the models is less than 30% at wavelengths $`\lambda >8\mu `$m. At shorter wavelengths the CG+ model has a bump that is not present in the SED of the vertical structure model. The reason can be traced back to the surface temperature assumed in the CG+ model, which is that of grains in an optically thin medium. In reality, however, only a fraction of the grains in the surface layer see the unextincted stellar light. More than 36 % of the material in the surface layer lies at optical depths larger than unity (along the grazing incident path of the stellar radiation). It is hard to assign a single temperature to the material in the surface layer. But a reasonable average temperature could be the temperature at an optical depth of $`\tau =0.7`$, which is the location where exactly half of the stellar radiation has been absorbed. The dashed line in Fig. 1 shows the SED as computed by the CG+ model in which the surface temperature is evaluated at an extinction of $`\mathrm{exp}(0.7)=0.5`$.

In Fig. 2 the temperature structure of the same annulus is shown as a function of the vertical coordinate $`Z/R`$, where $`Z`$ is the height above the midplane, and $`R`$ the radius from the star (in this case $`R=1\mathrm{A}\mathrm{U}`$). High above the surface of the disk the temperature is constant, and virtually equal to the optically thin dust grain temperature at that distance from the star. The slight difference is due to the fact that in the 1-D model the dust high above the disk is heated not only by the direct stellar radiation but also by the thermal emission of the disk below it.

As one goes towards smaller $`Z`$ the temperature starts to decrease due to extinction of the direct stellar flux by the material of the disk. This is the location of the disk’s “hot surface layer” which produces the dust emission features. The small tick mark on the temperature curve denotes the surface height $`H_s`$ where the grazing optical depth is unity, i.e. where most of the direct stellar radiation is absorbed. Continuing downwards from the surface layers, the temperature slope levels off again because the disk’s own infrared radiation becomes competitive with the direct stellar radiation in the thermal balance equation of the dust grains. This is where the disk’s photosphere is located and where most of the thermal continuum emission is produced (the midplane emission, in the terminology of CG models).

Over-plotted over the smooth temperature curve is the two-layer CG+ model, in its two variants mentioned above. The sharp step in the temperature profile is where the surface layer of the CG+ model starts: $`z=H_s`$. The main difference in the two variants of the model is the temperature of this surface layer. For one variant, the surface layer temperature equals that of optically thin dust grains. For the corrected variant this temperature is lower, and is closer to the temperature of the vertical structure model at the location of the tick mark.

Note that $`H_s`$ is different in CG+ and 1+1D models. Most of difference is due to the fact that in the 1+1D models the density in the disk surface deviates from the Gaussian shape predicted by isothermal models (see, for example, Fig. 3 of DN03). An additional, smaller difference is due to the fact that CG+ models use the mean opacity averaged over the stellar flux to define $`H_s`$, while the 1+1D models use the wavelength-dependent opacity.

Table 1 lists the parameters of two other annuli. Test A2 is the T Tauri analog of A1, with lower stellar mass and temperature, and the annulus taken at a smaller radius (to obtain similar temperatures in the disk). In Fig. 3 the SED of this test case is shown. The results are qualitatively similar to the case of the Herbig Ae star (test A1). In this particular example, however, the $`\mathrm{exp}(0.7)`$ correction of the surface temperature calculation in the CG+ model does not make a major improvement, although also in this case the overall SED of the 1+1D model is slightly better reproduced when this factor is taken into account.

As shown above, in many cases the CG+ model is reasonably well in agreement with the 1+1D vertical structure model. Yet, there are cases in which the differences are appreciable. We show an example in Fig. 4 which plots the SED of the A3 annulus. This annulus has $`\tau _V=300`$, and is at 50 AU from our Herbig Ae star. It is clearly seen that in this case the difference between CG+ and 1+1D predictions is rather large, though still within a factor of 2. The explanation for the breakdown of the CG+ model in this case is likely the use of mean opacities to compute the midplane temperature. As noted by DZN02, the midplane of the disk can cool down by leaking radiation at long wavelengths (where the optical depth of the disk becomes smaller than unity), more than expected on the basis of a mean-opacity analysis. This is clearly seen in Fig. 5, which shows the temperature structure of the A3 annulus. The 1+1D temperature profile, which goes from values lower than predicted by CG+ in the midplane to larger ones at intermediate $`z`$ results in lower fluxes around 100 $`\mu `$m and higher fluxes in the mid-infrared, so that the two-bump shape of CG+ practically disappears in the 1+1D model.

### 3.2 Full disk models

In this section we compare the SED and the physical structure of entire disks computed with the 1+1D vertical structure model and CG+ models. We present three illustrative cases: a disk around a Herbig Ae star, a T Tauri star and a Brown Dwarf. Table 2 lists the parameters of the three disk models.

We assume that the disk is in hydrostatic equilibrium and compute the flaring angle is self-consistently. This is done by evaluating the flaring index

$$\xi \frac{d\mathrm{lg}(H_s/R)}{d\mathrm{lg}R},$$
(1)
where $`H_s`$ is the height of the $`\tau _{\mathrm{grazing}}=1`$ surface of the disk. From the flaring index, the flaring angle can be computed directly:

$$\alpha =0.4\frac{R_{}}{R}+\xi \frac{H_s}{R}.$$
(2)
For numerical stability reasons this flaring index is always evaluated 2 grid points inwards of the point where it is used (see Chiang et al. 2001). In these models we assume that the stellar flux is not reduced due to occultation of the central star by the inner disk, as was assumed in the single-annulus models.

Fig. 6 shows the SED of model F1. Over-plotted are the two proposed CG+ models (the one with optically thin surface temperature and the one with the correction factor $`\mathrm{exp}(0.7)`$). One can see that overall the models agree, but that the CG+ model clearly shows a double-bump structure between 20 $`\mu `$m and 100 $`\mu `$m which is not seen in the full vertical structure model.

In Fig. 7 the same is plotted, but now for a T Tauri star (model F2). The disk is taken to be less massive than the Herbig star disk. The CG+ model works really well for this case, with errors mostly within 10%. Finally, in Fig. 8 the SED of a Brown Dwarf disk model is show (model F3). This disk is taken to have a much smaller outer radius and mass than the T Tauri case. Nevertheless, the disk remains very optically thick, as in the other models. Again, errors are reasonably small, though slightly larger than for model F2.

In Fig. 9 the height of the disk of model F1 (Herbig Ae star) and F2 (the T Tauri case) is shown as a function of radius, for both the 1+1D model and the CG+ model. The pressure scale height is derived from the temperature of the disk at the equatorial plane. At very small radii $`H_p`$ is virtually the same in CG+ and 1+1D models. At larger radii, $`H_p`$ in the 1+1D model drops below that of the CG+ model, because the equatorial plane temperature of the 1+1D model tends to be lower, as discussed in Section 3.1. The difference, however, is always very small.

The surface scale height in $`H_s`$ tends to be larger in 1+1D models, reflecting in part the deviation of the very upper layers from a Gaussian profile, as discussed of Section 3.1. The degree by which CG+ models can reproduce the results of 1+1D ones depends on the dust opacity in combination with other model parameters, such as the disk mass. Here we just note that the difference is larger for the HAe star, and for larger radii. However, these differences remain small and we can conclude that, in all cases we have considered, the disk shape can be derived with good accuracy using the simple approximations of the CG+ models.

Fig. 10 shows the flaring index as function of radius for model F2. This quantity gives information on the amount of radiation emitted by the disk at any given radius, because for $`RR_{}`$ the emission is proportional to $`\xi `$. Although the run of $`\xi `$ with $`R`$ is similar in the two models, we note that $`\xi `$ is a smooth function of radius in the CG+ models, but not in the 1+1D one. In particular, the kink in the flaring index 0.3 AU is not due to numerical inaccuracies but rather to the behavior of the midplane optical depth , which in this disk model drops below unity at long wavelengths at about that radius. The consequent cooling of the midplane below optically thick temperatures, discussed in §3.1, is seen as a slight kink in $`H_s`$ in Fig. 9. This effect, which depends on the disk and stellar parameters and on the exact wavelength dependence of the dust opacity, can have important consequences on the stability of 1+1D models.

Fig. 11 shows the brightness profiles of the two model approaches at two wavelengths, 10 $`\mu `$m and 1.3mm, where current instrumentation has the spatial resolution required to resolve the disks. The differences between the models are very small, probably within the observational accuracy of most experiments. The 10 $`\mu `$m profile has a clear exponential cut-off at large radii due to the fact that the dust becomes too cool to emit strongly at mid-infrared wavelengths. The CG+ and 1+1D models produce virtually identical emission at this wavelength, as seen also in the SED (Fig. 7). At 1.3 mm the differences are slightly larger (up to 30%), which reflects the fact that sometimes the 1+1D models can be slightly cooler in the disk’s midplane than the CG+ models, as discussed above.

## 4 Models with scattering

### 4.1 Effect of scattering on the disk

The effects of scattering of the stellar radiation on the structure and SED of circumstellar disks was discussed in detail in DN03. The main effect is that the thermal infrared flux from the disk is reduced at all wavelengths, with a larger effect at shorter than at longer wavelengths.

This reduction takes place because the scattering reflects away from the disk part of the stellar radiation, which therefore cannot be reprocessed into infrared radiation. The reflected fraction of the incident flux can be evaluated analytically using the H-functions of Chandrasekhar, which were discussed in this context in DN03. For a moderate value of the albedo ($`\omega 0.5`$) and isotropic scattering the reflected fraction is about 25% (of order $`\omega /2`$), and the reduction factor of the disk’s thermal emission is about 0.75, when averaged over the infrared spectrum. For higher albedo the reduction is larger, while it decreases if the scattering is more forward-peaked.

The reason why the reduction is stronger at short wavelength than at long wavelengths can be traced back to the fact the surface layers of the disk are affected more strongly by scattering than the interior. A photon that gets lost in the disk’s interior may scatter many times and still get absorbed rather than escape. Therefore the deep interior layers of the disk will have a temperature only slightly below the temperature of a disk without scattering.

The SEDs of a T Tauri disk (model F2) with increasing albedo, computed as in DN03, are shown in Fig. 12. The effect of reddening of the SED as $`\omega `$ increases is clearly seen. This suggests that CG+ model-predicted SEDs, which tend in general to overestimate the amount of flux at shorter wavelengths with respect to 1+1D models, may be significantly in error when realistic dust models are considered.

### 4.2 A simple recipe for scattering in CG+ models: is it possible?

Although the qualitative effect of scattering can be easily understood “a priori”, it is more difficult to estimate quantitatively how the infrared spectrum of the disk is affected at each individual wavelength. Still, it is important to explore the possibility of including scattering in the CG+ models in some simple way, and to assess the reliability of any such model by comparing its predictions to those of the more accurate 1+1D models.

A first recipe is based on the following considerations. In CG+ models, the stellar radiation is intercepted by the disk surface, which will re-emit (part as scattered light and part as reprocessed radiation) 1/2 of it toward the observer and 1/2 toward the midplane. This entire second 1/2 will be absorbed by the midplane, whose temperature, therefore, will not depend on $`\omega `$. Of the 1/2 emitted toward the observer, a fraction $`\omega `$ is in the form of scattered light, and consequently only a fraction (1-$`\omega `$) is in the form of thermal reprocessed flux at infrared wavelength. Since the temperature of this surface layer dust depends only on its absorption cross section, one can compute the SED by combining the emission of the midplane of a CG+ model with $`\omega =0`$ with the surface emission reduced by a factor (1-$`\omega `$/2). This recipe neglects multiple scattering and any wavelength-dependent effect. However, for low albedo it reproduces well the overall reduction of the infrared excess, and it is based on very simple arguments.

A second possibility, suggested also by CG97, is to apply the classical result that the fraction of incident photons absorbed by a plane-parallel slab is reduced by a factor $`\eta =\sqrt{1\omega }`$. The midplane will then have a temperature (in the optically thick case) lower than T($`\omega =0`$) by a factor $`\eta ^{1/4}=(1\omega )^{1/8}`$. The surface grain temperature is again unchanged, but the infrared emission of the surface needs to be reduced by a factor $`\eta `$. For not too high albedo, also this recipe will reproduce the overall reduction of the infrared excess found by DN03.

We show a comparison of these two prescriptions with the correct results of DN03 in Fig. 13. The top panel shows the comparison for the Herbig Ae star disk (F1), the bottom panel for the T Tauri star disk (F2). Each curve shows the fractional difference of the SED (reprocessed radiation only) of CG+ models, as labeled, to that of the comparison model, which is the 1+1D one with $`\omega =0.5`$ (isotropic scattering). The solid line shows the results for CG+ models where scattering is ignored. One can see that these models overestimate the flux at all wavelengths, and that the discrepancy can be larger than a factor of two (for the Herbig Ae disk) between about 2 and 5 $`\mu `$m. At longer wavelengths the agreement improves, and becomes quite good below about 100 $`\mu `$m.

The dashed and dotted lines show the results obtained with the CG+ model using recipe 1 and recipe 2, respectively. Both prescriptions seem to work equally well, and improve the SED description at wavelengths longer than about 2–3 $`\mu `$m, reducing the discrepancies below $``$ 20%. Note that the larger difference in the Herbig Ae SED around 40 $`\mu `$m is not related to scattering (see §3.2). Although we do not find any strong reason to prefer recipe 1 over recipe 2, we note that recipe 2 reproduces more correctly the overall reduction of the infrared excess for large values of $`\omega `$ (compare the run of $`\eta `$ with $`\omega `$ with the results of Fig.1 in DN03). Finally, it is interesting to notice that recipe 1 tends to suppress the 10 $`\mu `$m feature with respect to the continuum, while this is not the case for recipe 2.

## 5 Summary and conclusion

In this paper we analyze the reliability of two-layer models for passive irradiated flaring dusty disks around pre-main-sequence stars by comparing their SEDs with those computed from more realistic 1+1D vertical structure models. The two-layer model used in this paper is called CG+, and is an improved version of the Chiang & Goldreich model (see DDN01). The 1+1D model involves detailed 1-D vertical radiative transfer which treats scattering of the primary stellar photons penetrating the surface layers of the disk (see DN03). All the models presented in this paper can be downloaded as ASCII tables from http://www.mpa-garching.mpg.de/ PUBLICATIONS/DATA/radtrans/cgcompare/.

We start with a comparison of the SEDs for different disk models where scattering is ignored. In general, CG+ models over-predict the emission at short wavelengths. This is caused by the assumption that the surface grains have the temperature of grains in an optically thin medium. In many cases one can obtain a better agreement by computing the surface temperature with at an attenuation exp(-0.7) of the stellar flux. At longer wavelengths (mid and far-infrared), the CG+ SEDs sometimes shows two maxima, due to the two components of the model: one at mid-infrared wavelength (surface layer) and one at far infrared wavelengths (interior) (see, for example, Fig. 6). The 1+1D model, on the other hand, always produces a smooth SED (with the exception of dust features), which looks qualitatively somewhat different. However, a closer look shows that the agreement between different models is often quite good (within about 20%), with some potential problems deriving from the use of mean opacities in CG+ models. In our set of models, the largest discrepancies occur for the Herbig Ae star, where they reach $``$40% around 40 $`\mu `$m. The shape and strength of the 10 micron emission feature is in general well predicted by the CG+ models, and so is the peak-over-continuum ratio of this feature and its intensity profile. We therefore argue that also other dust features will be reasonably well reproduced by CG+ models. At millimeter wavelengths, both models predict similar fluxes (within 20% at most) and intensity profiles.

A comparison of the disk physical structure is also interesting. Both the pressure scale height $`H_p(R)`$, which depends on the midplane temperature, and the surface height $`H_s(R)`$ predicted by the CG+ model agree reasonably well with the 1+1D models.

Note that we have used in this paper an improved version of two-layer code, which does not make use of the analythical expressions for flaring angle and temperatures derived in CG97, but derives the disk structure self-consistently. A comparison of the CG97 analythical model with the results of 1D calculations can be found in Kraus & Krügel (in preparation).

CG+ models do not include scattering of the stellar light from dust grains. When compared to the more accurate 1+1D models which include scattering, even for quite low values of the albedo or rather peaked-forward scattering, we find that the differences in the SED can be larger than a factor of two in the near infrared. However, as the effect of scattering itself decreases at longer wavelengths, so do the differences between models. We have implemented two different simple recipes to our CG+ code to include scattering without changing the basic idea on which two-layer models are built. Both recipes seem to improve somewhat the agreement of CG+ models with the results of the 1+1D calculations, although recipe 2 is probably a better approximation for very large values of $`\omega `$. (see DN03).

In summary, the results from this paper show that CG+ models compare well with the much more detailed 1+1D vertical structure models. These two-layer models can therefore be used with reasonable confidence to fit observed SEDs of actual objects, provided that the parameters are not stretched too far. Problems could occur, for instance, when the disks are not optically thick enough, so that the annulus-by-annulus approach, on which both the CG+ and the 1+1D models are based, breaks down. In such a case one must resort to a full 2-D or even 3-D treatment of the problem. A future study should find out under which circumstances such a 2-D/3-D approach is unavoidable, or, oppositely, when the annulus-by-annulus approach (and thereby the CG+ approach) is justified.


format: text

text:
# “Water-free” computer model for fluid bilayer membranes

## I Introduction

When amphiphilic molecules such as lipids are brought into contact with water they tend to arrange so as to shield their “oily” hydrocarbon tail from the aqueous environment while exposing their hydrophilic head to the water. One of the simplest structures formed in this way is of a bilayer membrane - a double sheet of surfactants separating two aqueous phases israelachvili . Bilayer membranes are common in biological systems alberts . Living cells are separated from their extra-cellular surroundings by plasma membranes that control the transport of material into and out of the cell sackmann ; bloom . Most biological membranes are found in the fluid phase where the lipids comprising the bilayer can diffuse freely in the membrane plane. Another characteristic feature of lipid bilayers is their high flexibility which allows for large thermally-excited undulations safran ; lipowsky\_sackmann . The fluidity and low rigidity of membranes are important for many of their biological properties, such as their ability to change their shape easily and the possibility of proteins to insert themselves into the membrane gennis .

The thickness of membranes is comparable to the size of the constituting surfactant molecules (typically on the nanometer scale), while their lateral extension can greatly exceed their thickness and reach up to several micrometers. Consequently, coarse-grain phenomenological models, such as Ginzburg-Landau free energy functionals gompper\_schick  or the effective surface Hamiltonian safran ; lipowsky\_sackmann ; gompper\_schick ; helfrich , have been used in order to study the physical properties of membranes, as well as of other interfaces (like surfactants monolayer in microemulsions or vapor-fluid interfaces). In those theories the bilayer membrane is treated as a smooth continuous surface, and its elastic energy is related to the membrane area and the local curvatures. These theories have been very successful in describing the shape and phase diagrams of bilayer membranes seifert\_lipowsky ; wortis .

Phenomenological models describe the mesoscopic physical behavior of interfaces and membranes, but do not allow one to approach these systems on the molecular level. Many theories have been developed in an attempt to understand how the mesoscopic behavior emerges from the microscopic entities and the interactions between them. These theories include lattice “Ising-like” models gompper\_schick , molecular theories of the hydrocarbon chain packing szleifer , theories including the effect of electrostatic interactions andelman , and density functional theories chacon . The most microscopic detailed approach is employed in some computer simulations where the amphiphiles and water, and the interactions between them are modeled explicitly in full detail atomistic . Since these simulations require an enormously large computing time, they are restricted to fairly small systems consisting of 50-200 amphiphiles, and can be utilized to investigate phenomena occurring on short time scales of a few nanoseconds. In order to study mesoscale phenomena it is therefore necessary to dispense with some of the microscopic details in the simulations and use simplified models simplified\_review . A number of such simplified computer models have been devised by several groups. In these models the structure of the surfactant molecules is represented in a “coarse-grained” manner where a number of atoms are grouped together into a single site. The first level of coarse-graining is obtained by replacing the water molecules and the $`\mathrm{CH}_2`$ groups of the hydrocarbon chain by unified atoms simplified1a ; simplified1b ; simplified1c . This can reduce the number of atoms per lipids to about 50. Much more simplified models, in which the amphiphiles consists of only 5-10 atoms, were also presented smit ; goetz ; shelley . In these latter models the electrostatic potentials are usually ignored and the potentials of the chemical bonds are greatly simplified. At this level of simplification it is obviously impossible to address specific lipids systems, but rather the more general properties of self-assembling systems.

The size of amphiphilic systems which can be simulated using simplified models is constantly growing by virtue of the availability of inexpensive and powerful commodity PC hardware and due to the development of new simulation techniques such as Dissipative Particle Dynamics (DPD) dpd . Simulations of model systems consisting of $`N1000`$ lipids have been recently reported in the literature dpd ; lopez ; muller1 . The major restriction on the size of the systems in these simulation stems from the large number of atoms included in the simulation cell which is typically an order of magnitude larger than the number of amphiphilic molecules. The low ratio between the number of lipids and the total number of atoms is due to two factors. The first one is the number of atoms comprising each lipid molecule which, as discussed above, can vary from 50 to 5 depending on the level of simplification employed in the simulations. The second factor is the number of water molecules in the simulation cell. In bilayers simulations the typical number ratio of water to lipid molecules is in the range of 10 to 30 atomistic ; simplified1a ; simplified1b ; simplified1c ; smit ; goetz ; shelley . A great fraction of the computing time is, thus, “wasted” on the simulations of the water even when the water molecules are represented by a single (unified) atom. Only very few models have so far been proposed in which the amphiphiles are simulated without the presence of water. The major difficulty in establishing such “water-free” models is the need to mimic the hydrophobic effect that prevents the amphiphilic molecules from leaving the aggregate into the solvent. Drouffe at al. drouffe  and Noguchi at al. noguchi  have used ad hoc multibody potentials to overcome this problem. With the aid of these nonphysical potentials they have managed to observe the formation of fluid vesicles in their simulations. La Penna at al. lapenna  have studied a water-free flat bilayer model with Lennard-Jones (LJ) potentials that depend on the relative orientation of the lipids (and which are closely related to the Gay-Berne gay\_berne  potentials used in liquid crystals simulations). With this model they have been able to simulate bilayer membranes in both the solid and the gel phases. Fluid membranes, however, were found unstable against lipids evaporation from the membrane plane.

In this paper we present an exceptionally simple computer model of a fluid bilayer membrane. Our model has the following features: (a) It is a water-free model, i.e., we simulate the membrane without the presence of water. (b) The “lipids” forming the membranes consist of only three atoms, one representing the hydrophilic headgroup and the other two the hydrophobic tail. These three atoms are “glued” to each other to form a rigid linear trimer (the lipid), and have no additional interactions between them. (c) The different lipids interact through finite range (truncated) LJ interactions between their three sites. The parameters of the LJ potentials are fixed and do not depend neither on the relative orientation of the lipids (as in Ref. lapenna ), nor on their local density (i.e., there are no multibody interactions in our model). The above mentioned properties make our membrane model computationally very efficient (albeit a less “flexible” one in comparison to other simplified models with more interaction sites per amphiphile). To investigate the statistical mechanical properties of the membrane, we have performed a set of Monte Carlo (MC) simulations where for each MC run we have fixed the temperature, the number of lipids, and the projected area of the membrane. The projected area serves as the control parameter in our simulations, and we have investigated the phase behavior of the membrane as a function of it. We found that upon increasing the projected area (i.e., reducing the area density of the lipids) the membrane undergoes a solid-fluid phase transition. In the solid phase the lipids are not mobile and they pack in a hexagonal order. In the fluid phase the lipids are free to diffuse in the membrane plane. We have measured the spectrum of thermal undulations of the fluid membranes from which we have extracted the surface tension and the bending modulus that characterize the elastic behavior of the membrane. At low area densities we found another transition from negative to positive surface tension, accompanied by the formation of pores in the membranes. Such a behavior is indeed predicted by theoretical arguments litster ; netz ; muller2

The paper is organized as follows: In section II we present our computer model, and discuss the details of the simulations. In section III we describe the physical properties of the systems as obtained by the simulations. The section is divided into three subsections dealing, respectively, with the phase diagram of the system, its spectrum of thermal undulations and elastic properties, and the appearance of holes in membranes with large projected area. We summarize and discuss the results in section IV

## II Details of the Model and the Simulations

The lipids in our model system consist of three spherical atoms connected to form a linear trimer. The lipid molecules are rigid - they do not bend and the distance $`\sigma `$ between the center of the atoms is fixed (see fig. 1). We set $`\sigma =1`$ as our unit length scale throughout this paper. We shall label the three atoms forming each lipid as 1, 2, and 3. Atom 1 represents the hydrophilic head of the lipid, while atoms 2 and 3 represent its hydrophobic tail. The different lipids interact with each other via spherically symmetric pair potentials between their constituting atoms. The pair potential $`U_{ij}(r)`$ depicts the interactions between atom $`i`$ and atom $`j`$ of two different molecules separated a distance $`r`$ apart. The pair potentials $`U_{12}`$ and $`U_{13}`$ describe the interaction between hydrophobic and hydrophilic particles. They are given by the purely repulsive LJ potential

$$U_{12}^{\mathrm{LJ}}(r)=4ϵ_{12}\left(\frac{\sigma _{12}}{r}\right)^{12},$$
(1)
and

$$U_{13}^{\mathrm{LJ}}(r)=4ϵ_{13}\left(\frac{\sigma _{13}}{r}\right)^{18}.$$
(2)
The pair potentials $`U_{11}`$, $`U_{22}`$, and $`U_{33}`$ describe the interactions between two similar atoms, both either hydrophobic or hydrophilic. They are given by the attractive LJ potentials

$$U_{ii}^{\mathrm{LJ}}(r)=4ϵ_{ii}\left[\left(\frac{\sigma _{ii}}{r}\right)^{12}\left(\frac{\sigma _{ii}}{r}\right)^6\right],$$
(3)
where $`i=1,2,3`$. Finally, the interaction between the hydrophobic particles 2 and 3 is also depicted by an attractive LJ potential, but of the form

$$U_{23}^{\mathrm{LJ}}(r)=4ϵ_{23}\left[\left(\frac{\sigma _{23}}{r}\right)^2\left(\frac{\sigma _{23}}{r}\right)\right].$$
(4)
All pair potentials are truncated at the same cut-off $`r_c=2.5\sigma `$, and the discontinuity at $`r_c`$ is avoided by adding extra terms to the LJ potentials that ensure the vanishing of the potential, as well as of its first and second derivative, at $`r=r_c`$. The final form of the pair potentials is thus given by

$`U_{ij}(r)=U_{ij}^{\mathrm{LJ}}(r)`$ $``$ $`U_{ij}^{\mathrm{LJ}}(r_c){\displaystyle \frac{U_{ij}^{\mathrm{LJ}}(r)}{r}}|_{r=r_c}(rr_c)`$
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2U_{ij}^{\mathrm{LJ}}(r)}{r^2}}|_{r=r_c}(rr_c)^2.`$ (5)
The different pair potentials are depicted in fig. 2. The values of the parameters $`\sigma _{ij}`$ (in units of $`\sigma `$) and $`ϵ_{ij}`$ (in $`kT`$ units, where $`T`$ is the temperature and $`k`$ the Boltzmann constant) used in the simulations are summarized in the caption on fig. 2.

The pair potentials in our computer model have been designed to allow, on the one hand, the diffusion of molecules in the plane of the membrane but to restrict, on the other hand, their motion in the third direction. We have tested various models before we arrived to the one that we have used in the simulations. The original idea was to use dimers with one hydrophilic and one hydrophobic particles, and to describe the interactions between them by 6-12 LJ potentials \[Eq.(3)\] and a 12-power repulsive potential \[Eq.(1)\], depending on whether the atoms are of the same or different species. It turned out that the membranes depicted by such a model were unstable against the extraction of molecules from the membrane plane. To increase the membrane stability we added a third hydrophobic atom to the lipids. The pair interactions between this atom and the other two atoms are described by different forms of LJ potentials: For the interaction with the hydrophilic atom labeled 1 we use the more repulsive 18-power LJ potential $`U_{13}`$ \[Eq.(2)\], while for the interaction with the hydrophobic atom labeled 2 we use the 1-2 LJ potential $`U_{23}`$ \[Eq.(4)\]. The former potential establishes a strong repulsion between the hydrophobic and the hydrophilic parts of the lipids, thus reducing significantly (eliminating on the time scale of the simulations) the escape probability of molecules (more on this point in the next paragraph). The latter has a very shallow minimum which allows a greater mobility of the lipids in the membrane plane (by making small the energy changes due to a relative motion of the lipids with respect to each other). We have gone through a rather lengthy “trial and error” process of fine tuning the parameters $`\sigma _{ij}`$ and $`ϵ_{ij}`$ which control the range of pair repulsion and the depth of the attractive potential wells. Their values have been set to (a) make favorable the alignment of molecules next to each other at a distance slightly larger than $`\sigma `$, and (b) to make the attraction between molecules sufficiently strong to support the stability of the membrane, but not too strong to the extent that would entirely prevent the diffusion of the lipids.

It is not an easy task to form a fluid bilayer sheet in a model system that does not contain water. Membranes become fluid at low area densities and high temperatures, and under these conditions the lipids tend to escape quite easily from the membrane plane. It is the water that confines the lipids to the membrane. In the absence of water molecules this role has to be played by the hydrophobic heads which must form some sort of geometric or dynamic constraint for the extraction of lipids. In our model we establish such a constraint by making the excluded volume part of the pair potentials $`U_{ij}`$ non-additive, namely we make the size of a particle $`i`$ “seen” by another particle $`i`$ smaller than its size as seen by a particle $`j`$ of a different species. We can define the distance $`a_{ij}`$ at which the pair potential between them $`U_{ij}=kT`$ as a measure for the range of hard core repulsion between the two particles $`i`$ and $`j`$. (It is unlikely to find a pair $`i`$ and $`j`$ separated by a smaller distance.) It is customary to regard $`a_{ii}`$ as the diameter of atom $`i`$ and, with this interpretation, to expect for the additivity of the hard core diameters, i.e., to have $`a_{ij}(a_{ii}+a_{jj})/2`$ for $`ij`$. In our model we do not find this property (see fig. 2). The pair potentials in our system describe the effective interactions between the different atoms and they include the effect of the water molecules which are not simulated explicitly. Therefore, there is no a priory reason why the effective diameters associated with different particles should be strictly additive. The increased range of hard core repulsion between the hydrophilic atom 1 and the hydrophobic atoms 2 and 3 is designed to compensate for the absence of water from the simulation cell.

The simulations were performed with membranes consisting of $`N=1000`$ lipids (500 lipids in each layer) with periodic boundary conditions in the membrane $`(x,y)`$ plane, and with no boundaries in the normal $`z`$-direction. Subsequent MC configurations were generated by two types of move attempts: translations of lipids and rotations around the mid (second) atom. The MC unit time is defined as the time (measured in number of MC configurations) in which, on the average, we attempt to move and rotate each molecule once. The acceptance probability of both types of moves was approximately half. We performed a set of simulations of membranes with the same temperature $`T`$ and number of lipids $`N`$, and with varying projected areas. For each value of the projected area we studied 8 different membranes starting at different initial configurations. The initial configurations were created by randomly placing 500 lipids in two layers with a vertical (along the normal $`z`$ direction) separation $`\sigma `$ between the atoms labeled 3 in the two layers, and with all the lipids oriented normal to the membrane plane. The initial configurations were then “thermalized” over a period of $`510^5`$ MC time units, followed by a longer period of $`610^6`$ time units during which quantities of interest were evaluated. The duration of the MC runs is substantially larger than the relaxation time which we estimated in various ways: As a first approximation for the relaxation time we used the time it took the potential energy of the membrane to saturate from its high initial value (resulting from overlap of particles in the random initial configuration) to a final “typical” value. This time was of the order of $`10^4`$ MC time units. An independent estimate of the relaxation time was obtained from a study of the spectrum of thermal undulations of the membranes (see more details, later in the text). Inspection of the autocorrelation function of the amplitude of the longest wavelength mode led to a similar estimate of $`10^4`$ time units for the relaxation time. A more conservative estimate can be obtained from measurements of the self-diffusion constant of the lipids in the fluid phase (see, again, later in the text). The relaxation time can be associated with the time it takes a molecule to diffuse a distance equal to the pair potentials cut-off ($`2.5\sigma `$). The relaxation time obtained using this criterion was an order of magnitude larger ($`10^5`$ MC time units), still smaller than the equilibration time, and much smaller than the total length of the simulations.

## III Simulation Results

### III.1 Phase Diagram

The projected area of the membranes in the simulations ranges from $`A_p=L_p^2=(26.875)^2`$ to $`A_p=(30.625)^2`$ with intervals of $`\mathrm{\Delta }L_p=0.625`$. For all area densities we measured the self-diffusion constant of the molecules relative to the diffusion of the center of mass, defined by

$$D\underset{t\mathrm{}}{lim}\frac{\mathrm{\Delta }r^{}(t)^2}{4t}\underset{t\mathrm{}}{lim}\frac{1}{4Nt}\underset{i=1}{\overset{N}{}}\left[\left(\stackrel{}{r}_i(t)\stackrel{}{r}_{\mathrm{CM}}(t)\right)\left(\stackrel{}{r}_i(0)\stackrel{}{r}_{\mathrm{CM}}(0)\right)\right]^2,$$
(6)
where $`\stackrel{}{r}_i(t)`$ denotes the position of the $`i`$-th lipid (defined by the position of its mid atom) at time $`t`$, while $`\stackrel{}{r}_{\mathrm{CM}}(t)`$ denotes the position of the center of mass of lipids remark0 . We have also measured the self-diffusion coefficient in the membrane plane, defined by:

$$D_{xy}\underset{t\mathrm{}}{lim}\frac{1}{4Nt}\underset{i=1}{\overset{N}{}}\left\{\left[\left(x_i(t)x_{\mathrm{CM}}(t)\right)\left(x_i(0)x_{\mathrm{CM}}(0)\right)\right]^2+\left[\left(y_i(t)y_{\mathrm{CM}}(t)\right)\left(y_i(0)y_{\mathrm{CM}}(0)\right)\right]^2\right\},$$
(7)
where $`x`$ and $`y`$ denote Cartesian coordinates. (In all the simulations the membranes lied in the $`(x,y)`$ plane, while fluctuating in the normal $`z`$-direction.) As the lipids can only diffuse within the plane of the membrane, we found no difference between $`D`$ and $`D_{xy}`$.

At low projected area (high area density) we found the membrane in a solid phase characterized by two features: (a) The diffusion constant of the lipids is vanishingly small. \[The root mean square displacement $`\sqrt{\mathrm{\Delta }r^2}`$ has barely changed during the course of the simulations, and it has never exceeded the typical distance between neighbor molecules ($`\sigma `$).\] (b) The lipids arrange themselves in a hexagonal order in the membrane plane. A typical equilibrium configuration of a membrane with $`A_p=(26.875)^2`$ is shown in fig. 3 (a). A top view of the plane of mid (labeled 2) atoms of the membrane upper layer, revealing the hexagonal order of the lipids, is shown in fig. 3 (b). The lattice imperfections observed at fig. 3 (b) should be mainly attributed to the incommensurability of the 500 sites hexagonal lattice with the square simulation cell.

At larger values of the projected area \[$`A_p(28.125)^2`$\] we found the membranes in a fluid phase. The main feature that distinguishes fluid from solid membranes is the diffusion of the lipids. In fig. 4 we plot lipids mean square displacement $`\mathrm{\Delta }r_{}^{}{}_{}{}^{2}`$ \[see definition in Eq.(6)\] as a function of the simulation time $`t`$ for fluid membranes with different area densities. The slope of the asymptotically linear curves is four times larger than the self-diffusion constant $`D`$. One can observe the growth of $`D`$ with the increase of the projected area - a rather expected observation as the increase of the projected area means more room for the molecules to move. A typical equilibrium configuration of a fluid membrane with $`A_p=(28.75)^2`$ is depicted in fig. 5 (a). Another characteristic feature of the fluid membranes is the loss of in-plane hexagonal order, as demonstrated in fig. 5 (b) \[compare with fig. 3 (b)\].

The membranes with $`A_p=(30.0)^2`$ exhibited an interesting feature - they developed pores, as demonstrated in the configuration shown in fig. 6. These pores tended to appear irregularly in the membrane with a characteristic time scale $`\tau 210^5`$ for the formation of a pore, and a typical pore life time of a few thousand time units. Another interesting phenomenon which we observed for this value of $`A_p`$ and did not observe at lower projected areas was the occurrence of “flip-flops” - the transition of lipids from one layer to the other. In fig. 7 we look at the same membrane depicted in fig. 6. In this figure, however, we plot only the 500 lipids that were located in the upper layer in the initial configuration. About 30 of them have managed to diffuse from the upper to the lower layer during the course of the simulations. A similar (although not necessarily identical) number of lipids have moved in the opposite direction. Trans-bilayer diffusion is an important process in real bilayer membranes sackmann2 . To allow for uniform bilayer growth, some of the lipids must be transfered from one leaflet to the other during the self-assembly process. When a flat bilayer is bent to form a spherical vesicle, the area of the inner layer becomes smaller than the area of the outer layer, and it is the transition of lipids from the former to the latter that balances their area densities. It has been suggested, based on experiments matsuzaki  and computer simulations okazaki , that that the formation of pores and the flip-flop motion are closely interconnected. According to these studies, the pores provide a transverse diffusion conduit for the lipids, through which their hydrophilic headgroups cross the hydrophobic region of the membrane. Our study supports the conjecture about the possible relation between pores and flip-flops. In section III.3 we discuss the origin of the formation of the pores. We show that it is associated with the change in the sign of the membrane surface tension - from a negative to a positive value. Once the surface tension attains a sufficiently large positive value, the energy involved with the formation of the pores is compensated by the reduction in elastic energy.

While the formation of the pores allowed the diffusion of molecules between the two layers, we did not observe \[in membranes with $`A_p=(30.0)^2`$\] that pores also lead to the extraction of molecules from the membrane. It is, however, possible that the disassociation of pore-forming membranes occurs on time scales larger than the duration of our simulation. Fast disintegration of the membrane was observed when the projected area was increased to $`A_p=(30.625)^2`$, which was, therefore, the largest projected area set for the membranes in our study.

### III.2 Elasticity and Thermal Undulations

On length scales larger than the membrane thickness, the bilayer can be modeled as a smooth continuous sheet. The thermal undulations of the bilayer can be studied with Helfrich Hamiltonian helfrich  relating the elastic energy to the shape of the membrane:

$$=_S𝑑A\left[\gamma +\frac{1}{2}\kappa \left(c_1+c_22c_0\right)^2+\kappa _Gc_1c_2\right].$$
(8)
The integration in the above equation is carried over the whole surface of the membrane. Three elastic moduli are involved with the Helfrich Hamiltonian: the surface tension $`\gamma `$, the bending modulus $`\kappa `$, and the saddle-splay modulus $`\kappa _G`$. The quantities $`c_1`$ and $`c_2`$ appearing in the above equation are the local principle curvatures of the surface (see a rigorous definition in Ref.fdavid ) which are surface invariants with respect to similarity transformations (translations and rotations), while $`c_0`$ is the spontaneous curvature of the surface. For flat bilayers $`c_0=0`$. It is customary to dispense with the use of the local curvatures in favor of two other (local) invariants: the mean curvature $`H(c_1+c_2)/2`$, and the Gaussian curvature $`Kc_1c_2`$. If one only considers fluctuations which do not change the topology of the membrane, then the total energy associated with the last term in Eq.(8) is a constanthelfrich . We, thus, arrive to the following more simplified form of Eq.(8):

$$=_S𝑑A\left(\gamma +2\kappa H^2\right).$$
(9)
There are various ways to parameterize the surface. One of them is the Monge representation, where the surface is represented by a height function, $`z=h(x,y)`$, above a reference $`xy`$ plane. For a nearly flat surface, i.e., when the derivatives of the height function with respect to $`x`$ and $`y`$ are small – $`h_x,h_y1`$, one obtains the following approximation for Eq.(9):

$$=𝑑x𝑑y\left[\frac{1}{2}\gamma \left(h_x^2+h_y^2\right)+\frac{1}{2}\kappa \left(^2h\right)^2\right].$$
(10)
Note that unlike Eq.(9), the integral in Eq.(10) runs over the reference $`(x,y)`$ surface rather than over the actual surface of the membrane.

Equations (8)–(10) are expected to be valid only on length scales larger than the thickness of the membrane. The undulatory motion on smaller length scales (which we did not investigate in this study) is dominated by the so called “protrusion modes” protrusion . In our simulations the profile of the bilayers was defined by mapping the system with linear size (of the projected area) $`L`$ onto an $`8\times 8`$ grid whose mesh size $`l=L/8`$ is indeed larger than the typical width of the membrane. The local height of the bilayer was then defined as the average of the local heights of the two layers. The latter were evaluated by the mean height of the lipids (whose positions were identified with the coordinates of their mid atoms) belonging to each layer, which were instantaneously located inside the local grid cell. The discretized form of Hamiltonian (10) is

$$=a_0\underset{\stackrel{}{r}}{}\left[\frac{1}{2}\gamma \left(h_x^2+h_y^2\right)+\frac{1}{2}\kappa \left(^2h\right)^2\right],$$
(11)
where summation goes over the discrete grid coordinates, and $`a_0=l^2`$ is the area of the grid cells. In Fourier coordinates we define

$$h\left(\stackrel{}{r}\right)=\frac{l}{L}\underset{\stackrel{}{q}}{}h_\stackrel{}{q}e^{i\stackrel{}{q}\stackrel{}{r}},$$
(12)
and

$$h_\stackrel{}{q}=\frac{l}{L}\underset{\stackrel{}{r}}{}h\left(\stackrel{}{r}\right)e^{i\stackrel{}{q}\stackrel{}{r}},$$
(13)
where the two-dimensional wavevector $`\stackrel{}{q}`$ has $`8^2=64`$ discrete values satisfying $`\{q_x,q_y=2\pi n/L`$, $`n=4,3,\mathrm{},2,3\}`$. In Fourier space the different modes decouple:

$$=\frac{a_0}{2}\underset{\stackrel{}{q}}{}\left[\gamma |\stackrel{}{q}|^2+\kappa |\stackrel{}{q}|^4\right]|h_\stackrel{}{q}|^2,$$
(14)
and, by invoking the equipartition theorem, we find that the mean square amplitude of the mode $`\stackrel{}{q}`$ (the “spectral intensity”)

$$a_0|h_\stackrel{}{q}|^2=\frac{kT}{\left(\gamma |\stackrel{}{q}|^2+\kappa |\stackrel{}{q}|^4\right)}.$$
(15)
The instantaneous amplitudes of the different $`q`$-modes were evaluated using Eq.(13) once in every 100 MC time units, and were averaged over the course of the simulations. To extract the values of $`\gamma `$ and $`\kappa `$ we used the inverse from of Eq.(15)

$$\frac{1}{a_0|h_\stackrel{}{q}|^2}=\frac{\left(\gamma |\stackrel{}{q}|^2+\kappa |\stackrel{}{q}|^4\right)}{kT},$$
(16)
and plotted $`1/a_0|h_\stackrel{}{q}|^2`$ and a function of $`|\stackrel{}{q}|^2`$. The results of the spectral analysis of the undulations for fluid membranes with $`(A_p=28.75)^2`$ (squares) and $`A_p=(29.375)^2`$ (circles) are presented in fig. 8. The error bars represent one standard deviation in the estimates of the averages. The curves depict the best fit of the numerical data to Eq.(16), obtained when $`\gamma `$ and $`\kappa `$ take the following values:

$`\gamma `$ $`=`$ $`1.4\pm 0.2{\displaystyle \frac{kT}{\sigma ^2}}`$
$`\kappa `$ $`=`$ $`54\pm 2kT,`$ (17)
for $`A_p=(28.75)^2`$, and,

$`\gamma `$ $`=`$ $`0.6\pm 0.2{\displaystyle \frac{kT}{\sigma ^2}}`$
$`\kappa `$ $`=`$ $`42\pm 2kT,`$ (18)
for $`A_p=(29.375)^2`$. We verified the validity of Eq.(16) by attempting to fit our data to other polynomial functional forms, including a constants and a $`|\stackrel{}{q}|^6`$ terms. The contributions of these terms to the fit were small, and did not result a significant change in our estimates of $`\gamma `$ and $`\kappa `$, based on Eq.(16).

The above values of the bending modulus $`\kappa `$ are somewhat larger than the values commonly reported in experiments in phospholipids: $`\kappa 1020kT`$ seifert\_lipowsky . We should, in principle, correct our values of $`\kappa `$ by considering the reduction in the effectively measured bending modulus due to long-wavelength thermal undulations. The correction term depends logarithmically on the size of the system safran ; peliti\_leibler :

$$\frac{\mathrm{\Delta }\kappa }{\kappa }=\frac{3kT}{4\pi \kappa }\mathrm{ln}\left(\frac{L}{a}\right),$$
(19)
where $`a`$ is some microscopic length. Setting $`al=L/8`$ (the mesh size in our simulations), and using the values of $`\kappa `$ in Eqs.(17) and (18), we find that this correction amounts to about $`1\%`$ of the value of $`\kappa `$ and, therefore, falls within the uncertainty in our estimates of the bending modulus.

The fact that the surface tension $`\gamma `$ is negative has an interesting implication: It means that the size of the membrane cannot grow indefinitely, but an upper bound exists

$$L_c2\pi \sqrt{\frac{\kappa }{|\gamma |}}$$
(20)
above which there are small $`q`$ modes with $`|\stackrel{}{q}|<q_c=2\pi /L_c`$ that make the system unstable \[see Eq.(14)\]. One can also understand the origin of this instability in “real” space, rather in $`q`$ space: When $`\gamma <0`$ the elastic energy of the membrane decreases by increasing its area, and it is the bending energy that stabilizes the system. Modes with larger wavelength (smaller $`|\stackrel{}{q}|`$) require smaller curvatures to increase the area of the membrane and, thus, cost less bending energy. Membranes with linear size $`L>L_c`$ have long wavelength modes that reduce the elastic area energy more than they increase the elastic bending energy. For systems whose linear size is smaller than, but close to, $`L_c`$ the amplitude of the small $`q`$ modes become large \[see Eq.(15)\], and the approximation in Eqs.(11) and (14) is no longer valid. It is then necessary to include higher order terms in the Hamiltonian (14), and to consider their influence on the spectrum of thermal undulations and on the stability of the membrane.

### III.3 Pore Formation

The projected area $`\stackrel{~}{A}_p`$ at which the surface tension vanishes is called the saturated (Schulman) area schulman ; degennes\_taupin . One can evaluate the saturated area of our model membranes using a linear approximation for the relation between the surface tension and the excess area $`\delta A_p=A_p\stackrel{~}{A}_p`$ remark

$$\gamma =K_A\left(\frac{\delta A_p}{\stackrel{~}{A}_p}\right).$$
(21)
The coefficient $`K_A`$ appearing in the above equations is the area compressibility modulus of the membrane. Using the values of $`\gamma `$ obtained from the simulations for $`A_p=(28.75)^2826.6`$, and $`A_p=(29.375)^2862.9`$ \[Eqs.(17) and (18)\] in Eq.(21), we derive the following estimates:

$$K_A=19.6\pm 6.6\frac{kT}{\sigma ^2},$$
(22)
and,

$$\stackrel{~}{A}_p=890\pm 17,$$
(23)
for $`K_A`$ and $`\stackrel{~}{A}_p`$. Equation (23) suggests that our membrane with $`A_p=(30.0)^2=900`$ might be found above the saturated area $`\stackrel{~}{A}_p`$ and, therefore, may have a positive surface tension. The attempt to verify this conjecture by analyzing the spectrum of the membranes, as done for the fluid membranes with lower projected areas, is hampered by two technical difficulties. The first one is related to the flip-flop motion of lipids between the two leaflets. The trans-bilayer diffusion makes it computationally complicated to determine which lipids are related to which layer of the membrane and, therefore, it becomes difficult to calculate the profile of the layers. The other difficulty results from the holes which are created in the membrane. These pores may have an area larger than $`a_0`$, the area of the grid cells. In such a case we find an empty cell with no lipids inside, and the height of the membrane at the corresponding grid point cannot the evaluated (unless one interpolates this value using the height of the membrane at the adjacent grid points). We have taken advantage of the fact that the typical time for the appearance of the pores and for the flip-flop motion which accompanies their formation, is larger than the relaxation time of the spectrum, and used short MC runs (during which pores were not observed) to estimate the surface tension of the membrane. We found a positive surface tension with a magnitude of the order of $`\gamma 1kT/\sigma ^2`$, which is roughly half an order of magnitude larger than the value anticipated by Eqs.(21)-(23). While the spectral analysis supports our conjecture that $`\gamma `$ is positive for $`A_p=(30.0)^2`$, one should not attempt to use many independent short runs to achieve a more accurate estimate of $`\gamma `$. It is unclear how well equilibrated the membranes in these short MC runs are. Moreover, it is incorrect to base such an estimate on statistical averaging restricted to membranes without pores. The creation of the pores tends to reduce the surface tension since they make the effective area of the membrane smaller.

Membranes with a positive surface tension can reduce their elastic energy by decreasing their area, and the formation of pores is obviously one of the mechanisms to achieve that. Other ways of reducing the membrane area which are not possible in our model is to decrease the projected area or to increase the area density by adsorbing lipids from the solvent. For the case of a pore formation, one has to consider the line tension energy price involved with the creation of the hole. The simplest theoretical model discussing pore formation was suggested by Litster litster . In this zero-temperature model, the contribution of a circular hole of radius $`R_{\mathrm{pore}}`$ to the free energy of a membrane with a positive surface tension $`\gamma `$ is given by

$$F_{\mathrm{pore}}=\gamma \pi R_{\mathrm{pore}}^2+\lambda 2\pi R_{\mathrm{pore}},$$
(24)
where $`\lambda `$ is the line tension of the hole. According to this model a pore with a radius larger than the critical value of $`\lambda /\gamma `$ is predicted to grow without bound. Such a thermodynamically large circular hole can be created only if the critical energy barrier $`\pi \lambda ^2/\gamma `$ is accessible by thermal fluctuations. At a finite temperature it is necessary to take into account the entropy of the pores and the picture becomes more complicated. Recent computer simulations schillcock  have demonstrated that the typical shape of thermally induced pores is non-circular but rather of a self-avoiding ring or a branched polymer. The most striking feature predicted by this study was, however, the fact that pores can appear at zero, and even at small negative surface tension.

The major drawback of the above model is the fact that while it predicts the expansion of the pore without limit, the first term in Eq.(24), assuming a linear relation between the reduction in elastic energy and the area of the pore, applies to small pores only. An improved model can be obtained by assuming other forms of the free energy dependence on the pore area. We first consider a zero temperature model where the membrane does not fluctuate in the normal direction. The free energy of the membrane (which at zero temperature coincides with the potential energy) has a minimum at the saturated area which we shall now denote by $`\stackrel{~}{A}^E`$ \[compare this notation with the one used in Eq.(21)\] to indicate that it is determined by energy consideration. The subscript $`p`$ has been omitted since the projected area is also the total area of the membrane in this case. Close to $`\stackrel{~}{A}^E`$ we can use the quadratic approximation to describe the dependence of the free energy density $`f`$ on the excess area $`\delta A^E=A\stackrel{~}{A}^E`$

$$f\frac{F}{\stackrel{~}{A}^E}=\frac{1}{2}K_A^E\left(\frac{\delta A^E}{\stackrel{~}{A}^E}\right)^2,$$
(25)
where, as in the case of the saturated area, we use the superscript $`E`$ in the notation of the area compressibility $`K_A^E`$. If a pore of area $`A_{\mathrm{pore}}`$ is formed then the area of the membrane is reduced by $`A_{\mathrm{pore}}`$ and, consequently, the pore contribution to the free energy density is given by:

$`f_{\mathrm{pore}}(A,A_{\mathrm{pore}})`$ $`=`$ $`{\displaystyle \frac{1}{2}}K_A^E\left({\displaystyle \frac{\delta A^EA_{\mathrm{pore}}}{\stackrel{~}{A}^E}}\right)^2`$
$``$ $`{\displaystyle \frac{1}{2}}K_A^E\left({\displaystyle \frac{\delta A^E}{\stackrel{~}{A}^E}}\right)^2`$
$`+`$ $`{\displaystyle \frac{2\lambda \sqrt{\pi }}{\stackrel{~}{A}^E}}\sqrt{A_{\mathrm{pore}}}.`$ (26)
As in Eq.(24), we consider a circular hole and, thus, its perimeter and area are related by $`\mathrm{\Gamma }=\sqrt{4\pi A}`$. The equilibrium size of the pore $`A_{\mathrm{pore}}^{}`$ is found by solving the equation $`f_{\mathrm{pore}}/A_{\mathrm{pore}}=0`$, and in addition by verifying that $`f_{\mathrm{pore}}(A,A_{\mathrm{pore}}^{})<f_{\mathrm{pore}}(A,0)=0`$. While in Litster’s model a membrane with positive surface tension can be only metastable against the formation of a pore, the model presented here yields a different scenario: Pores are thermodynamically unfavorable as long as the line tension satisfies

$$\lambda >\lambda ^{}=\sqrt{\frac{2}{27\pi }}\frac{(\delta A^E)^{3/2}}{\stackrel{~}{A}^E}K_A^E.$$
(27)
At this value a first order first transition occurs, and a pore of size $`A_{\mathrm{pore}}^{}=2/3\delta A^E`$ is created. The pore grows gradually as $`\lambda `$ is decreased below this value. When $`\lambda 0`$, $`A_{\mathrm{pore}}^{}\delta A^E`$, and the effective area of membrane attains the optimal (Schulman) value $`\stackrel{~}{A}^E`$. As in Litster’s model, there exists a free energy barrier for the formation of the pore. At the transition ($`\lambda =\lambda ^{}`$) the height of the barrier is

$$\mathrm{\Delta }F\lambda ^{4/3}(\delta A^E)^{1/3}/K_{A}^{E}{}_{}{}^{1/3}.$$
(28)
A theory for the entropic contribution to the free energy of the pore has been recently presented by Sens and Safran sens\_safran . According to this theory, hole formation is one of the mechanisms to “redeem” the degrees of freedom associated with the long wavelength modes in the fluctuation spectrum which are eliminated by the surface tension. To a first approximation, the effect of this entropic surface tension can be easily incorporated into Eqs.(25) and (26). Let us assume for a moment that $`K_A^E=0`$, and that the projected area of the membrane $`A_p`$ is fixed. Because of the thermal fluctuations, the total area of the membrane will be larger than $`A_p`$. As has been explained in Ref. sens\_safran , there exists an optimal total area at which the membrane is tensionless

$$\stackrel{~}{A}^SA_p\left[1+\frac{kT}{8\pi \kappa }\mathrm{ln}\left(\frac{A_p}{l^2}\right)\right],$$
(29)
where $`\kappa `$ is the bending modulus, and $`l`$ is some molecular cut-off length. The superscript $`S`$ denotes the fact that the optimal area discussed here is entropic in nature, and does not need to be equal to $`\stackrel{~}{A}^E`$ in Eq.(25). A membrane with $`A\stackrel{~}{A}^S`$ will experience a surface tension. The free energy associated with this entropic surface tension can be calculated analytically. Here, however, we shall use the quadratic approximation in $`\delta A^SA\stackrel{~}{A}^S`$

$$f\frac{F}{A_p}=\frac{1}{2}K_A^S\left(\frac{\delta A^S}{A_p}\right)^2,$$
(30)
which is valid only close to the minimum of the free energy at $`\stackrel{~}{A}^S`$. The entropic area compressibility in Eq.(30) is given by sens\_safran

$$K_A^S=\frac{32\pi ^3\kappa ^2}{A_pkT}.$$
(31)
Combining the energetic (25) and the entropic (30) contributions to the free energy, we find another quadratic form for the total free energy of the membrane

$$f\left(A\right)\frac{F}{A_p}=\frac{1}{2}K_A\left(\frac{\delta A}{A_p}\right)^2,$$
(32)
where the excess area $`\delta A=A\stackrel{~}{A}`$ is defined with respect to the minimum at

$$\stackrel{~}{A}=\frac{K_A^EA_p+K_A^S\stackrel{~}{A}^S}{\frac{A_p}{\stackrel{~}{A}^E}K_A^E+K_A^S},$$
(33)
and the effective area compressibility is equal to

$$K_A=K_A^E\frac{A_p}{\stackrel{~}{A}^E}+K_A^S.$$
(34)
The optimal area and the area compressibility appearing in the above two equations (and which include both energetic and entropic contributions) should replace their purely energetic counterparts in Eq.(26) for the pore free energy density and in Eqs.(27) and (28) for the critical line tension and the free energy barrier. For typical values of phospholipids: $`\kappa =10kT5\times 10^{13}`$ ergs, and $`A_p=(10\mu m)^2=10^6\mathrm{cm}^2`$, we get upon substitution in Eq.(31), $`K_A^S5\times 10^3\mathrm{ergs}/\mathrm{cm}^2`$. This value of $`K_A^S`$ is several orders of magnitude smaller than the area compressibility typically found in stretching experiments $`K_A10^2\mathrm{ergs}/\mathrm{cm}^2`$ seifert\_lipowsky  and, therefore, the entropic contribution to $`K_A`$ and $`\stackrel{~}{A}`$ can be neglected remark2 . The effect of the thermal fluctuations can be felt only in small membranes (such as in this paper), as $`K_A^S`$ increases by decreasing the projected area. Unfortunately, the membranes in our computer study are too small to allow the description of this effect by Eqs.(29)–(31). The derivation of these equations is based on the assumption that the long wavelength behavior of the membrane is dominated by the surface tension. This requires that $`\kappa (2\pi /\sqrt{A_p})^2\gamma `$ \- a criterion which is not satisfied in our case. The long wavelength fluctuations in our membranes are mainly controlled by the curvature elasticity.

Is the appearance of pores in the simulations in accord with the model described by Eq.(26)? In order to answer this question we need to evaluate the line tension $`\lambda `$ of the pore. The line tension $`\lambda `$ has the dimensions of energy per unit length. Its magnitude can be estimated by noting that the lipids on the rim of the pore have 1-2 less neighbors compared to the other lipids. Therefore, the energy cost associated with each such lipid is of the order of the interaction energy between two adjacent molecules which is roughly $`kT`$. The length occupied by each lipid along the perimeter of the hole is of the order of $`\sigma `$, and so $`\lambda kT/\sigma `$. This value of $`\lambda `$ should be smaller than the critical value $`\lambda ^{}`$ given by Eq.(27). Using the values of $`K_A`$ and $`\stackrel{~}{A}`$ provided by Eqs.(22) and (23), we arrive at the estimate $`\lambda ^{}kT/\sigma `$ for the membrane with $`A_p=30^2=900`$ ($`\delta A10`$). This means that $`\lambda `$ and $`\lambda ^{}`$ are of the same order of magnitude and, thus, may obey the criterion given by Eq.(27) for the thermodynamic stability of membranes with holes. The fact the pores in our simulations appear for only short time intervals, before they close up, may indicate that $`\lambda `$ is, in fact, slightly larger than $`\lambda ^{}`$, and that the pores are only metastable. In addition to to the values of $`\lambda `$ and $`\lambda ^{}`$, we also need to check the free energy barrier for the formation of the pores, as given by Eq.(28). We find $`\mathrm{\Delta }FkT`$, and so the opening of a pore can be nucleated by thermal fluctuations.

## IV Summary and Discussion

We have introduced a new simple computer model for bilayer membranes whose main feature is the fact that the system is simulated in vacuum rather than in aqueous environment. The elimination of the solvent from the simulations greatly improves computational efficiency. Devising a “water-free” model is a great challenge as the water molecules, via the electrostatic interactions between them and the lipids, play a central role in the aggregation and the stabilization of the membrane through the resulting hydrophobic effects. The self-assembly of the system has not been investigated in this paper. (The reader is referred to the simulations presented in Refs. drouffe  and noguchi  in which this issue has been addressed). We did, however, demonstrated that bilayers, once they are formed, can be stable without the surrounding solvent. One only needs to modify the interactions between the lipids, and use effective potentials that compensate for the absence of water by producing a barrier against the disintegration of the membrane. In this model we have been able to stabilize the membrane using pair-wise short range interactions only - another feature that reduces the computational effort involved with the simulations. To the best of our knowledge, this is the first water-free computer model in which fluid membranes are being observed without the need of multibody interactions.

We have found that our simple model reproduces many known features of bilayer membranes, such as the transition from a high density solid phase to a low density fluid phase. We have inspected the spectrum of thermal fluctuations of the fluid membranes, and found it to be well described by the Helfrich Hamiltonian. From the analysis of the spectral intensity of the different modes, we have extracted the surface tension and bending modulus of the system. Based on our numerical results for the surface tension, we have attempted to determine the optimal area of the membrane at which the surface tension vanishes. Indeed, for areas larger than the saturated area, we have found evidences that the surface tension becomes positive. Fluid membranes with positive surface tension can develop pores, and the creation of pores allows the diffusion of lipids from one layer to the other (flip-flops). The opening of holes in our membranes is in agreement with a simple model that takes into account the contributions to the area compressibility of both the inter-particle forces and the thermal fluctuations.

In order to make a closer contact with biological systems, it is necessary to extend the model presented here to include the other elements found in biomemebranes such as the membranes proteins and the cytoskeleton. It would be interesting to see whether these additional components can also be modeled in a coarse-grain manner that would minimize both the computational and the conceptual complexity. Such a model may shed light on an abundance of challenging problems like the effect of the cytoskeleton on the elastic properties of the bilayer, or the role played by the membrane proteins in transport processes across the membrane.

We thank David Andelman, Grace Brannigan, Daan Frenkel, and Thomas Powers for useful discussions, and to Niels Grønbech Jensen and Kurt Kremer for their comments on the manuscript. Special thanks to Philip Pincus and Claus Jeppesen for their substantial contribution to the work, and to Jeffrey Barteet for technical support. This work was supported by the National Science Foundation under Award No. DMR-0203755. The Materials Research Laboratory at UC Santa Barbara is supported by NSF No. DMR-0080034.


text:
# A sum-product estimate in fields of prime order

Let $`q`$ be a prime, $`A`$ be a subset of a finite field $`F:=/q`$, $`|A|<\sqrt{|F|}`$. We prove the estimate $`\mathrm{max}(|A+A|,|AA|)c|A|^{1+\epsilon }`$ for some $`\epsilon >0`$ and $`c>0`$. This extends the result of \[BKT\].

Key words: subsets of finite fields, groups.

MSC: 11B75, 11T30.

§1. Introduction

Let $`q`$ be a prime, $`F=/q`$, $`F^{}=F\{0\}`$, and let $`A`$ be a nonempty subset of $`F`$. We consider the sum set

$$A+A:=\{a+b:a,bA\}$$

and the product set

$$AA:=\{ab:a,bA\}.$$

Let $`|A|`$ denote the cardinality of $`A`$. We have the obvious bounds

$$|A+A|,|AA||A|.$$

The bounds are clearly sharp if $`A=F`$ or if $`|A|=1`$. We can expect some improvement in other cases. However, good lower estimates for $`\mathrm{max}(|A+A|,|AA|)`$ were not known for a long time. Recently a breakthrough was made by J. Bourgain, N. Katz, and T. Tao\[BKT\] who proved the following result.

###### Theorem A

Let $`A`$ be a subset of $`F`$ such that

$$|F|^\delta <|A|<|F|^{1\delta }$$

for some $`\delta >0`$. Then one has a bound of the form

$$\mathrm{max}(|A+A|,|AA|)c(\delta )|A|^{1+\epsilon }$$

for some $`\epsilon =\epsilon (\delta )>0`$ and $`c(\delta )>0`$.

Also, in \[BKP\] a reader can find the history of the problem, generalizations and applications of Theorem A.

However, Theorem A does not estimate $`\mathrm{max}(|A+A|,|AA|)`$ if $`|A|`$ is small comparatively to $`|F|`$. The aim of this paper is to give such an estimate.

###### Theorem 1

Let $`A`$ be a subset of $`F`$ such that

$$|A|<|F|^{1/2}.$$

Then one has a bound of the form

$$\mathrm{max}(|A+A|,|AA|)c|A|^{1+\epsilon }$$

for some $`\epsilon >0`$ and $`c>0`$.

Clearly, Theorem A and Theorem 1 immediately imply uniform estimates for $`|A|<|F|^{1\delta }`$.

###### Corollary 1

Let $`A`$ be a subset of $`F`$ such that

$$|A|<|F|^{1\delta }$$

for some $`\delta >0`$. Then one has a bound of the form

$$\mathrm{max}(|A+A|,|AA|)c(\delta )|A|^{1+\epsilon }$$

for some $`\epsilon =\epsilon (\delta )>0`$ and $`c(\delta )>0`$.

To prove Theorem A, the authors associated with a set $`AF`$ the following set

$$I(A):=\{a_1(a_2a_3)+a_4(a_5a_6):a_1,\mathrm{},a_6A\}.$$

They found lower bounds for $`|I(A)|`$ and applied those bounds for estimation of $`\mathrm{max}(|A+A|,|AA|)`$. Using the main idea of \[BKT\] we give new lower estimates for $`|I(A)|`$.

We denote

$$AA:=\{ab:a,bA\}.$$

Throughout the paper $`c`$ and $`C`$ will denote absolute positive constants. If $`f`$ and $`g`$ are functions, we will write $`fg`$ or $`gf`$ if $`|A|CB`$ for some constant $`C`$ (uniformly with respect to variables in $`f`$ and $`g`$).

###### Theorem 2

Let $`A`$ be a subset of $`F`$ such that

$$|A|<|F|^{1/2}.$$

Then one has a bound of the form

$$|AA|\times |I(A)|c|A|^{5/2}.$$

Observing that $`|I(A)||AA|`$ we deduce from Theorem 2 an estimate for $`|I(A)|`$.

###### Corollary 2

Let $`A`$ be a subset of $`F`$ such that

$$|A|<|F|^{1/2}.$$

Then one has a bound of the form

$$|I(A)|c|A|^{5/4}.$$

Also, one can get a good lower bound for $`|I(A)|`$ if $`|A|>|F|^{1/2}`$.

###### Theorem 3

Let $`A`$ be a subset of $`F`$ such that

$$|A|>|F|^{1/2}.$$

Then one has a bound of the form

$$|I(A)||F|/2.$$

§2. A version of a lemma of J. Bourgain, N. Katz, and T. Tao

and the proof of Theorem 3

Let $`A`$ be a subset of $`F`$ and $`\xi F`$. Denote

$$S_\xi (A):=\{a+b\xi :a,bA\}.$$

We will use the following analog of Lemma 4.2 from \[BKT\].

###### Lemma 1

Let $`\xi F^{}`$ and

$$|S_\xi (A)|<|A|^2.$$

Then

$$|I(A)||S_\xi (A)|.$$

###### Demonstration Proof

By the supposition on $`|S_\xi (A)|`$, the surjection

$$A^2F:(a,b)a+b\xi $$

cannot be one-to-one. Thus there are $`(a_1,b_1)(a_2,b_2)`$ with

$$(a_1a_2)+(b_1b_2)\xi =0.$$
$`1`$
We observe that $`b_1b_2`$. Denote

$$S:=(b_1b_2)S_\xi (A):=\{(b_1b_2)s:sS_\xi (A)\}.$$

We have

$$|S|=|S_\xi (A)|,$$

and every element $`sS`$ can be represented in a form

$$s=(b_1b_2)a+(b_1b_2)b\xi ,a,bA.$$

Substituting $`(b_1b_2)\xi `$ from (1) we get

$$s=(b_1b_2)a+(a_2a_1)b.$$

Therefore, $`SI(A)`$, and Lemma 1 is proved.

To prove Theorem 3 with a slightly weaker form we can use Lemma 2.1 from \[BKP\]. However, we shall need a following generalization of that lemma.

###### Lemma 2

Let $`AF`$, $`GF^{}`$. Then there exists $`\xi G`$ such that

$$|S_\xi (A)||A|^2|G|/(|A|^2+|G|).$$

###### Demonstration Proof

Denote for $`\xi G`$ and $`sF`$

$$f_\xi (s):=\{(a,b):a,bA,a+b\xi =s\}.$$

We have

$$\begin{array}{c}\underset{sF}{}f_\xi (s)^2=|\{(a_1,b_1,a_2,b_2):a_1,b_1,a_2,b_2A,a_1+b_1\xi =a_2+b_2\xi \}|\\ =|A|^2+|\{(a_1,b_1,a_2,b_2):a_1,b_1,a_2,b_2A,a_1a_2,a_1+b_1\xi =a_2+b_2\xi \}|.\end{array}$$

Taking the sum over $`\xi G`$ and observing that for any $`a_1,b_1,a_2,b_2A`$ with $`a_1a_2`$ there is at most one $`\xi `$ such that $`a_1+b_1\xi =a_2+b_2\xi `$ we obtain

$$\underset{gG}{}\underset{sF}{}f_\xi (s)^2|A|^2|G|+|A|^4.$$

Therefore, we can fix $`\xi G`$ so that

$$\underset{sF}{}f_\xi (s)^2|A|^2+|A|^4/|G|.$$
$`2`$
(Clearly we can assume that $`G\mathrm{}`$.) By Cauchy—Schwartz inequality,

$$\left(\underset{sF}{}f_\xi (s)\right)^2|S_\xi (A)|\underset{sF}{}f_\xi (s)^2.$$

Moreover,

$$\underset{sF}{}f_\xi (s)=|A|^2.$$

Therefore, by (2),

$$|S_\xi (A)||A|^4/(|A|^2+|A|^4/|G|)=|A|^2|G|/(|A|^2+|G|).$$

###### Demonstration The proof of Theorem 3

Take $`G=F^{}`$. Taking into account the supposition on $`|A|`$ we deduce from Lemma 2 that for some $`\xi `$

$$|S_\xi (A)||A|^2|G|/(|A|^2+|G|)>|A|^2|G|/(2|A|^2)=(|F|1)/2$$

Therefore, $`|S_\xi (A)||F|/2`$. Also, we have $`|S_\xi (A)|p<|A|^2`$. Thus, Theorem 3 follows from Lemma 1.

§3. Some preparations

Let $`AF^{}`$ and

$$H:=\{sF:|\{(a,b):a,bA,s=a/b\}||A|^2/(5|AA|)\}.$$

Denote by $`G`$ the multiplicative subgroup of $`F^{}`$ generating by $`H`$.

###### Lemma 3

There is a coset $`G_1`$ of $`G`$ such that

$$|AG_1||A|/3.$$
$`3`$
###### Demonstration Proof

Assume the contrary. Let $`A_1,A_2,\mathrm{}`$ be the nonempty intersections of $`A`$ with cosets of $`G`$. Take a minimal $`k`$ so that

$$\left|\underset{i=1}{\overset{k}{}}A_i\right||A|/3$$

and denote

$$A^{}=\underset{i=1}{\overset{k}{}}A_i,A^{\prime \prime }=AA^{}.$$

We have

$$|A^{}||A|/3.$$

On the other hand,

$$|A^{}|\left|\underset{i=1}{\overset{k1}{}}A_i\right|+|A_k|<2|A|/3.$$

Hence,

$$|A|/3|A^{}|2|A|/3$$

and

$$|A^{}|\times |A^{\prime \prime }|=|A^{}|(|A||A^{}|)2|A|^2/9.$$
$`4`$
Denote for $`sF^{}`$

$$f(s):=\{(a,b):aA^{},bA^{\prime \prime },a/b=s\}.$$

Note that if $`aA^{},bA^{\prime \prime },`$ then $`a/bH`$. Therefore, for any $`s`$ we have the inequality $`f(s)<|A|^2/(5|AA|)`$. Hence,

$$\underset{sF^{}}{}f(s)^2\frac{|A|^2}{5|AA|}\underset{sF^{}}{}f(s)=\frac{|A|^2|A^{}|\times |A^{\prime \prime }|}{5|AA|}.$$
$`5`$
Denote for $`sF^{}`$

$$g(s):=\{(a,b):aA^{},bA^{\prime \prime },ab=s\}.$$

By Cauchy—Schwartz inequality,

$$\left(\underset{sF}{}g(s)\right)^2|AA|\underset{sF}{}g(s)^2.$$

Therefore,

$$\underset{sF^{}}{}g(s)^2\left(\underset{sF}{}g(s)\right)^2/|AA|=\frac{(|A^{}|\times |A^{\prime \prime }|)^2}{|AA|}.$$
$`6`$
Now observe that both the sums $`_{sF^{}}f(s)^2`$ and $`_{sF^{}}g(s)^2`$ are equal to the number of solutions of the equation $`a_1^{}a_1^{\prime \prime }=a_2^{}a_2^{\prime \prime }`$, $`a_1^{},a_2^{}A^{}`$, $`a_1^{\prime \prime },a_2^{\prime \prime }A^{\prime \prime }`$. Thus, comparing (5) and (6) we get

$$|A^{}|\times |A^{\prime \prime }||A|^2/5.$$

But the last inequality does not agree with (4), and the proof is complete.

We will use the function $`S_\xi (A)`$ defined in the beginning of §2.

###### Lemma 4

Let $`AF^{}`$ and $`|A|>1`$. Then there exists $`\xi G`$ such that

$$\mathrm{min}(|A|^3/(5|AA|),|A|^2|G|/(|A|^2+|G|))S_\xi (A)<|A|^2.$$
$`7`$
###### Demonstration Proof

We consider two cases.

1. Case 1: there exists $`gG`$ such that $`S_g(A)=|A|^2`$. We claim that

$$\xi G|A|^3/(5|AA|)S_\xi (A)<|A|^2.$$
$`8`$
Assume that (8) does not hold. Take an arbitrary $`gG`$ satisfying $`S_g(A)=|A|^2`$ (this means that the elements $`a+bg`$, $`a,bA`$ are pairwise distinct) and an arbitrary $`hH`$. Denote

$$A_h=\{bA:bhA\}.$$

We have $`|A_h||A|^2/(5|AA|)`$ because $`hH`$. By our supposition on $`g`$, all the sums $`a+b(gh)=a+(bh)g`$, $`aA`$, $`bA_h`$, are distinct. Therefore, $`S_{gh}(A)|A|^3/(5|AA|)`$. Our supposition that (8) does not hold implies that $`S_{gh}(A)=|A|^2`$.

So, we see that if an elements $`gG`$ satisfies the condition $`S_g(A)=|A|^2`$ then for any $`hH`$ the elements $`gh`$ also satisfies this condition. Since $`H`$ generates $`G`$ we deduce that the condition $`S_g(A)=|A|^2`$ holds for all elements $`gG`$. But this is impossible because $`S_1(A)|A|(1+|A|)/2<|A|^2`$, and (8) is proved.

2. Case 2: for all $`gG`$ we have $`S_g(A)<|A|^2`$. Then the existence of a required $`\xi G`$ immediately follows from Lemma 2, and the proof of Lemma 4 is complete.

###### Lemma 5

Let $`G`$ be a subgroup of $`F^{}`$, $`BG`$, $`|B|<\sqrt{|F|}`$. Then

$$|BB||A|^{5/2}/|G|.$$

###### Demonstration Proof

We use arguments from \[HBK\]. Consider the cosets $`G_1,G_2,\mathrm{}`$ of $`G`$ in $`F^{}`$. For any coset $`G_t`$ and $`sG_t`$ denote

$$N_t:=|\{(g_1,g_2):g_1,g_2G,g_1g_2=s\}|.$$

(It is clear that the definition is correct, namely, it does not depend on the choice of $`sG_t`$.) We order $`G_1,G_2,\mathrm{}`$ in such a way that

$$N_1N_2\mathrm{}.$$
$`9`$
Also, we consider $`N_t=0`$ if $`t`$ exceeds the number of the cosets of $`G`$ in $`F^{}`$. Lemma 5 from \[HBK\] claims that if

$$|G|^4T<|F|^3,$$
$`10`$
then

$$\underset{t=1}{\overset{T}{}}N_t(|G|T)^{2/3}.$$
$`11`$
Therefore, if (10) holds, then, by (9) and (11),

$$N_T|G|^{2/3}T^{1/3},$$
$`12`$
and, moreover,

$$\underset{t=1}{\overset{T}{}}N_t^2\underset{t=1}{\overset{T}{}}\left(|G|^{2/3}t^{1/3}\right)^2|G|^{4/3}T^{1/3}.$$
$`13`$
Now denote for $`t=1,2,\mathrm{}`$

$$L_t:=|\{(b_1,b_2):b_1,b_2B,b_1b_2G_t\}|,$$

$$M_t:=|\{(b_1,b_2,b_3,b_4):b_1,b_2,b_3,b_4B,b_1b_2=b_3b_4G_t\}|.$$

We consider $`L_t=M_t=0`$ if $`t`$ exceeds the number of the cosets of $`G`$. For every element $`b_1B`$ there are at most $`N_t`$ elements $`b_2B`$ such that $`b_1b_2G_t`$. Therefore,

$$L_tN_t|B|.$$
$`14`$
Further, for every element $`sG_t`$ there are at most $`N_t`$ elements $`b_1B`$ such that $`b_1sB`$. Therefore,

$$M_tN_tL_tN_t^2|B|.$$
$`15`$
Also,

$$\underset{t}{}L_t=|\{(b_1,b_2):b_1,b_2B,b_1b_2\}|=|B|(|B|1).$$
$`16`$
The statement of Lemma 5 is trivial if $`|G||B|^{3/2}`$. Thus, we assume that

$$|G|<|B|^{3/2}.$$
$`17`$
Let us take

$$T:=[|B|^{3/2}/|G|].$$

By the supposition on $`|B|`$ and (17), we have

$$|G|^4T|G|^3|B|^{3/2}<|B|^6<|F|^3,$$

and (10) holds. By (15) and (13),

$$\underset{t=1}{\overset{T}{}}M_t|B|\underset{t=1}{\overset{T}{}}N_t^2|B||G|^{4/3}T^{1/3}.$$

By (15), (12), and (16),

$$\underset{t>T}{}M_t\underset{t>T}{}N_tL_tN_T\underset{t}{}L_t|G|^{2/3}T^{1/3}|B|^2.$$

Thus, taking into account the definition of $`T`$, we get

$$\underset{t}{}M_t|B||G|^{4/3}T^{1/3}+|B|^2|G|^{2/3}T^{1/3}|B|^{3/2}|G|.$$
$`18`$
Denote for $`sF^{}`$

$$f(s):=|\{(b_1,b_2):b_1,b_2B,b_1b_2=s\}|.$$

By Cauchy—Schwartz inequality,

$$\left(\underset{sF^{}}{}f(s)\right)^2\underset{sF^{}}{}f(s)^2|\{b_1b_2:b_1,b_2G,b_1b_2\}|.$$

Also,

$$\underset{sF^{}}{}f(s)=\underset{t}{}L_t,\underset{sF^{}}{}f(s)^2=\underset{t}{}M_t,$$

and, by (16) and (18), we find

$$|BB|1+(|B|(|B|1))^2/(|B|^{3/2}|G|)|B|^{5/2}/|G|.$$

Lemma 5 is proved.

§4. The proofs of Theorems 1 and 2

###### Demonstration Proof of Theorem 2

In the case $`|AA|>|A|^{3/2}`$ the assertion is trivial because $`|AA||A|`$, $`|I(A)||AA|`$. Thus, we will consider that $`|AA||A|^{3/2}`$ and, therefore,

$$|A|^3/(|AA|)|A|^{3/2}.$$
$`19`$
We take the group $`G`$ defined in the beginning of §3. If $`|G|>|A|^{3/2}`$, then

$$|A|^2|G|/(|A|^2+|G|)|A|^2|A|^{3/2}/(|A|^2+|A|^{3/2})|A|^{3/2}/2.$$

Hence, by Lemma 4 and (19), there exists $`\xi G`$ such that

$$|A|^{3/2}S_\xi (A)<|A|^2.$$

Lemma 1 claims that $`|I(A)||A|^{3/2}`$. Thus,

$$|AA|\times |I(A)||A|\times |A|^{3/2},$$

and we have the required estimate.

Now it suffices to consider the case if both the conditions (19) and

$$|G||A|^{3/2}$$
$`20`$
are satisfied. We see from (20) that

$$|A|^2|G|/(|A|^2+|G|)|A|^2|G|/(2|A|^2)=|G|/2.$$

Hence, by Lemma 4 and (19), there exists $`\xi G`$ such that

$$|G|S_\xi (A)<|A|^2,$$

and, by Lemma 1,

$$|I(A)||G|.$$
$`21`$
We take a coset $`G_1`$ of $`G`$ in accordance with Lemma 3. Fix an arbitrary $`g_1G_1`$. Let

$$B:=\{bG:bg_1A\}.$$

By (3), we have

$$|B|=|AG_1||A|/3.$$

We get from Lemma 5 that

$$|AA||AG_1AG_1|=|BB|(|A|/3)^{5/2}/|G||A|^{5/2}/|G|.$$
$`22`$
Now Theorem 2 follows from (21) and (22).

Note that $`|AA||I(A)|`$. Thus, Theorem 2 implies the inequality

$$|I(A)||A|^{5/4}$$
$`23`$
provided that $`|A|<|F|^{1/2}`$. Theorem 1 is a corollary of (23) and Lemma 2.4 from \[BKT\].

The research was fulfilled during the author’s visit to the University Aroma ARE and was supported by the Mathematical Department of the University and the Institute Nazionale di Alta Matematica.

References

\[BKT\] J. Bourgain, N. Katz, and T. Tao, A sum-product estimate in finite fields and their applications. ArXiv:math.CO/0301343 v1, January 29, 2003.

\[HBK\] D. R. Heath-Brown, S. V. Konyagin, New bounds for Gauss sums derived from $`k`$th powers, and for Heilbronn’s exponential sums. Quart. J. Math. 51 (2000), 221–235.

Department of Mechanics and Mathematics, Moscow State University,

Moscow, 119992, Russia.

E-mail address: konyaginok.ru